An algorithm that I'm trying to implement requires finding the roots of a 10th degree polynomial, which I created with sympy, it looks like this:
import sympy
import numpy as np
det = sympy.Poly(1.3339507303385e-16*z**10 + 6.75390067076469e-14*z**9 + 7.18791227134351e-12*z**8 + 2.27504286959352e-10*z**7 + 2.37058998324426e-8*z**6 + 1.63916629439745e-6*z**5 + 3.0608041245671e-5*z**4 + 4.83564348562906e-8*z**3 + 2.0248073519853e-5*z**2 - 4.73126177166988e-7*z + 1.1495883206077e-6)
For finding the roots of the polynomial, I use the following code:
coefflist = det.coeffs()
solutions = np.roots(coefflist)
print(coefflist)
[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
print(solutions)
[-3.70378229e+02+0.00000000e+00j -1.18366138e+02+0.00000000e+00j
2.71097137e+01+5.77011644e+01j 2.71097137e+01-5.77011644e+01j
-3.59084863e+01+1.44819591e-02j -3.59084863e+01-1.44819591e-02j
2.60969082e-03+7.73805425e-01j 2.60969082e-03-7.73805425e-01j
1.42936329e-02+2.49877948e-01j 1.42936329e-02-2.49877948e-01j]
However, when I substitute z with a root, lets say the first one, the result is not zero, but some number:
print(det.subs(z,solutions[0]))
-1.80384169514123e-6
I would have expected, that the result probably isn't the integer 0, but 1e-6 is pretty bad (it should be zero, right?). Is there a mistake in my code? Is this inaccuracy normal? Any thoughts/suggestions would be helpful. Is there a more accurate alternative to compute the roots of a 10th degree polynomial?
You do not need sympy, the methods in numpy are completely sufficient. Define the polynomial by its list of coefficients and compute the roots
p=[1.33395073033850e-16, 6.75390067076469e-14, 7.18791227134351e-12, 2.27504286959352e-10, 2.37058998324426e-8, 1.63916629439745e-6, 3.06080412456710e-5, 4.83564348562906e-8, 2.02480735198530e-5, -4.73126177166988e-7, 1.14958832060770e-6]
sol= np.roots(p); sol
giving the result
array([ -3.70378229e+02 +0.00000000e+00j, -1.18366138e+02 +0.00000000e+00j,
2.71097137e+01 +5.77011644e+01j, 2.71097137e+01 -5.77011644e+01j,
-3.59084863e+01 +1.44819592e-02j, -3.59084863e+01 -1.44819592e-02j,
2.60969082e-03 +7.73805425e-01j, 2.60969082e-03 -7.73805425e-01j,
1.42936329e-02 +2.49877948e-01j, 1.42936329e-02 -2.49877948e-01j])
and evaluate the polynomials at these approximate roots
np.polyval(p,sol)
giving the array
array([ 2.28604877e-06 +0.00000000e+00j, 1.30435230e-10 +0.00000000e+00j,
1.05461854e-11 -7.56043461e-12j, 1.05461854e-11 +7.56043461e-12j,
-3.98439686e-14 +6.84489332e-17j, -3.98439686e-14 -6.84489332e-17j,
1.18584613e-20 +1.59976730e-21j, 1.18584613e-20 -1.59976730e-21j,
6.35274710e-22 +1.74700545e-21j, 6.35274710e-22 -1.74700545e-21j])
Obviously, evaluating a polynomial close to a root involves lots of catastrophic cancellations, that is, the intermediate terms are large of opposite sign and cancel out, but their errors are proportional to their original sizes. To get an estimate of the combined error size, replace the polynomial coefficients with their absolute values and also the evaluation points.
np.polyval(np.abs(p),np.abs(sol))
resulting in
array([ 1.81750423e+10, 8.40363409e+05,
8.08166359e+03, 8.08166359e+03,
2.44160616e+02, 2.44160616e+02,
2.50963696e-05, 2.50963696e-05,
2.65889696e-06, 2.65889696e-06])
In the case of the first root, the scale multiplied with the machine constant gives an error scale of 1e+10*1e-16=1e-6, which means that the value at the root is as good as zero within the framework of double-precision floating point.