# Solving for the roots of a trilinear system of polynomials

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.

• You might consider feeding your system to GroebnerBasis[] as a preprocessing step before using Solve[]. – J. M. will be back soon Aug 10 '15 at 1:36
• GroebnerBasis doesn't seem to be too quick on spitting out a result either... – LLlAMnYP Aug 10 '15 at 12:21
• Yes, I read somewhere before that what Solve[] does first of all, is call GroebnerBasis[]. I guess that is the time consuming step here. 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: isc.tamu.edu/publications-reports/tr/9602.pdf – Kagaratsch Aug 10 '15 at 15:42
• Regarding the U-resultant, why don't you give it a go yourself, then show us how far you got, and we'll go from there? – MarcoB Aug 10 '15 at 16:28
• Sure thing! That is what I am going to do. – Kagaratsch Aug 10 '15 at 17:27

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}}
*)

• So does this mean that you found analytic results for all the solutions and you are checking them numerically? Or are the solutions to the last step only available numerically? – Kagaratsch Nov 25 '15 at 19:23
• What if you do GroebnerBasis[eqsys, x, {y, z}] instead, and feed the result of that to Solve[]? – J. M. will be back soon Nov 25 '15 at 19:23
• @J.M. The GroebnerBasis is still running....8GB and counting...I think I'll kill it. Shouldn't I just get eqx, though, more or less? – Michael E2 Nov 25 '15 at 20:53
• @Kagaratsch Yes, I found symbolic solutions, and the last step was to check them numerically. As you can see, for a given set of values for the 21 parameters, one of each pair is a solution (or appears to be, on numerical grounds and consistent with theory). – Michael E2 Nov 25 '15 at 21:00
• Michael, yes, I had thought you could get eqx/solx more directly with that method; I guess not. – J. M. will be back soon Nov 25 '15 at 21:13