Is this a bug in NSolve in mathematica?

Question

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 root 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?

Answer 1

Unfortunately your expectations are unfounded.

As you say, Solve[] gives an exact solution and N[] introduces a small error, but only once, when you evaluate it. NSolve[] on the other hand uses numerical approximations from the get go and thus significant rounding error can accumulate.

You are also limited by the default precision for the calculation and this can lead to complete failure of the numerical method, e.g., missing roots (see Wilkinson's polynomial). You can counteract this by telling NSolve[] to use greater precision like so:

In[1]  := Total[q^4 /. NSolve[poly == 0, q, WorkingPrecision -> 50]] - 
          Total[q^4 /. N[Solve[poly == 0, q]]]
Out[1] := 0. - 3.66374*10^-15 I

When using numerical methods it is always important to keep the errors in mind. Since this is true for for a large variety of numerical analysis problems from solving long polynomials to diagonalizing large matrices to integration of weird functions etc. etc. there is no one correct approach, and Mathematica needs to be told to, e.g., raise the WorkingPrecision, or to apply a different numerical technique.