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.)