Tweaking Solve for systems of polynomial equations

Question

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.

Answer 1

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.

Answer 2

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.