# Roots given by Solve are not satisfied by the equation

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?

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} *) • @MacroB thank you so much. Jan 9, 2021 at 17:18 • @MichealE2 Thank you for explaining it so nicely. Jan 11, 2021 at 17:49 • @GauravMaurya You're welcome. I guess you meant to comment on my post, though. :) Jan 11, 2021 at 18:03 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. • To address your last remark, often it is not needed. Although it can be useful for assessing whether a lower-precision result is or is not plausible (as you and @MarcoB have in effect shown). Jan 10, 2021 at 0:04 • @DanielLichtblau Yes, I do that all the time and on many answers on the site. My main objective is to dispel the impression that Solve or NSolve is faulty, by giving a reason why 10^11 can be considered a small error in$y$in this case. Another source of error I omitted is that the error in$y$is$dy \approx f'(x)\,dx$(e.g. eqn'[l]*l*$MachineEpsilon /. sol // Abs), which gives similar bounds. Maybe a third thing, which I should have emphasized, is that round-off error in the coefficients leads to a large $y$ residual. Questions on root-finding on this site often ignore these factors. Jan 11, 2021 at 17:58
• When I first looked at the question, I was going to remark about using an exact solver with approximate input, but @Marco beat me to it. Then I realized it wasn't an issue here, and that "a loss of numerical precision" did not explain what seemed to be wrong. Jan 11, 2021 at 18:00
• One way I assess residuals is to go to scale-free coordinates. (1) Rescale equations so that max coefficient is 1 (or norm of coeffs is 1, take your pick). (2) Homogenize the equations, (3) Homogenize the results. (4) Rescale results so max coordinate is 1 (or norm...). (5) Now assess residuals. This approach helps to avoid issues with scaling of the equations and also with relative size of solutions. So "large" solutions will be less likely to give outsized residuals. Jan 12, 2021 at 0:50