I am trying to solve the following polynomial equations in Mathematica:

e= {-1 + c[1]^2 + s[1]^2, -1 + c[2]^2 + s[2]^2, -1 + c[3]^2 + s[3]^2, -1 +     
  (-0.70873 c[1] + c[2] - 0.70548 s[1]) y1[ 1], -1 + (-0.916596 c[2] + c[3] +    
  0.399814 s[2]) y1[ 2], -1 + (c[1] + 0.808085 c[3] + 0.589066 s[3]) y1[ 3],   
  -1 + (c[1] + 0.808085 c[3] + 0.589066 s[3]) y2[ 1], (-0.70548 c[1] +   
  0.70873 s[1] + s[2]) y1[ 1] - (-0.589066 c[3] + s[1] + 0.808085 s[3])  
  y2[1], -1 + (-0.70873 c[1] + c[2] - 0.70548 s[1]) y2[ 2], (0.399814 c[2] + 
  0.916596 s[2] + s[3]) y1[ 2] - (-0.70548 c[1] + 0.70873 s[1] + s[2]) y2[ 2],   
 -1 + (-0.916596 c[2] + c[3] + 0.399814 s[2]) y2[ 3], (0.589066 c[3] - s[1] - 
 0.808085 s[3]) y1[ 3] - (0.399814 c[2] + 0.916596 s[2] + s[3]) y2[3]}

using NSolve command. i.e.,

sol = NSolve[e==0, Variables[e]]

I am only interested in real solutions, so in principle I could also use

NSolve[e==0, Variables[e], Reals].

The command indeed solves this system and gives me 8 solutions, all real. However, when I insert them back to the equations to check if the solutions satisfy the equations, they don't!

Do[ Print[eatsol[i] = e /. Thread[var -> sol[[i]]]], {i, lrealsol}]



Thus, the first 6 solutions don't even remotely satisfy the equations! I would be fine with some numerical errors etc. but it seems that the infinite solutions also appear here.

Do you know what is going on here? Is there any problem with NSolve?

Update: Please run the below code in Mathematica:

e = {-1 + c[1]^2 + s[1]^2, -1 + c[2]^2 + s[2]^2, -1 + c[3]^2 + 
    s[3]^2, -1 + (-0.5646917067485537` c[1] + c[2] - 
       0.8253019304045068` s[1]) y1[1], -1 + (0.5033628291371058` c[2] + c[3] - 
       0.8640751484929357` s[2]) y1[2], -1 + (c[1] - 0.47880555467239877` c[3] + 
       0.877920976406679` s[3]) y1[3], (-0.8253019304045068` c[1] + 
       0.5646917067485537` s[1] + s[2]) y1[1] - (-0.877920976406679` c[3] + s[1] - 
       0.47880555467239877` s[3]) y2[1], -1 + (c[1] - 0.47880555467239877` c[3] + 
       0.877920976406679` s[3]) y2[1], -1 + (-0.5646917067485537` c[1] + c[2] - 
       0.8253019304045068` s[1]) y2[2], (-0.8640751484929357` c[2] - 
       0.5033628291371058` s[2] + s[3]) y1[2] - (-0.8253019304045068` c[1] + 
       0.5646917067485537` s[1] + s[2]) y2[2], -1 + (0.5033628291371058` c[2] + c[3] - 
       0.8640751484929357` s[2]) y2[3], (0.877920976406679` c[3] - s[1] + 
       0.47880555467239877` s[3]) y1[3] - (-0.8640751484929357` c[2] - 
       0.5033628291371058` s[2] + s[3]) y2[3]};
var = Variables[e];
sol = NSolve[e == 0, var];
lsol = Length[sol];
 Print[Chop[e /. sol[[i]]]],
 {i, lsol}]

Basically, the problem is that the solutions given by the NSolve command are not actually the solutions of the equations!

  • $\begingroup$ Try e /. sol // Chop. $\endgroup$ Oct 31, 2012 at 20:58
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2)Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign` $\endgroup$
    – chris
    Oct 31, 2012 at 20:59
  • $\begingroup$ or if you want to be super sure, 'ee=e// Rationalize[#, 0.00000001] &; sol = Solve[Thread[ee == 0], Variables[ee]]; ee /. sol // N[#, 20] & // Chop' $\endgroup$
    – chris
    Oct 31, 2012 at 21:05
  • $\begingroup$ b.gatessucks and chris, thanks for your comments. My problem is different actually. The problem is that some of the equations at the first 6 solutions is -1, and not zero (or some small number). So Chop etc. is not going to help.d $\endgroup$
    – dbm
    Oct 31, 2012 at 21:12
  • 5
    $\begingroup$ Folks, please stop downvoting this. it has actual merit and at this point is more or less clear. $\endgroup$ Oct 31, 2012 at 23:36

3 Answers 3


The issue is numeric instability. If you rationalize the equations and solve at higher precision you will get small residuals.

re = Rationalize[e, 0];
var = Variables[e];
sol50 = NSolve[re == 0, var, WorkingPrecision -> 50];
lsol = Length[sol];

Max[Abs[re /. sol50]]

(* Out[208]= 4.8358090074041686763439090398`1.0824091485306464*^-49 *)

If you check sizes of variables you will notice many orders of magnitude separate them in some solutions. This is an indication that high precision may be needed in order to get significant digits for the smaller ones.

--- edit ---

Okay, I missed this subtlety of "infinite" solutions. They are indeed infinite, in an important sense (and this gives the numeric instability). In a certain way the system is not "generic". Specifically, coefficients satisfy one or more (unknown to me) relevant algebraic conditions that put the system on what is called the "discriminant variety" for the family of systems having the same Newton polytope. Remark: Since the system is generic for your purposes, that means your entire family of systems lives on this variety.

In such a situation the number of solutions can drop from the generic (in this case) 8. But that does not quite happen for the example above, because we have coefficients given only approximately, via machine doubles. So we actually have a system NEAR but not ON the discriminant variety.

So what happens to the 8 solutions when we are on this variety? In this case 6 have, loosely speaking, wandered off to infinity. Continuity of solutions as a function of coefficients implies that systems near the variety will have 6 "large" solutions. That's what NSolve is showing. When I round and solve to high precision I get a perfectly good solution set, albeit residuals of the "bad" solutions will be orders of magnitude larger than those of the good ones.

Homotopy continuation methods will not show behavior this for either of two reasons.

(1) They detect "diverging" solutions via heuristics in the tracking phase. In this case 6 would be dropped once they got too large in the tracking.

(2) In some cases one already has a correct result for a related system from the same part of the discriminant variety. In this case a common tactic is to use what is called the "Cheater's Homotopy", using these old solutions to go to solutions for the new system. Here we know in advance that we have the right number of starting points and the tracking will work just fine (barring other numeric issues).

NSolve has a system option that can be used, in a heuristic manner, to elicit better behavior. Or worse, depending on mood. The idea is implementation dependent: we compute a numerical Groebner basis, and a nondefault setting will tell that code to set to zero any coefficients that are below some threshold smaller than the average size in the rest of the polynomial. I show how to use this below, with 10^(-6) for that tolerance. Remark: if I go to 10^(-4) it will hang. That's one way the heuristic can mess up. Another is that it can give garbage results. It is very much a situation where one needs to make sure the results (when they appear) seem plausible. Even finding a viable tolerance is a trial-and-error process, informed perhaps by the nature of the problem. Here I know we have approximated coefficients (only) to several digits, so I expect that 10^(-8) or so will work for my purpose. And it did, as did 10^(-6).

SetSystemOptions["NSolveOptions" -> {"Tolerance" -> 10^(-6)}];
sol = NSolve[e == 0, var]

(* {{c[1] -> -0.840799, c[2] -> 0.0280179, c[3] -> 0.877839, 
  s[1] -> 0.541347, s[2] -> -0.999607, s[3] -> -0.478956, 
  y1[1] -> 17.8457, y1[2] -> 0.56958, y1[3] -> -0.594672, 
  y2[1] -> -0.594672, y2[2] -> 17.8457, 
  y2[3] -> 0.56958}, {c[1] -> 0.840799, c[2] -> -0.0280179, 
  c[3] -> -0.877839, s[1] -> -0.541347, s[2] -> 0.999607, 
  s[3] -> 0.478956, y1[1] -> -17.8457, y1[2] -> -0.56958, 
  y1[3] -> 0.594672, y2[1] -> 0.594672, y2[2] -> -17.8457, 
  y2[3] -> -0.56958}} *)

Here is the example where I first ran into this near-discriminant-variety computational mishap four years ago.


It took some head-scratching between myself and the web page author to figure out what had gone wrong.

--- end edit ---

  • $\begingroup$ That worked! However, if I still think there is more fundamental problem here: NSolve gives solutions at infinity. If you see the actual solutions in sol, some variables take values of the order 10^16. These are solutions at infinity for a reasonable numerical precision say 10^13 or so. So there is a catch22 here. If we choose WorkingPrecision->50, the solutions at infinity are considered as affine solutions, if we choose the default WorkingPrecision then these solutions still appear but don't satisfy the equations! May be I am explaining it not properly?! $\endgroup$
    – dbm
    Oct 31, 2012 at 23:04
  • $\begingroup$ I'll edit my response tomorrow to address that issue of "infinite" solutions. Meanwhile you might mention where the system came from or the nature of the problem. That might help when I add an explanation. $\endgroup$ Oct 31, 2012 at 23:39
  • $\begingroup$ Thanks for your comment and your request to others to stop down-voting. The problem comes from particle theory where I have a nonlinear parametric system of equations and I want to find all real solutions of it. In particular, I want to see if the no. of real solutions is constant for the parametric system. The present system is at a particular parameter-point (a generic point). I started having doubt about the infinite solutions was from my cross-checking the solutions with homotopy continuation method which gave only two solutions always instead of 8 by NSolve. Thanks again. $\endgroup$
    – dbm
    Oct 31, 2012 at 23:57
  • $\begingroup$ Edited to add explanation. $\endgroup$ Nov 1, 2012 at 14:19
  • $\begingroup$ @Licthblau, Excellent! Thanks a lot for your detailed explanation. In short, we need to be extra careful while using NSolve for systems with approximate coefficients. I already knew discriminant variety and homotopy continuation method. But I didn't know about the near discriminant variety issue, and specially that it could affect solving systems so drastically. Thanks a lot again! dbm368 $\endgroup$
    – dbm
    Nov 1, 2012 at 16:40
ss = NSolve[Thread[e == 0], Variables@e, WorkingPrecision -> 50];
Select[Abs /@ Join @@ (e /. ss ), # > 10^-9 &]
  • $\begingroup$ Please see the update in my original massage. $\endgroup$
    – dbm
    Oct 31, 2012 at 22:38
  • $\begingroup$ @dbm Please see my answer $\endgroup$ Oct 31, 2012 at 22:55

I wanted to comment a "Thank You" to @dbm for asking and @Daniel Lichtblau for answering, but don't have the necessary reputation yet.

There's also a way to bypass this problem in certain cases that wasn't explicitly stated here, so I figured I'd mention it.

I ran into this same issue when randomly generating momenta (also for a particle physics problem), and so had $k_1 + k_2 + k_3 + k_4 = 0$, where $k_i$ are four-component vectors also satisfying $k_i^2 = 0$. These were the explicit algebraic constraints that were putting me on or off (not sure which) the discriminant variety, per @Daniel Lichtblau's answer:

Specifically, coefficients satisfy one or more (unknown to me) relevant algebraic conditions that put the system on what is called the "discriminant variety" for the family of systems having the same Newton polytope.

The tiny numerical error I got from solving the constraint equations $$ k_1 + k_2 + k_3 + k_4 = 0 \quad \textrm{and} \quad k_i^2 = 0 $$ numerically was enough to result in the infinite solutions described by @dbm. The solution in my specific case (and likely in many cases where the algebraic constraints can be handled explicitly), was to impose the constraints analytically instead of numerically. That reduced the numerical error enough that NSolve returned only the non-infinite solutions, all of which satisfied the original equations.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.