Roots given by Solve are not satisfied by the equation

Question

I am finding roots of a polynomial equation using Solve. However, two of the five roots given by Solve do not satisfy the equation. Here how I do it:

delta = 10.0;
lund = 10^4;
alpha = 10^(-7);
g = 10.0;
eta = 1.0;
VA = lund*eta/delta;
k0 = Sqrt[6]/(1.0*delta);

ks = 1/Sqrt[2] k0;
kz = Sqrt[k0^2 - ks^2];
k = Sqrt[ks^2 + kz^2];

beta = 10^(-2)*VA^2*kz^2*k^2/(1.0*alpha*g*ks^2);

Omega = VA/(2.0*0.0003*delta);
omgM = VA*kz;
omge = eta*k^2;
omgA2 = (Sqrt[-g*alpha*beta]*ks/k)^2;
omgO = 2*Omega*kz/k;

eqn[l_] = 
  I l^5 + 2 l^4 omge + omgA2 omge omgM^2 - 
   2 l^2 omge (omgA2 + omgM^2 + omgO^2) - 
   I l^3 (omgA2 + omge^2 + 2 omgM^2 + omgO^2) + 
   I l (omgM^4 + omgA2 (omge^2 + omgM^2) + omge^2 omgO^2);

Solve[eqn[l] == 0, l]

The roots are:

{{l -> -235702. + 3.22907*10^-8 I}, {l -> -0.126759 + 
    0.0602467 I}, {l -> 0. - 0.000493465 I}, {l -> 
   0.126759 + 0.0602467 I}, {l -> 235702. + 3.23999*10^-8 I}}

The first and last roots do not satisfy the equation:

eqn[-235702 + 3.23*10^(-8) I]

gives

-1.21572*10^15 + 2.38912*10^21 I

Are the roots not correct? or Am I missing something?

Answer 1

You are running into a problem because of a loss of numerical precision in your calculations.

Let me change your code to use arbitrary precision throughout:

delta = 10;
lund = 10^4;
alpha = 10^(-7);
g = 10;
eta = 1;
VA = lund*eta/delta;
k0 = Sqrt[6]/delta;

ks = 1/Sqrt[2] k0;
kz = Sqrt[k0^2 - ks^2];
k = Sqrt[ks^2 + kz^2];

beta = 10^(-2)*VA^2*kz^2*k^2/(alpha*g*ks^2);

Omega = VA/(6/10000*delta);
omgM = VA*kz;
omge = eta*k^2;
omgA2 = (Sqrt[-g*alpha*beta]*ks/k)^2;
omgO = 2*Omega*kz/k;

eqn[l_] = 
  I l^5 + 2 l^4 omge + omgA2 omge omgM^2 - 
   2 l^2 omge (omgA2 + omgM^2 + omgO^2) - 
   I l^3 (omgA2 + omge^2 + 2 omgM^2 + omgO^2) + 
   I l (omgM^4 + omgA2 (omge^2 + omgM^2) + omge^2 omgO^2);

Let's now Solve at arbitrary precision, then obtain the numerical values of the roots with N, using arbitrary-precision math with a number of digits of precision equivalent to machine-precision numbers, but all the while keeping track of precision using the arbitrary precision machinery:

N[
  eqn[l] /. Solve[eqn[l] == 0, l], 
  $MachinePrecision
]

(* Out: {0.*10^-43 + 0.*10^-42 I, 
         0.*10^-61 + 0.*10^-61 I, 
         0.*10^-64 + 0.*10^-65 I, 
         0.*10^-61 + 0.*10^-61 I, 
         0.*10^-43 + 0.*10^-42 I}  *)

Chop[%]

(* Out: {0, 0, 0, 0, 0} *)

Answer 2

The roots returned by Solve are nearly as accurate as possible within rounding error at machine precision. Even if you take @Marco's exact solution and convert the roots to machine precision (N[sol]), they won't be significantly more accurate. The reason is that given a polynomial $y=\sum a_k x^k$ expanded in powers, a rounding error in $y$ would be up to about the max of $|a_k x^k| \, \epsilon$, where $\epsilon \approx 2\cdot10^{-16}$ is given by $MachineEpsilon in Mathematica. For the largest monomial term, we have $|a_k x^k| \approx 7\cdot10^{26}$ for the root $x \approx 2\cdot10^5$, and the error in $y$ could be up to $1.6\cdot10^{11}$.

Here are the rounding-error bounds of the largest monomial terms and the residual errors in "y" for each root, which are all about the same size:

Abs[List @@ eqn[l]*$MachineEpsilon] /. sol // Map@Max
(*  {1.6153*10^11, 3.4102*10^-8, 1.1990*10^-10, 3.4102*10^-8, 1.6153*10^11}  *)

eqn[l] /. sol // Abs
(*  {6.8869*10^11, 3.0719*10^-8, 2.0145*10^-11, 7.4505*10^-9, 1.3451*10^11}  *)

If you're interested in determining the roots to an ordinary number of digits, which might be 3, or 6, or 10, depending on the field, then the results of Solve are satisfactory. You won't get a significantly different result for the roots using higher precision. If you're interested in having enough precision to make the residual be "small," whatever that means, then you will need a lot of extra precision. To get the residual to be less than 1, you will need at least 12 more digits, or a precision of 16+12 = 28:

SetPrecision[eqn[l], 16 + 12];
% /. Solve[% == 0, l] // Abs
(*  {0.*10^1, 0.*10^-18, 0.*10^-22, 0.*10^-18, 0.*10^1}  *)

The 0.*10^1 indicates that there was some extra rounding error and 12 extra digits was not quite enough: The rounding error bound is somewhere around 10. So 29-30 digits should be enough, but I wonder how often a 30-digit answer is needed.