Polynomial roots¶

For this example we will use a publicly available package for working with polynomials. It should be available using the following line, if you have followed installation instructions for these scripts.

using Polynomials

Our first step is to construct a polynomial with six known roots.

r = [-2.0,-1,1,1,3,6]
@show p = fromroots(r);

p = fromroots(r) = Polynomial(36.0 - 36.0*x - 43.0*x^2 + 44.0*x^3 + 6.0*x^4 - 8.0*x^5 + 1.0*x^6)

Now we use a standard numerical method for finding those roots, pretending that we don’t know them already.

@show r_computed = sort(roots(p));

r_computed = sort(roots(p)) = [-2.0000000000000013, -0.999999999999999, 0.9999999902778504, 1.0000000097221495, 2.9999999999999996, 5.999999999999992]

Here are the relative errors in each of the computed roots. The @. notation at the start means essentially to do the given operations on each element of the given vectors.

@. abs(r - r_computed) / r

6-element Array{Float64,1}:
 -6.661338147750939e-16
 -9.992007221626409e-16
  9.722149640900568e-9
  9.722149529878266e-9
  1.4802973661668753e-16
  1.3322676295501878e-15

It seems that the forward error is acceptably close to machine epsilon for double precision in all cases except the double root at \(x=1\). This is not a surprise, though, given the poor conditioning at such roots.

Let’s consider the backward error. The data in the rootfinding problem are the polynomial coefficients. We can apply fromroots to find the coefficients of the polynomial (that is, the data) whose roots were actually computed by the numerical algorithm.

@show p_computed = fromroots(r_computed);

p_computed = fromroots(r_computed) = Polynomial(35.99999999999993 - 35.999999999999915*x - 42.99999999999997*x^2 + 43.99999999999998*x^3 + 5.999999999999973*x^4 - 7.999999999999991*x^5 + 1.0*x^6)

We find that in a relative sense, these coefficients are very close to those of the original, exact polynomial:

cp = coeffs(p)
cpc = coeffs(p_computed)
@. abs(cp-cpc)/cp

7-element Array{Float64,1}:
  1.973729821555834e-15
 -2.3684757858670005e-15
 -6.609699867535816e-16
  4.844609562000683e-16
  4.440892098500626e-15
 -1.1102230246251565e-15
  0.0

In summary, even though there are some computed roots relatively far from their correct values, they are nevertheless the roots of a polynomial that is very close to the original.

Fundamentals of Numerical Computation

Polynomial roots¶