Re: Is this a bug in NSolve in mathematica?

*To*: mathgroup at smc.vnet.net*Subject*: [mg113237] Re: Is this a bug in NSolve in mathematica?*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Wed, 20 Oct 2010 04:07:10 -0400 (EDT)

Scott Morrison wrote: > One would expect and hope that if you ask Mathematica to find the > roots of a polynomial, it should give the same (approximate) answers > whether you do this symbolically, then find numerical approximations > to these exact answers, or whether you do it numerically. Here's an > example which (in Mathematica 7, running on OS X) where this fails > badly: > > --------------------------------------------------------------------------------------------------------------------------- > poly = -112 + 1/q^28 + 1/q^26 - 1/q^24 - 6/q^22 - 14/q^20 - 25/q^18 - > 38/q^16 - 52/q^14 - 67/q^12 - 81/q^10 - 93/q^8 - 102/q^6 - 108/ > q^4 - 111/q^2 - 111 q^2 - 108 q^4 - 102 q^6 - 93 q^8 - 81 q^10 - > 67 q^12 - 52 q^14 - 38 q^16 - 25 q^18 - 14 q^20 - 6 q^22 - q^24 + > q^26 + q^28; > > Total[q^4 /. NSolve[poly == 0, q]] - Total[q^4 /. N[Solve[poly == 0, > q]]] > --------------------------------------------------------------------------------------------------------------------------- > > (Note: this is actually a Laurent polynomial, and if you multiply > through by a large power of q the problem goes away.) > > The last line here is just a demonstration that the solutions found > are very different; in fact it's the quantity we were trying to > compute in the problem we were working on. > > If you look closely at the output of NSolve[poly == 0, q] and of > N[Solve[poly == 0, q], you'll see that NSolve only gives 54 roots > instead of the expected 56. It's not that it just missed a repeated or > anything; it's missing the two largest roots in magnitude > (approximately +/- 1.59) > > Is this a bug in Mathematica? Does anyone have an explanation for why > this is happening? > > (ps -- This is crossposted on http://stackoverflow.com/questions/3965667/is-this-a-bug-in-nsolve-in-mathematica, > if anyone would prefer to answer there.) This is in the nether zone between bug and (numeric) feature. It arises because NSolve has a need to verify the result (due to presence of denominators), coupled with modest conditioning problems that lead to bad residuals for the two roots in question. Using machine precision, those roots are {{q -> 1.5922762814263167}, {q -> -1.5922421375870808}} If you plug these into poly you will see that the residuals are large. Indeed, the results shown are only correct to four or so places, and this accounts for much of the problem. You can work around this using NSolve as below. One simply specifies a working precision higher than MachinePrecision. rts1 = Sort[q/.NSolve[poly == 0, q, WorkingPrecision->25]] In this case the large roots are {{q -> -1.59223869382733814141818983865564857786`24.69897000433602}, {q -> 1.59223869382734244429358140525375579006`24.69897000433602}} and they will be retained in the verification process. I will look into ways to improve on this for a future release. Daniel Lichtblau Wolfram Research