Polynomial roots - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
from scipy.linalg import companion

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

from executable_engineering import plot_polynomial_roots_companion

Explicit formulae¶

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.

Certain polynomials of order $> 4$ may have analytically solvable roots but there is no such general formula per the Abel-Ruffini theorem.

Polynomial roots via the companion matrix¶

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.

The companion matrix is formed by placing the negated, normalized coefficients of the polynomial in the top row, and placing 1s along the first subdiagonal. (If the polynomial’s leading coefficient is not 1, all coefficients are first divided by the leading coefficient).

For example, consider a simple monic cubic polynomial (leading coefficient is 1):

P(x) = x^3 + c_2 x^2 + c_1 x + c_0

(1)

The companion matrix $\mathbf{C}$ is constructed as:

\mathbf{C} = \begin{bmatrix} -c_2 & -c_1 & -c_0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}

(2)

Recovery of the original polynomial¶

To see why this works, consider the eigenvalues $\lambda$ , which we can solve from the characteristic equation $\det(\lambda \mathbf{I} - \mathbf{C}) = 0$ :

\det(\lambda \mathbf{I} - \mathbf{C}) = \det \begin{bmatrix} \lambda + c_2 & c_1 & c_0 \\ -1 & \lambda & 0 \\ 0 & -1 & \lambda \end{bmatrix}

(3)

Expanding the determinant along the third column:

\begin{aligned} \det(\lambda \mathbf{I} - \mathbf{C}) &= c_0 \det \begin{bmatrix} -1 & \lambda \\ 0 & -1 \end{bmatrix} - 0 + \lambda \det \begin{bmatrix} \lambda + c_2 & c_1 \\ -1 & \lambda \end{bmatrix} \\ &= c_0 \big( (-1)(-1) - (\lambda)(0) \big) + \lambda \big( (\lambda+c_2)(\lambda) - (c_1)(-1) \big) \\ &= c_0 + \lambda(\lambda^2 + c_2 \lambda + c_1) \\ &= \lambda^3 + c_2 \lambda^2 + c_1 \lambda + c_0 \end{aligned}

(4)

Because the characteristic polynomial of this matrix precisely recovers the original polynomial $P(\lambda)$ , finding the eigenvalues of the companion matrix is mathematically equivalent to finding the roots of the polynomial!

Example: Polynomial roots via the companion matrix¶

The following tools generates a random matrix of degree $n$ , shows the associated companion matrix and calculates its eigenvalues.

N.B. in many software environments, $j$ is used as the imaginary number instead of $i$ .

plot_polynomial_roots_companion(4)

Polynomial:
       4         3         2
8.978 x + 9.313 x + 6.168 x - 3.908 x - 8.047

Companion Matrix:
[[-1.04 -0.69  0.44  0.9 ]
 [ 1.    0.    0.    0.  ]
 [ 0.    1.    0.    0.  ]
 [ 0.    0.    1.    0.  ]]

Roots from eigenvalues of companion matrix:
[-0.44+1.03j -0.44-1.03j -0.92+0.j    0.77+0.j  ]

Loading...

Roots of a Fourier Series¶

The companion matrix method can also be cleverly adapted to find the roots of a Fourier series. The core idea is to use Euler’s formula and a change of variables to transform a finite Fourier series into a standard algebraic polynomial.

1. Express as a Complex Exponential Series Consider a finite Fourier series of degree $N$ :

f(x) = a_0 + \sum_{k=1}^N \left( a_k \cos(kx) + b_k \sin(kx) \right)

(5)

Using Euler’s formula ( $\cos(kx) = \frac{e^{ikx} + e^{-ikx}}{2}$ and $\sin(kx) = \frac{e^{ikx} - e^{-ikx}}{2i}$ ), this can be rewritten as a sum of complex exponentials:

f(x) = \sum_{k=-N}^N c_k e^{ikx}

(6)

(where the complex coefficients $c_k$ are derived from $a_k$ and $b_k$ ).

2. Change of Variables We introduce a new complex variable, $z = e^{ix}$ . This means $e^{ikx} = z^k$ , and $e^{-ikx} = z^{-k}$ . Substituting this into our series yields a Laurent polynomial:

F(z) = \sum_{k=-N}^N c_k z^k = 0

(7)

3. Convert to a Standard Polynomial To eliminate the negative exponents and turn this into a standard polynomial, we simply multiply the entire equation by $z^N$ (since $z = e^{ix}$ , we know $|z| = 1 \neq 0$ ):

z^N F(z) = \sum_{k=-N}^N c_k z^{k+N} = \sum_{m=0}^{2N} c_{m-N} z^m = 0

(8)

We now have a standard polynomial of degree $2N$ .

4. Solve and Map Back to the Original Domain Because it is a standard polynomial, we can construct its $2N \times 2N$ companion matrix and compute its eigenvalues to find all $2N$ complex roots for $z$ .

Finally, we must map our roots $z_j$ back to the original variable $x$ . Since $z = e^{ix}$ , any strictly real root $x$ will result in a root $z$ that lies exactly on the unit circle in the complex plane ( $|z| = 1$ ).

We filter the eigenvalues for those where $|z| \approx 1$ (allowing for minor numerical floating-point errors). The real roots of the original Fourier series are simply the angles of these complex numbers:

x = \arg(z)

(9)

Summary of the companion matrix approach¶

This method has the benefits that 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 that are expressed in a linear basis:

standard polynomials (including taylor expansions!)
orthogonal polynomials (Legendre, Hermite, Laguerre,... )
Fourier series expansions