Solving for the roots of a trilinear system of polynomials

Question

I have been trying to solve for the roots of the following system of trilinear polynomials:

eq1 = a6 x y + a5 x z + a4 y z + a3 x + a2 y + a1 z + a0;
eq2 = x y z + b6 x y + b5 x z + b4 y z + b3 x + b2 y + b1 z + b0;
eq3 = c6 x y + c5 x z + c4 y z + c3 x + c2 y + c1 z + c0;

In the hopes that Mathematica might know the answer, I try:

Solve[eq1 == 0 && eq2 == 0 && eq3 == 0, {x, y, z}]

Mathematica then behaves quite ambiguously. It does not return an empty set to signal that the solution cannot be obtained, but instead starts calculating - which is encouraging. However, the calculation seems to go on indefinitely without ever finishing. Does this mean the solution cannot be obtained? Or should I tweak the equations somehow to help Mathematica find the result quicker? Thanks for any suggestion.

EDIT

Looking through the web, it seems that the U-resultant technique might give the solutions here. Maybe someone could help me implement it in Mathematica? The technique is described here.

Answer 1

One can get solutions by solving in stages. We can eliminate y and z from the system to solve for the x coordinates. Second, we can eliminate z from two of the equations, plug in x, and solve for y. Finally, we can plug in x and y into the third equation and solve for z. The second step introduces extraneous solutions. We get a value for y for each branch lying above a given value for x. (In the OP's system there are always two because the elimination of z results in a quadratic equation.) One improvement would be to determine which branch is correct. It turns out that which of the two solutions for y is correct depends on the values of the parameters. To resolve this is probably what is taking Mathematica so long. At best it seems a matter of hour(s) or more, instead of the seconds and minutes the code below takes.

eq1 = a6 x y + a5 x z + a4 y z + a3 x + a2 y + a1 z + a0;
eq2 = x y z + b6 x y + b5 x z + b4 y z + b3 x + b2 y + b1 z + b0;
eq3 = c6 x y + c5 x z + c4 y z + c3 x + c2 y + c1 z + c0;

eqsys = eq1 == 0 && eq2 == 0 && eq3 == 0;

The first step, solve for x:

eqx = Eliminate[eqsys, {(*x,*)y, z}]    
solx = Solve[eqx, x];
{Length[solx], LeafCount[solx]}
(*  {6, 91519 } *)

Given the size so far, one has to wonder about the value of solving the whole system. In for a penny....

The second step, solve the first two equations for y using the solutions for x:

eqxy = Eliminate[eq1 == 0 && eq2 == 0, {(*x,y,*)z}];
solxy = Flatten[Function[xsol, Join[xsol, #] & /@ Solve[eqxy /. xsol, y]] /@ solx, 1];

solxy /. _Rule -> Rule
(*
  {{Rule, Rule}, {Rule, Rule}, {Rule, Rule}, {Rule, Rule},
   {Rule, Rule}, {Rule, Rule}, {Rule, Rule}, {Rule, Rule},
   {Rule, Rule}, {Rule, Rule}, {Rule, Rule}, {Rule, Rule}}
*)

GatherBy[solxy, First] // Dimensions
(*  {6, 2, 2}  *)

We can see we have six pairs of solutions, two for each solution for x found above. The LeafCount of solxy is up to nearly 6.5 million, which no person can read, so we'll have to be satisfied with this check.

Now, the third and final step, solve the third equation for z. This is relatively easy, since the equation is linear in z, and the LeafCount of the final solution stays well under 26 million.

solxyz = Flatten[
   Function[xysol, Join[xysol, #] & /@ Solve[eq3 == 0 /. xysol, z]] /@ solxy,
   1];

solxyz /. _Rule -> Rule
(*
  {{Rule, Rule, Rule}, {Rule, Rule, Rule}, {Rule, Rule, Rule},
   {Rule, Rule, Rule}, {Rule, Rule, Rule}, {Rule, Rule, Rule},
   {Rule, Rule, Rule}, {Rule, Rule, Rule}, {Rule, Rule, Rule},
   {Rule, Rule, Rule}, {Rule, Rule, Rule}, {Rule, Rule, Rule}}
*)

We still have twelve solutions in six pairs. Below we can check them numerically. Out of each pair, one of the answers solves the original system. We can see in the three examples below that if the parameters are changed, which of the pairs solves the equation changes. It is not clear if all 2^6 = 64 combinations are possible; if they are, one should be able to check in a finite amount of time.

params = DeleteCases[Variables[{eq1, eq2, eq3}], x | y | z];

Clear[check];
(* check[n] returns the differences between the two sides of the equations in the system
    for a random seed n.  A true solution will result in an element of the form
    0. && 0. && 0.  The numerical complexity of the solutions requires high precision
    for accurate computation.  *)
check[n_] := check[n] = BlockRandom[SeedRandom[n];
    Partition[
     eqsys /. Equal -> Subtract /. (solxyz /. #) /. # &[
      Thread[params -> RandomReal[{-10, 10}, Length@params, WorkingPrecision -> 20]]],
     2]];

check[1] /. x_?NumericQ :> x == 0
check[2] /. x_?NumericQ :> x == 0
check[3] /. x_?NumericQ :> x == 0
(*
  {{True, False}, {False, True}, {False, True},
   {False, True}, {False, True}, {False, True}}

  {{True, False}, {True, False}, {True, False},
   {True, False}, {False, True}, {False, True}}

  {{True, False}, {True, False}, {True, False},
   {False, True}, {True, False}, {True, False}}
*)