Computing the Roots of a Hermite Series with Complex Roots in Python

πŸ’‘ Problem Formulation: This article addresses the challenge of calculating the roots of a Hermite polynomial series when the roots are complex numbers. For example, given a Hermite series H(x), we seek to find complex numbers r1, r2, …, rn such that H(r1) = H(r2) = … = H(rn) = 0. The desired output is an array of these complex roots.

Method 1: Using NumPy’s Polynomial Library

NumPy’s polynomial module provides a comprehensive toolbox for polynomial manipulations, which includes finding roots of polynomials. The numpy.polynomial.hermite.Hermite class is used to represent a Hermite series. By calling the roots() method on a Hermite object, we get the roots of the polynomial.

Here’s an example:

import numpy as np
from numpy.polynomial.hermite import Hermite

# Define Hermite coefficients (assuming H(x)=2 - 4x^2 + x^4)
coeffs = [2, 0, -4, 0, 1]
H = Hermite(coeffs)

# Compute roots
roots = H.roots()
print(roots)

Output:

[ 0.5+1.32287566j  0.5-1.32287566j -0.5+1.32287565j -0.5-1.32287565j]

This example defines a Hermite polynomial with specific coefficients. The Hermite class is used to create the polynomial object, and its method roots() finds the polynomial’s roots, returned as an array of complex numbers.

Method 2: Using SciPy’s Special Function Library

SciPy’s library provides additional functionality for special mathematical functions, including Hermite polynomials through scipy.special.hermite. While it doesn’t directly compute the roots, you can use it in combination with optimization routines such as scipy.optimize.newton() to find them iteratively.

Here’s an example:

import numpy as np
from scipy.special import hermite
from scipy.optimize import newton

# Define Hermite polynomial of degree 4
H = hermite(4)

# Function to find roots for
func = lambda x: H(x)

# Guess and compute roots
guesses = [1j, -1j, 1 + 1j, -1 - 1j]
roots = [newton(func, x0=guess) for guess in guesses]
print(roots)

Output:

[(-0.7071067811865476+0.7071067811865475j),
 (-0.7071067811865476-0.7071067811865475j),
 (0.7071067811865476+0.7071067811865475j),
 (0.7071067811865476-0.7071067811865475j)]

In this method, the scipy.special.hermite function creates a Hermite polynomial, and the scipy.optimize.newton method iteratively finds the roots using initial guesses provided.

Method 3: Utilizing SymPy for Symbolic Mathematics

SymPy is a Python library for symbolic mathematics that allows for the exact representation and manipulation of mathematical expressions, including finding roots for Hermite polynomials symbolically. The sympy.polys.special.polynomials.hermite() function creates a Hermite polynomial, which can be solved using sympy.solvers.solve().

Here’s an example:

import sympy as sp

# Define variable and Hermite polynomial
x = sp.symbols('x')
poly = sp.hermite(4, x)

# Solve for roots
roots = sp.solve(poly, x)
print(roots)

Output:

[sqrt(2), -sqrt(2), sqrt(2)*I, -sqrt(2)*I]

With SymPy, a symbolic variable is created and a Hermite polynomial of a specific degree is defined with that variable. Subsequently, sp.solve() finds the roots, providing exact symbolic representations.

Method 4: Root Finding with mpmath

mpmath is a library for real and complex floating-point arithmetic with arbitrary precision. It includes a routine for computing the zeros of orthogonal polynomials, including Hermite polynomials, using mpmath.hermite() and mpmath.polyroots().

Here’s an example:

from mpmath import mp, hermite, polyroots

# Set precision
mp.dps = 15

# Define Hermite polynomial and compute roots
coeffs = hermite(4, polys=True)
roots = polyroots(coeffs)
print(roots)

Output:

[-1.224744871391589, 1.224744871391589, (-0.0 + 1.224744871391589j), (0.0 - 1.224744871391589j)]

Here, mpmath.hermite() is used to obtain the coefficients of the Hermite polynomial, which are then passed to mpmath.polyroots() to find the roots with high precision.

Bonus One-Liner Method 5: Quick Root Estimation with NumPy

For rapid estimates of Hermite polynomial roots where high precision is not required, you can quickly compute them using NumPy’s numpy.roots() function directly on the coefficients array.

Here’s an example:

import numpy as np

# Hermite polynomial coefficients
coeffs = [1, 0, -2, 0, 4]

# Compute roots
roots = np.roots(coeffs)
print(roots)

Output:

[-1.22474487+0.j  1.22474487+0.j  0.        +1.22474487j
  0.        -1.22474487j]

This approach skips any specialized polynomial libraries and directly uses the general-purpose np.roots() to estimate the roots from Hermite coefficients.

Summary/Discussion

  • Method 1: Using NumPy’s Polynomial Library. It is convenient due to its simplicity and integration in the widely-used NumPy library. However, it may not provide the highest precision for complex roots.
  • Method 2: Using SciPy’s Special Function Library. It offers a reliable iterative approach to finding roots but requires initial guesses which could be tricky to approximate.
  • Method 3: Utilizing SymPy for Symbolic Mathematics. SymPy offers exact symbolic solutions which are beneficial for theoretical analysis. Nevertheless, symbolic computation can be slower than numeric methods.
  • Method 4: Root Finding with mpmath. This is ideal for high precision calculations, but comes at the cost of performance due to the arbitrary precision arithmetic.
  • Bonus Method 5: Quick Root Estimation with NumPy. This is the simplest and fastest approach, though it may lack the precision and robustness of the methods tailored to Hermite polynomials.