Polynomial roots#

Polynomial roots are nice because we can solve for them directly. Polynomials of degree \(n\) have exactly \(n\) roots, including multiplicity (coincident roots). The roots may be complex, even in polynomials with real coefficients. If the root is complex (\(a+bi\)), then its complex conjugate (\(a-bi\)) is also a root.

The roots of polynomials of order \(\le 4\) can be solved analytically.

E.g.:

\(a x^2 +b x + c = 0\)

has roots:

\(x = \frac{-b \pm \sqrt{b^2-4ac}} {2a}\)

per the quadratic formula.

Certaint polynomials of order \(\gt 4\) may have analytically solvable roots but there is no such general formula per the Abel-Ruffini theorem.

All the roots of polynomials may be found numerically by solving for the eigenvalues of the polynomial companion matrix. This is the basis for modern numerical methods which we will cover when we get to matrix eigenvalue calculations.

Example usage of numerical root finding tools.#

import numpy as np
from scipy.linalg import companion

# Define the coefficients of the polynomial
coefficients = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]

# Print the polynomial
print("Polynomial:")
print(np.poly1d(coefficients))

# Calculate the companion matrix
companion_matrix = companion(coefficients)

# Print the companion matrix
print("\nCompanion Matrix:")
print(np.round(companion_matrix, 2))

# Find the eigenvalues (roots) of the companion matrix
roots = np.linalg.eigvals(companion_matrix)

# Print the roots
print("\nRoots from eigenvalues of companion matrix:")
print(np.round(roots, 2))

print("\nRoots from 'roots' function:")
print(np.round(np.roots(coefficients),2))
Polynomial:
   10     9     8     7     6     5     4     3     2
1 x  + 2 x + 3 x + 4 x + 5 x + 6 x + 7 x + 8 x + 9 x + 10 x + 11

Companion Matrix:
[[ -2.  -3.  -4.  -5.  -6.  -7.  -8.  -9. -10. -11.]
 [  1.   0.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   1.   0.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   1.   0.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   1.   0.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   1.   0.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   1.   0.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   1.   0.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   1.   0.   0.]
 [  0.   0.   0.   0.   0.   0.   0.   0.   1.   0.]]

Roots from eigenvalues of companion matrix:
[-1.26+0.36j -1.26-0.36j -0.88+0.96j -0.88-0.96j -0.25+1.26j -0.25-1.26j
  0.44+1.17j  0.44-1.17j  0.95+0.73j  0.95-0.73j]

Roots from 'roots' function:
[-1.26+0.36j -1.26-0.36j -0.88+0.96j -0.88-0.96j -0.25+1.26j -0.25-1.26j
  0.44+1.17j  0.44-1.17j  0.95+0.73j  0.95-0.73j]

Note: in many software environments, \(j\) is used as the imaginary number instead of \(i\).

Finding roots of the companion matrix via the eigenvalue is a powerful and standard approach. In particular it:

  • finds all the roots simultaneously

  • does not require initial guesses

  • is numerically robust thanks to robust eigenvalue algorithms

However, it is limited to appropriate scalar functions, including:

  • standard polynomials (including taylor expansions!)

  • orthogonal polynomials (Legendre, Hermite, Laguerre,… )

  • Fourier series expansions (via complex representation and a change of variables to make it a polynomial)