Torus-geometry algebraic equations using NSolve and Reduce

Question

Somehow a set of naive-looking equations cannot be solved by using NSolve. Mathematica returns a message like this:

NSolve::nsmet: This system cannot be solved with the methods available to NSolve. >>

Try using Reduce, Mathematica still cannot return an answer after some lone time, but it prompts a message after aborting the calculation:

Still waiting for a safe time to evaluate $Inspector[]

They are a set of six Cosine/Sine equations with 6 to-be-solved variables. The six equations are named below as F1, F2, F3, G1, G2, G3. And the six variables are angles a1, b1, b2, c2, a3, c3, all are Reals, and all are in the range of 0 to 2 Pi.

Here are the equations and the code:

F1[a1_, b1_, b2_, c2_] = -(41/4) + 
   Cos[a1] (10 + 3 Cos[b1]) - Cos[b2] (4 + (7 Cos[c2])/2) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
F2[a1_, b1_, b2_, c2_] = (10 + 3 Cos[b1]) Sin[a1] + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
F3[a1_, b1_, b2_, c2_] = 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   3 Sin[b1] - (4 + (7 Cos[c2])/2) Sin[b2] + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
G1[a1_, b1_, a3_, c3_] = 
  Cos[a1] (10 + 3 Cos[b1]) - Cos[a3] (33/2 + (47 Cos[c3])/8) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);
G2[a1_, b1_, a3_, c3_] = (10 + 3 Cos[b1]) Sin[a1] - (33/2 + (47 Cos[c3])/8) Sin[
     a3] + 1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);
G3[a1_, b1_, a3_, c3_] = 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 3 Sin[b1] + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);

    NSolve[F1[a1, b1, b2, c2] == 0 && 
  F2[a1, b1, b2, c2] == 0 && F3[a1, b1, b2, c2] == 0 &&
   G1[a1, b1, a3, c3] == 0 && 
  G2[a1, b1, a3, c3] == 0 && G3[a1, b1, a3, c3] == 0 &&
   0 <= a1 < 2 Pi && 0 <= b1 < 2 Pi && 0 <= b2 < 2 Pi && 
  0 <= c2 < 2 Pi && 0 <= a3 < 2 Pi && 0 <= c3 < 2 Pi, {a1, b1, b2, c2,
   a3, c3}, Reals]

Reduce[F1[a1, b1, b2, c2] == 0 && 
  F2[a1, b1, b2, c2] == 0 && F3[a1, b1, b2, c2] == 0 &&
   G1[a1, b1, a3, c3] == 0 && 
  G2[a1, b1, a3, c3] == 0 && G3[a1, b1, a3, c3] == 0 &&
   0 <= a1 < 2 Pi && 0 <= b1 < 2 Pi && 0 <= b2 < 2 Pi && 
  0 <= c2 < 2 Pi && 0 <= a3 < 2 Pi && 0 <= c3 < 2 Pi, {a1, b1, b2, c2,
   a3, c3}, Reals]

My challenge: Is there any resolution why this set of equations cannot be evaluated/solved numerically by NSolve? And the Reduce? How can the equations be solved numerically?

For approximation, it is totally fine to have only about a 5 digit precision. The numerical solution needs not be very precise or exact.

Is there a built-in Weierstrass substitution (Tangent half-angle substitution), or should this just be substituted manually? — mystery, Nov 11 '15 at 14:58
I show one approach below. As for the challenge as stated, I will observe that it borders on obvious. The system is far from trivial algebraically, and were it but somewhat more complicated (being intentionally vague there), I would expect the minute or so of computational time could become an hour or a day or more. Conciseness of formulation does not equate to ease of solving. — Daniel Lichtblau, Nov 11 '15 at 16:20
New Update: There seem to have bugs in the two answers below. Both cases the Math Kernal closes and jump out of the session. — mystery, Nov 12 '15 at 02:03
New Update: we may get the code working if we try to modify the method via Method->"EndomorphismMatrix" in NSolve. — mystery, Nov 12 '15 at 18:06

Daniel Lichtblau · Answer 1 · 2015-11-12T23:35:23.120

This should get you started. First the basic definitions to keep this self contained.

F1[a1_, b1_, b2_, c2_] = -(41/4) + Cos[a1] (10 + 3 Cos[b1]) - 
   Cos[b2] (4 + (7 Cos[c2])/2) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
F2[a1_, b1_, b2_, c2_] = (10 + 3 Cos[b1]) Sin[a1] + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
F3[a1_, b1_, b2_, c2_] = 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   3 Sin[b1] - (4 + (7 Cos[c2])/2) Sin[b2] + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);
G1[a1_, b1_, a3_, c3_] = 
  Cos[a1] (10 + 3 Cos[b1]) - Cos[a3] (33/2 + (47 Cos[c3])/8) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);
G2[a1_, b1_, a3_, 
   c3_] = (10 + 3 Cos[b1]) Sin[a1] - (33/2 + (47 Cos[c3])/8) Sin[
     a3] + 1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);
G3[a1_, b1_, a3_, c3_] = 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 3 Sin[b1] + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 
132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);

I convert from trigonometric to ordinary polynomials in the "sin,cos" variables. Also add the usual trig identities.

tpolys = TrigExpand[{F1[a1, b1, b2, c2], F2[a1, b1, b2, c2], 
    F3[a1, b1, b2, c2], G1[a1, b1, a3, c3], G2[a1, b1, a3, c3], 
    G3[a1, b1, a3, c3]}];

subs = {Cos[a_] :> ToExpression[StringJoin["c", ToString[a]]], 
   Sin[a_] :> ToExpression[StringJoin["s", ToString[a]]]};

polys = Numerator[Together[tpolys /. subs]];
vars = Variables[polys];
cvars = Cases[vars, ca_ /; StringMatchQ[ToString[ca], "c" ~~ __]];
svars = Map[ToExpression[StringReplacePart[ToString[#], "s", 1]] &, 
   cvars];
extrapolys = cvars^2 + svars^2 - 1;

system = Join[polys, extrapolys]

(*Out[457]= {-123 + 120 ca1 + 36 ca1 cb1 - 48 cb2 - 42 cb2 cc2 - 
  40 sa1 - 12 cb1 sa1 - 12 sb1 + 16 sb2 + 14 cc2 sb2 + 14 sc2, 
 41 - 40 ca1 - 12 ca1 cb1 + 16 cb2 + 14 cb2 cc2 + 120 sa1 + 
  36 cb1 sa1 - 12 sb1 + 16 sb2 + 14 cc2 sb2 + 14 sc2, 
 41 - 40 ca1 - 12 ca1 cb1 + 16 cb2 + 14 cb2 cc2 - 40 sa1 - 
  12 cb1 sa1 + 36 sb1 - 48 sb2 - 42 cc2 sb2 + 14 sc2, 
 240 ca1 + 72 ca1 cb1 - 396 cb2 - 141 cb2 cc2 - 80 sa1 - 24 cb1 sa1 - 
  24 sb1 + 132 sb2 + 47 cc2 sb2 + 47 sc2, -80 ca1 + 132 ca3 - 
  24 ca1 cb1 + 47 ca3 cc3 + 240 sa1 + 72 cb1 sa1 - 396 sa3 - 
  141 cc3 sa3 - 24 sb1 + 47 sc3, -80 ca1 + 132 ca3 - 24 ca1 cb1 + 
  47 ca3 cc3 - 80 sa1 - 24 cb1 sa1 + 132 sa3 + 47 cc3 sa3 + 72 sb1 + 
  47 sc3, -1 + ca1^2 + sa1^2, -1 + ca3^2 + sa3^2, -1 + cb1^2 + 
  sb1^2, -1 + cb2^2 + sb2^2, -1 + cc2^2 + sc2^2, -1 + cc3^2 + sc3^2} *)

Now we can solve as an explicitly polynomial system and discard solutions we know we don't want.

AbsoluteTiming[solns = NSolve[system, vars];]
realsolns = Select[solns, FreeQ[#, Complex] &]

(* Out[458]= {70.711000, Null}

Out[459]= {{ca1 -> 0.990630862403695, ca3 -> 0.9979298911476929, 
  cb1 -> 0.934822592828836, cb2 -> 0.7084047167567334, 
  cc2 -> -0.8662454504796822, cc3 -> -0.9999646433353357, 
  sa1 -> 0.136566813442575, sa3 -> 0.064311224360259, 
  sb1 -> 0.3551150814680469, sb2 -> -0.7058064588857618, 
  sc2 -> -0.4996186746971847, 
  sc3 -> -0.008408833303751029}, {ca1 -> 0.972091493240758, 
  ca3 -> 0.998009033155138, cb1 -> -0.21452025398996, 
  cb2 -> 0.8166902165390839, cc2 -> -0.7788950631571803, 
  cc3 -> -0.8246152396274358, sa1 -> -0.2346020638556038, 
  sa3 -> 0.0630711585948249, sb1 -> -0.9767195392697885, 
  sb2 -> -0.5770763310610337, sc2 -> 0.6271542722956184, 
  sc3 -> 0.5656940010378115}, {ca1 -> 0.971507343886548, 
  ca3 -> 0.9821460285499906, cb1 -> 0.9547232610884067, 
  cb2 -> 0.9261788653705688, cc2 -> 0.4910652194470647, 
  cc3 -> -0.8573410874746823, sa1 -> -0.2370094527740707, 
  sa3 -> -0.1881201193479485, sb1 -> -0.2974953766665983, 
  sb2 -> 0.3770844875111226, sc2 -> 0.8711228116865195, 
  sc3 -> 0.5147487333197118}, {ca1 -> 0.9934049728071653, 
  ca3 -> 0.9999696133112789, cb1 -> 0.8608241235336918, 
  cb2 -> 0.9946210910361825, cc2 -> -0.9110950988182576, 
  cc3 -> -0.9739805097467558, sa1 -> 0.1146584399772491, 
  sa3 -> -0.007795696710271489, sb1 -> 0.5089025424471438, 
  sb2 -> 0.1035803316990695, sc2 -> -0.4121961898790074, 
  sc3 -> -0.2266317667193753}, {ca1 -> 0.990630862243453, 
  ca3 -> 0.9991361112906542, cb1 -> 0.9348225888440881, 
  cb2 -> 0.7084047175344481, cc2 -> -0.8662454499292845, 
  cc3 -> -0.009747051357440676, sa1 -> 0.1365668118363446, 
  sa3 -> 0.04155756756853254, sb1 -> 0.3551150765172533, 
  sb2 -> -0.7058064571601258, sc2 -> -0.4996186701131406, 
  sc3 -> -0.9999524921381223}, {ca1 -> 0.9934049733388737, 
  ca3 -> 0.9999846585487638, cb1 -> 0.8608241316947258, 
  cb2 -> 0.9946210912529978, cc2 -> -0.9110951047761441, 
  cc3 -> -0.2266528053669567, sa1 -> 0.1146584446193901, 
  sa3 -> -0.005539205115196547, sb1 -> 0.5089025584468246, 
  sb2 -> 0.1035803317166068, sc2 -> -0.4121962042393866, 
  sc3 -> -0.9739756187235958}, {ca1 -> 0.9715073439149983, 
  ca3 -> 0.9938480439189615, cb1 -> 0.954723195368781, 
  cb2 -> 0.9261788594717576, cc2 -> 0.4910651488053675, 
  cc3 -> 0.5056808891042688, sa1 -> -0.2370094528804785, 
  sa3 -> -0.1107522660767406, sb1 -> -0.2974954149639326, 
  sb2 -> 0.3770844348460939, sc2 -> 0.871122798741875, 
  sc3 -> -0.8627206024468279}, {ca1 -> 0.972091493240887, 
  ca3 -> 0.9993117833812833, cb1 -> -0.2145200534910795, 
  cb2 -> 0.8166902361952415, cc2 -> -0.7788948435769054, 
  cc3 -> 0.5647276928243201, sa1 -> -0.23460206466597, 
  sa3 -> 0.03709387881898878, sb1 -> -0.9767194217775088, 
  sb2 -> -0.5770761644367379, sc2 -> 0.6271543154910285, 
  sc3 -> -0.8252773053512772}} *)

You can get candidate angles from this by taking inverse trigs, say arccos's of the cosine values.

ArcCos[cvars /. realsolns]

(* Out[462]= {{0.1369949224052047, 0.06435562868860302, 
  0.363037167803944, 0.7835609165398749, 2.61843413925789, 
  3.133183507020114}, {0.2368091879519997, 0.06311303912836261, 
  1.786996949192178, 0.6151442423709425, 2.463698385131359, 
  2.54031803280406}, {0.2392864347266495, 0.1892477383549457, 
  0.3020681676620585, 0.3866463759054244, 1.057484182694673, 
  2.600878105372413}, {0.1149111746751269, 0.007795747894100299, 
  0.5339094523934118, 0.1037664493253039, 2.716729378155218, 
  2.912974564098388}, {0.1369949235785647, 0.04156953504571647, 
  0.3630371790249482, 0.783560915437994, 2.618434138156255, 
  1.580543532495388}, {0.1149111700378031, 0.005539222050615004, 
  0.5339094363568763, 0.1037664472320943, 2.716729392609224, 
  1.799435996760292}, {0.239286434606611, 0.110979946027804, 
  0.3020683885717234, 0.3866463915486326, 1.057484263787367, 
  1.040625309080896}, {0.2368091879514502, 0.03710244125552175, 
  1.78699674391434, 0.6151442083093166, 2.463698035009807, 
  0.970693068874044}} *)

These could in principle require further adjustment by multiples of 2Pi but it appears all these actually meet the required ranges from the original formulation. So we might be done at this point.

There is a caveat to how I selected "real" solutions. In general there could be smallish complex parts that should just be smacked with Chop. This is more prone to happening in version 10 (I used 9 above), due to introduction of a different default method for NSolve. So you might need to alter the selection criterion.

--- edit ---

I may have had a cut-and-paste or other error yesterday. When I repeated this in version 10 I get the results below. Notice I also use the nondefault Method setting for NSolve.

tpolys = TrigExpand[{F1[a1, b1, b2, c2], F2[a1, b1, b2, c2], 
    F3[a1, b1, b2, c2], G1[a1, b1, a3, c3], G2[a1, b1, a3, c3], 
    G3[a1, b1, a3, c3]}];

subs = {Cos[a_] :> ToExpression[StringJoin["c", ToString[a]]], 
   Sin[a_] :> ToExpression[StringJoin["s", ToString[a]]]};

polys = Numerator[Together[tpolys /. subs]];
vars = Variables[polys];
cvars = Cases[vars, ca_ /; StringMatchQ[ToString[ca], "c" ~~ __]];
svars = Map[ToExpression[StringReplacePart[ToString[#], "s", 1]] &, 
   cvars];
extrapolys = cvars^2 + svars^2 - 1;
system = Join[polys, extrapolys]

(* Out[177]= {-123 + 120 ca1 + 36 ca1 cb1 - 48 cb2 - 42 cb2 cc2 - 
  40 sa1 - 12 cb1 sa1 - 12 sb1 + 16 sb2 + 14 cc2 sb2 + 14 sc2, 
 41 - 40 ca1 - 12 ca1 cb1 + 16 cb2 + 14 cb2 cc2 + 120 sa1 + 
  36 cb1 sa1 - 12 sb1 + 16 sb2 + 14 cc2 sb2 + 14 sc2, 
 41 - 40 ca1 - 12 ca1 cb1 + 16 cb2 + 14 cb2 cc2 - 40 sa1 - 
  12 cb1 sa1 + 36 sb1 - 48 sb2 - 42 cc2 sb2 + 14 sc2, 
 240 ca1 - 396 ca3 + 72 ca1 cb1 - 141 ca3 cc3 - 80 sa1 - 24 cb1 sa1 + 
  132 sa3 + 47 cc3 sa3 - 24 sb1 + 47 sc3, -80 ca1 + 132 ca3 - 
  24 ca1 cb1 + 47 ca3 cc3 + 240 sa1 + 72 cb1 sa1 - 396 sa3 - 
  141 cc3 sa3 - 24 sb1 + 47 sc3, -80 ca1 + 132 ca3 - 24 ca1 cb1 + 
  47 ca3 cc3 - 80 sa1 - 24 cb1 sa1 + 132 sa3 + 47 cc3 sa3 + 72 sb1 + 
  47 sc3, -1 + ca1^2 + sa1^2, -1 + ca3^2 + sa3^2, -1 + cb1^2 + 
  sb1^2, -1 + cb2^2 + sb2^2, -1 + cc2^2 + sc2^2, -1 + cc3^2 + sc3^2} *)

AbsoluteTiming[
 solns = NSolve[system, vars, Method -> "EndomorphismMatrix"];]
realsolns = Select[solns, FreeQ[#, Complex] &]

(* Out[178]= {17.763143, Null}

Out[179]= {{ca1 -> 0.979750003005, ca3 -> 0.996336766907, 
  cb1 -> -0.429137780626, cb2 -> 0.0317189516497, 
  cc2 -> -0.866934864031, cc3 -> -0.887281421052, 
  sa1 -> -0.200224702778, sa3 -> 0.0855163545761, 
  sb1 -> -0.903239040968, sb2 -> -0.999496827673, 
  sc2 -> 0.498421449716, 
  sc3 -> 0.461228446452}, {ca1 -> 0.993540741044, 
  ca3 -> 0.999856850531, cb1 -> 0.844669056691, cb2 -> 0.97345748387, 
  cc2 -> -0.913706877747, cc3 -> -0.961917766587, 
  sa1 -> 0.113475970725, sa3 -> -0.0169197646192, 
  sb1 -> 0.535288880379, sb2 -> 0.2288679264, sc2 -> -0.406373894372, 
  sc3 -> -0.273339002747}, {ca1 -> 0.9768581776, ca3 -> 0.90728502041,
   cb1 -> 0.746082752419, cb2 -> 0.683597785102, 
  cc2 -> 0.663822633362, cc3 -> -0.940421761982, 
  sa1 -> -0.213888056756, sa3 -> -0.420516220492, 
  sb1 -> 0.665853231927, sb2 -> 0.72985893691, sc2 -> 0.74789003939, 
  sc3 -> -0.340010160984}, {ca1 -> 0.959065483504, 
  ca3 -> 0.909334697489, cb1 -> -0.631858449557, cb2 -> -0.7184612274,
   cc2 -> 0.75499757274, cc3 -> -0.918342134509, 
  sa1 -> 0.283184389373, sa3 -> 0.416065389022, 
  sb1 -> -0.775083801765, sb2 -> -0.695567009487, 
  sc2 -> -0.655727584595, sc3 -> 0.395787473242}} *)

ArcCos[cvars /. realsolns]

(* Out[180]= {{0.201587262496, 0.0856209299892, 2.01433430099, 
  1.53907205404, 2.61981567378, 2.66221344643}, {0.113720927089, 
  0.0169205723658, 0.564849713432, 0.230914581164, 2.72311067655, 
  2.86473013207}, {0.21555340255, 0.434014217947, 0.728636850982, 
  0.81811557438, 0.844877861826, 2.79466495135}, {0.28711279366, 
  0.429114103112, 2.25474494198, 2.37238385409, 0.715145908928, 
  2.7346674664}} *)

Hoping this is an improvement over what I showed earlier.

--- end edit ---

@ Daniel: Very nice convenient changing of variables, I just wonder whether the "NSolve[system, vars]," you have there, should it be modified to NSolve[system={0,0,0,.....}, vars]? There should be a series of zeros are the RIGHT hand side of equations? Thank you Daniel. — mystery, Nov 11 '15 at 18:46
NSolve is forgiving in terms of whether it receives a set of explicit equations vs. a set (List) of expressions explicitly set to zero (that is, {expr1,...,exprn}==0) vs just a set of expressions. In that last case NSolve infers that they are to be set to zero. — Daniel Lichtblau, Nov 11 '15 at 19:10
Somehow there are some typos, the correct input should be: TrigExpand[{F1[a1, b1, b2, c2], F2[a1, b1, b2, c2], F3[a1, b1, b2, c2], G1[a1, b1, b2, c2], G2[a1, b1, a3, c3], G3[a1, b1, a3, c3]}] — mystery, Nov 11 '15 at 19:44
And your calculation output seem not to be matched to the input... ? Not sure why. — mystery, Nov 11 '15 at 19:44
I may have made a cut-and-paste mistake somewhere. If I get a chance I'll try to replicate the computation. — Daniel Lichtblau, Nov 11 '15 at 20:16
If you run the new-input I edited for you above, I find, at least for my Mathematica 10, it will jump out of the evaluation and turn off the Kernal. — mystery, Nov 11 '15 at 21:00
I later also had that experience, in version 10. I filed a bug report and am even trying to learn a bit more (it's not actually mine) in a debugger. Suffice it to say I believe the crash will be taken seriously. Sorry about that. Aso for getting back v9 behavior, I will suggest in NSolve setting Method->"CompanionMatrix" and see if matters improve. — Daniel Lichtblau, Nov 12 '15 at 00:48

mystery · Answer 2 · 2015-11-12T01:53:08.360

I post an answer based on the comment from J.M. using the Weierstrass substitution and parallel to Daniel Lichtblau set-up.

Pre-set-up:

    F1[a1_, b1_, b2_, c2_] := -(41/4) + Cos[a1] (10 + 3 Cos[b1]) - 
   Cos[b2] (4 + (7 Cos[c2])/2) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);

F2[a1_, b1_, b2_, c2_] := (10 + 3 Cos[b1]) Sin[a1] + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);

F3[a1_, b1_, b2_, c2_] := 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   3 Sin[b1] - (4 + (7 Cos[c2])/2) Sin[b2] + 
   1/16 (-159 + 16 Cos[b2] + 14 Cos[b2] Cos[c2] + 16 Sin[b2] + 
      14 Cos[c2] Sin[b2] + 14 Sin[c2]);

G1[a1_, b1_, a3_, c3_] := 
  Cos[a1] (10 + 3 Cos[b1]) - Cos[a3] (33/2 + (47 Cos[c3])/8) + 
   1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);

G2[a1_, b1_, a3_, 
   c3_] := (10 + 3 Cos[b1]) Sin[a1] - (33/2 + (47 Cos[c3])/8) Sin[
     a3] + 1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);

G3[a1_, b1_, a3_, c3_] := 
  1/4 (50 - 10 Cos[a1] - 3 Cos[a1] Cos[b1] - 10 Sin[a1] - 
      3 Cos[b1] Sin[a1] - 3 Sin[b1]) + 3 Sin[b1] + 
   1/32 (-400 + 132 Cos[a3] + 47 Cos[a3] Cos[c3] + 132 Sin[a3] + 
      47 Cos[c3] Sin[a3] + 47 Sin[c3]);

Code:

tpolys = TrigExpand[{F1[a1, b1, b2, c2], F2[a1, b1, b2, c2], 
    F3[a1, b1, b2, c2], G1[a1, b1, a3, c3], G2[a1, b1, a3, c3], 
    G3[a1, b1, a3, c3]}];

subs = {Cos[a_] :> (
    1 - ToExpression[StringJoin["t", ToString[a]]]^2)/(
    1 + ToExpression[StringJoin["t", ToString[a]]]^2), 
   Sin[a_] :> (2 *ToExpression[StringJoin["t", ToString[a]]])/(
    1 + ToExpression[StringJoin["t", ToString[a]]]^2)};

polys = Numerator[Together[tpolys /. subs]];

vars = Variables[polys];



NSolve[polys, vars,Reals]
Reduce[polys == {0, 0, 0, 0, 0, 0}, vars,Reals]

However, if you really run the code at the last line, I find, at least for my Mathematica 10, it will jump out of the evaluation and turn off the Kernal. This can be a bug, but quotient out this bug problem --- it can be a way to NSolve the problem if the Mathematica bug is fixed!

(Tempted to downvote for use of "kernal"...). I'll try this when I get a chance. Are you saying that the Reduce[...] formulation crashes the Mathematica kernel? If so that would be unrelated to the bug we both ran against in version 10 NSolve, and it would also warrant further investigation. — Daniel Lichtblau, Nov 12 '15 at 16:29
New Update: I get the code work if I tried to modify the method via Method->"EndomorphismMatrix." Namely, try: NSolve[system,vars,Reals,{Method->"EndomorphismMatrix"}] — mystery, Nov 12 '15 at 18:05
How long does it take for that to run? I'm getting no result after several minutes (better than a crash though). — Daniel Lichtblau, Nov 12 '15 at 20:23

Torus-geometry algebraic equations using NSolve and Reduce

2 Answers2

Linked