# Tweaking Solve for systems of polynomial equations

I am trying to solve the following system of $6$ quadratic equations in $6$ variables:

system =
{ 1/144 (Sqrt[3] x (Sqrt[5] + 1) + 2Sqrt[3] (t1 + t2 - 1) + 4 Sqrt[3] t1
- 2Sqrt[3] t2 + 2 Sqrt[3] (x - 1))^2 + 1/16 (x (Sqrt[5] + 1) - 2t1 - 4t2 - 2x + 4)^2
+ 1/9 (x (Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6]) + 3(x - 1) Sqrt[1/2 Sqrt[5] + 7/6])^2 == s,
1/144 (Sqrt[3] x (Sqrt[5] + 1) + 2Sqrt[3] (x - 1) + 2Sqrt[3] (Sqrt[5] + 1))^2
+ 1/16 (x (Sqrt[5] + 1) - 2x + 2)^2 + 1/9 (x (Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])
+ 3*(x - 1) Sqrt[1/2 Sqrt[5] + 7/6] - Sqrt[3] + 3 Sqrt[1/2 Sqrt[5] + 7/6])^2 == s,
1/36 (Sqrt[3] (t1 + t2 - 1) + 2 Sqrt[3] t1 - Sqrt[3] t2 - Sqrt[3] (Sqrt[5] + 1))^2
+ 1/4 (t1 + 2t2 - 1)^2 + 1/9 (Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])^2 == s,
1/144 (Sqrt[3] (u1 + u2 - 1) (Sqrt[5] + 1) - 2Sqrt[3] (t1 + t2 - 1)
- 4Sqrt[3] t1 + 2 Sqrt[3] t2 + 2Sqrt[3]u1 - 4Sqrt[3] u2)^2 +
1/16 ((u1 + u2 - 1) (Sqrt[5] + 1) - 2t1 - 4t2 - 2u1 + 2)^2 +
1/9 ((u1 + u2 - 1)(Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])
+ 3u1 Sqrt[1/2 Sqrt[5] + 7/6] + 3u2 Sqrt[1/2 Sqrt[5] + 7/6])^2 == s,
1/144 (Sqrt[3] (u1 + u2 - 1) (Sqrt[5] + 1) + Sqrt[3] x (Sqrt[5] + 1)
+ 2Sqrt[3] u1 - 4 Sqrt[3] u2 + 2Sqrt[3] (x - 1))^2 +
1/16 ((u1 + u2 - 1)(Sqrt[5] + 1) - x (Sqrt[5] + 1) - 2u1 + 2x - 2)^2
+ 1/9 ((u1 + u2 - 1)(Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])
+ x (Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6]) + 3u1 Sqrt[1/2 Sqrt[5] + 7/6]
+ 3u2 Sqrt[1/2 Sqrt[5] + 7/6] + 3(x - 1)Sqrt[1/2 Sqrt[5] + 7/6])^2 == s,
1/144 (Sqrt[3] (u1 + u2 - 1)(Sqrt[5] + 1) + 2Sqrt[3] u1
- 4Sqrt[3] u2 - 2Sqrt[3](Sqrt[5] + 1))^2 + 1/16 ((u1 + u2 - 1)(Sqrt[5] + 1) - 2u1)^2
+ 1/9 ((u1 + u2 - 1)(Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])
+ 3u1 Sqrt[1/2 Sqrt[5] + 7/6] + 3u2 Sqrt[1/2 Sqrt[5] + 7/6] +
Sqrt[3] - 3Sqrt[1/2 Sqrt[5] + 7/6])^2 == s};


Unfortunately Mathematica does not give a solution:

Solve[ system, {s, t1, t2, u1, u2, x}]


It ran for about an hour and it didn't finish.
Numerically the system can be solved with NSolve.
I am looking for the solution, where all the variables ({s, t1, t2, u1, u2, x}) are positve.
Numerically it should be this one:

{s -> 1.8156, t1 -> 0.290762, t2 -> 0.352453,
u1 -> 0.332044,  u2 -> 0.0729072, x -> 0.645495}


Is there a way to obtain the exact solutions?
Things I tried that didn't work:

• setting the domain to Reals. (Then even Nsolve doesn't find a solution)
• adding an extra inequality s > 0.
• eliminating the variable s by hand, which gives a system of $5$ quadratic equations with $5$ variables.
• With "Mathematica does not give a solution" you mean you let it run for a couple of hours and it didn't finish or did it return but unevaluated? Commented Jan 15, 2014 at 13:04
• @Sjoerd C. de Vries: I let it run for about an hour and it didn't finish. Commented Jan 15, 2014 at 13:05
• If the coefficients in your quadratic equations were rationals (so you could write equivalent equations with integer coefficients), it would be easier for Mathematica to provide an exact answer. Commented Mar 22, 2014 at 15:41

The issue we encounter here is an instance of utility of update of the Root object in Mathematica 9.

We should provide the assumptions you've mentioned:

assumptions = {s > 0, x > 0, t1 > 0, t2 > 0, u1 > 0, u2 > 0};


Now we can expect that the new supplemented system Join[system, assumptions] might be solved in terms of the root objects involving multivariable polynomials, in fact we can get the solution instantly,
let's solve the system symbolically e.g. with respect to s:

s /. ToRules @ Reduce[ Join[ system, assumptions], {s, t1, t2, u1, u2, x}]

 Root[{-5 + #1^2 &,
103432987 + 46256631 #1 - 609434392 #2 - 272547320 #1 #2 + 1523870292 #2^2
+ 681494812 #1 #2^2 - 2136466980 #2^3 - 955452228 #1 #2^3 + 1858641502 #2^4
+ 831194522 #1 #2^4 - 674274382 #2^5 - 301515206 #1 #2^5 + 54637540 #2^6
+ 24400884 #1 #2^6 - 1318386 #2^7 - 569734 #1 #2^7 + 10082 #2^8 &}, {2, 1}]


We can transform the root object of the form: Root[ {f1, f2}, {k1, k2}] (new capability of Mathematica 9) to the form Root[ f, k] with RootReduce:

RootReduce @ %

Root[  160801 - 1536272 #1 + 21904868 #1^2 - 217962112 #1^3 + 1256858478 #1^4
- 4717349124 #1^5 + 12495147544 #1^6 - 24298372458 #1^7 + 35173457839 #1^8
- 37805732980 #1^9 + 29448527368 #1^10 - 15736320636 #1^11 + 5246771058 #1^12
- 924552262 #1^13 + 60348584 #1^14 - 1318386 #1^15 + 5041 #1^16 &, 5]


Although the full solution is rather involved its numerical values are quite simple:

N @ Reduce[ Join[ system, assumptions], {s, t1, t2, u1, u2, x}]

 s == 1.8156    && t1 == 0.290762 && t2 == 0.352453 &&
u1 == 0.332044 && u2 == 0.0729072 && x == 0.645495


Warning: Reduce[ Join[ system, assumptions], {s, t1, t2, u1, u2, x}] in Mathematica 8 takes rather a long time (it's been running ~ $30$ minutes without solving the system.)

The main improvement in Mathematica 9 is that symbolic solutions can be
written in terms of multipolynomial root objects:

Root[$\{f_1, f_2,\dots \},\{k_1, k_2,\dots \}$ ] represents the last coordinate of the exact vector $\{a_1, a_2,\dots \}$ such that $a_i$ is the $k_i$-th root of the polynomial equation $f_i( a_1,\dots a_{i-1}, x)=0$.

unlike in the former versions of Mathematica. We can't get immediately the solution with Solve in Mathematica 9 since its transformation to radicals might be rather complicated if possible at all and most likly we encounter the latter case. Nontheless we can get from root objects almost as much as we like, see e.g. How do I work with Root objects?. We can guess that in the next versions of the system this issue will be improved.

• I'm going to show a different approach. Regardless, this shows a simple method and is moreover a very nice explanation of the new Root capabilities in version 9. Commented Jan 15, 2014 at 18:20
• @DanielLichtblau Thanks, I was a bit puzzled with the fact that eliminating variables capabilities od Solve or Eliminate hadn't seemed to provide a direct approach to this system in Mathematica 8. Therefore one has to discover somewhat nonstandard solution. Commented Jan 15, 2014 at 18:51
• nice solution! I used Mathematica 8 before, now I use Mathematica 9. Having a Root object is nice. I wonder whether you can express the solution in terms of the finite algebraic extension that is defined by the coefficients of my system?! Commented Jan 16, 2014 at 12:37

An alternative route is to solve numerically to "adequate" precision (can involve trial and error), and then convert to algebraic numbers.

vars = {x, s, u1, u2, t1, t2};In[289]:= Timing[
sols1 = NSolve[system, vars, WorkingPrecision -> 300];]

(* Out[289]= {1.884000, Null} *)


Now select the one of interest.

sols2 = vars /. sols1;
sols3 = Cases[sols2, {_?(Element[#, Reals] && # > 0 &) ..}];
sol = sols3[[1]];


Here it is, to low approximation.

sol // N

(* Out[293]= {0.645495309224, 1.81560223491, 0.332043561796, \
0.0729071548476, 0.290761734788, 0.352452740635} *)


Now convert.

Timing[exvals = RootApproximant[sol]]

(* {5.768000, {Root[-41 + 2406 #1 - 2804 #1^2 - 11920 #1^3 +
35854 #1^4 - 48088 #1^5 + 51560 #1^6 - 92310 #1^7 +
136499 #1^8 - 95276 #1^9 + 46202 #1^10 - 61466 #1^11 +
64290 #1^12 - 19898 #1^13 - 2730 #1^14 + 972 #1^15 + 71 #1^16 &,
5], Root[
160801 - 1536272 #1 + 21904868 #1^2 - 217962112 #1^3 +
1256858478 #1^4 - 4717349124 #1^5 + 12495147544 #1^6 -
24298372458 #1^7 + 35173457839 #1^8 - 37805732980 #1^9 +
29448527368 #1^10 - 15736320636 #1^11 + 5246771058 #1^12 -
924552262 #1^13 + 60348584 #1^14 - 1318386 #1^15 + 5041 #1^16 &,
5], Root[-3231 - 20880 #1 + 32430 #1^2 + 205752 #1^3 +
361976 #1^4 - 31034 #1^5 - 2388194 #1^6 - 3457962 #1^7 -
4152301 #1^8 - 816290 #1^9 + 10299276 #1^10 + 2267938 #1^11 -
3187810 #1^12 + 167604 #1^13 + 382842 #1^14 - 92718 #1^15 +
5751 #1^16 &, 4],
Root[-369 + 9072 #1 + 1053588 #1^2 - 18740400 #1^3 +
51387632 #1^4 - 35043358 #1^5 - 80109508 #1^6 + 205675176 #1^7 -
123966325 #1^8 - 43454770 #1^9 + 50328778 #1^10 -
3388962 #1^11 - 3522208 #1^12 + 1789782 #1^13 - 194904 #1^14 -
46224 #1^15 + 5751 #1^16 &, 5],
Root[80739 - 501966 #1 - 116490 #1^2 + 5939274 #1^3 -
11929078 #1^4 + 92852 #1^5 + 43561778 #1^6 - 101000582 #1^7 +
141911857 #1^8 - 139202168 #1^9 + 102269216 #1^10 -
60817892 #1^11 + 26786870 #1^12 - 7817358 #1^13 +
1642392 #1^14 - 170964 #1^15 + 5751 #1^16 &, 4],
Root[64251 + 243738 #1 - 504834 #1^2 - 2032764 #1^3 - 977942 #1^4 +
2220422 #1^5 + 2767104 #1^6 - 2704036 #1^7 + 2877417 #1^8 +
4747162 #1^9 + 743040 #1^10 - 959712 #1^11 - 901126 #1^12 -
342012 #1^13 - 16650 #1^14 + 24732 #1^15 + 5751 #1^16 &, 5]}} *)


One can verify, more or less, by plugging these into the equations and checking that they give zero to very high precision (much more than what was used to find the roots, say). This is not a proof, unless you precompute a validating precision (which will be impossibly high). An alternative might be to use RootReduce though that could be slow going.