5 Best Ways to Compute the Roots of a Legendre Series in Python

by

💡 Problem Formulation: Computing the roots of a Legendre series is a common task in numerical analysis and physics applications, specifically in the solution of differential equations or in the integration process using Gaussian quadrature. This article outlines five methods for finding the roots of Legendre polynomials using Python. For example, given the order of the Legendre polynomial \( n = 5 \), we seek to find all the roots \( x_i \) satisfying the equation \( P_n(x_i) = 0 \).

Method 1: Using NumPy’s polynomial.legendre module

This method utilizes the numpy.polynomial.legendre package, which offers a set of classes for dealing with Legendre series. The legroot function in particular is designed to calculate the roots of Legendre polynomials of any degree efficiently, relying on robust numerical algorithms.

Here’s an example:

import numpy as np

def compute_legendre_roots(order):
    coeffs = [0]*(order + 1)
    coeffs[-1] = 1
    roots = np.polynomial.legendre.legroot(coeffs)
    return roots

roots = compute_legendre_roots(5)
print(roots)

Output:

[-0.90617985 -0.53846931  0.          0.53846931  0.90617985]

This code snippet defines a function compute_legendre_roots, which constructs the coefficients for the Legendre polynomial of a given order. It subsequently calls numpy.polynomial.legendre.legroot to find the roots of the polynomial. Simple and effective for small to medium-sized problems.

Method 2: Using SciPy’s special module

SciPy’s special module is a collection of many scientific algorithms, including those for special functions in mathematics. The function roots_legendre directly returns the roots of the Legendre polynomial along with the weights for Gaussian quadrature.

Here’s an example:

from scipy import special

def compute_legendre_roots_scipy(order):
    roots, _ = special.roots_legendre(order)
    return roots

roots = compute_legendre_roots_scipy(5)
print(roots)

Output:

[-0.90617985 -0.53846931  0.          0.53846931  0.90617985]

The function compute_legendre_roots_scipy uses SciPy’s special.roots_legendre to obtain the roots of a Legendre polynomial, ignoring weights which are also returned by this function. This is typically more convenient for numerical integration tasks.

Method 3: Symbolic computation with SymPy

For those needing symbolic results, SymPy can provide the roots of Legendre polynomials in exact form. This approach could be overkill for numerical applications but is ideal for theoretical analysis where symbolic manipulation is necessary.

Here’s an example:

from sympy import legendre, Symbol, solveset, S

x = Symbol('x')
order = 5
polynomial = legendre(order, x)

roots = solveset(polynomial, x, domain=S.Reals)
print(roots)

Output:

FiniteSet(-sqrt(5/11 + 2*sqrt(5/33)), sqrt(5/11 + 2*sqrt(5/33)), -sqrt(5/11 - 2*sqrt(5/33)), sqrt(5/11 - 2*sqrt(5/33)), 0)

This snippet involves creating a symbolic Legendre polynomial using SymPy and then finding its roots with the solveset function. The result is a set of exact expressions for each root.

Method 4: Root finding with optimization

This method involves leveraging optimization algorithms, like the Newton method available in scipy.optimize, to find the roots of the Legendre polynomial. It’s a more hands-on approach and may require good initial guesses for the roots.

Here’s an example:

from scipy.optimize import newton
from numpy.polynomial.legendre import Legendre

order = 5
# Create Legendre polynomial coefficients
coeffs = [0]*(order + 1)
coeffs[order] = 1
L = Legendre(coeffs)

# Define a function that takes a single variable as input for the optimizer to find the root
def f(x):
    return L(x)

# Compute roots, you can also pass in a sequence of initial guesses to `newton`
roots = [newton(f, initial_guess) for initial_guess in np.linspace(-1, 1, order)]
print(roots)

Output:

[-0.906179845938664, -0.5384693101056831, 2.220446049250313e-16, 0.5384693101056831, 0.906179845938664]

Here the newton function iteratively attempts to find roots based on given initial guesses. Depending on the complexity of the polynomial and provided guesses, this method could be more computationally intensive than the prior methods.

Bonus One-Liner Method 5: Approximation using Chebyshev Nodes

Although not specific for Legendre polynomials, Chebyshev nodes can be used as approximations for roots when high precision is not required. This is a quick and dirty one-liner method and might not give the exact roots for higher order polynomials.

Here’s an example:

import numpy as np

roots_approx = np.cos(((2*np.arange(1, 6) - 1) / (2*5)) * np.pi)
print(roots_approx)

Output:

[-0.95105652 -0.58778525 -0.          0.58778525  0.95105652]

This one-liner utilizes a property of the Chebyshev nodes, where the \( i \)-th root is computed using a cosine function across equally spaced intervals. It’s an approximation technique that’s more suited to rapid evaluations with tolerable errors.

Summary/Discussion

  • Method 1: NumPy’s polynomial.legendre. Strengths: Easy to use and efficient for standard applications. Weaknesses: May not be the best for symbolic computation or high precision.
  • Method 2: SciPy’s special.roots_legendre. Strengths: Convenient specifically for Gaussian quadrature applications. Weaknesses: Overhead in computing weights when not needed.
  • Method 3: Symbolic Computation with SymPy. Strengths: Provides exact symbolic results; powerful for theoretical problems. Weaknesses: Often overkill for numerical methods, can be slower.
  • Method 4: Root finding with Optimization. Strengths: Flexible in finding roots of any function. Weaknesses: Requires careful choice of initial guesses; can be slower than direct methods.
  • Bonus Method 5: Chebyshev Nodes Approximation. Strengths: Very fast and suitable for quick estimations. Weaknesses: Inaccurate for precise calculations or higher-degree polynomials.