Re: Is this a bug in NSolve in mathematica?

*To*: mathgroup at smc.vnet.net*Subject*: [mg113240] Re: Is this a bug in NSolve in mathematica?*From*: Bob Hanlon <hanlonr at cox.net>*Date*: Wed, 20 Oct 2010 04:07:44 -0400 (EDT)

Just use higher precision 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; prec = 25; Total[q^4 /. NSolve[poly == 0, q, WorkingPrecision -> prec]] - Total[q^4 /. N[Solve[poly == 0, q], prec]] // Chop 0 $Version "7.0 for Mac OS X x86 (64-bit) (February 19, 2009)" Length[Solve[poly == 0, q]] 56 Length[N[Solve[poly == 0, q]]] 56 Length[NSolve[poly == 0, q, WorkingPrecision -> 25]] 56 Bob Hanlon ---- Scott Morrison <scott.morrison at gmail.com> 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.)