Evaluating 2D Legendre Series with 1D Coefficients Arrays in Python

by

💡 Problem Formulation: Evaluating a 2D Legendre series involves a scenario where you have two sets of variables, usually x and y , and you want to compute the series sum based on the Cartesian product of these variables using a given set of coefficients stored in a 1D array. For example, given 1D coefficients array c , the goal is to find the value of the series at the points (x, y) produced by the Cartesian product of two separate ranges of x and y values.

Method 1: Using NumPy with Legendre Polynomial Functions

This method involves employing NumPy and the NumPy polynomial Legendre module, which provides a convenient set of functions to handle Legendre polynomials. By using numpy.polynomial.legendre.leggrid2d , you can easily evaluate the 2D Legendre series over a grid defined by the Cartesian product of x and y vectors based on the coefficients array.

Here’s an example:

import numpy as np
from numpy.polynomial.legendre import leggrid2d

# Define the 1D coefficients array
coeffs = np.array([1, 2, 3])

# Define the 1D arrays for x and y
x = np.linspace(-1, 1, 5)
y = np.linspace(-1, 1, 5)

# Evaluate the 2D Legendre series
result = leggrid2d(x, y, coeffs)

print(result)

The output will be the evaluated 2D Legendre series on the Cartesian product of x and y.

This snippet demonstrates creating a grid using np.linspace() for x and y values, then evaluating the legendre series with leggrid2d(). It’s concise and leverages NumPy’s efficiency for performance.

Method 2: Custom Implementation Using Recursive Legendre Polynomial Calculation

A more hands-on approach involves manually implementing the recursive definition of Legendre polynomials. This method might be attractive if you want full control over the calculation process or if you want to work without external libraries.

Here’s an example:

import numpy as np

def legendre_poly(n, x):
    if n == 0:
        return x*0 + 1
    elif n == 1:
        return x
    else:
        return ((2*n-1)*x*legendre_poly(n-1, x)-(n-1)*legendre_poly(n-2, x))/n

# Define coefficients array and x, y ranges
coeffs = [1, 2, 3]
x_vals = np.linspace(-1, 1, 5)
y_vals = np.linspace(-1, 1, 5)

# Compute the cartesian product of x and y
X, Y = np.meshgrid(x_vals, y_vals)
Z = np.zeros_like(X)

# Evaluate the legendre series for each order on the grid
for n, c in enumerate(coeffs):
    Z += c * np.outer(legendre_poly(n, x_vals), legendre_poly(n, y_vals))

print(Z)

The output is a 5×5 matrix representing the evaluated Legendre series on the grid.

Here, a custom function legendre_poly() implements the recursion for Legendre polynomials, which is then used for each coefficient in the array over the Cartesian product of x_vals and y_vals.

Method 3: Exploiting SymPy for Symbolical Legendre Polynomial Expansion

SymPy is a Python library for symbolic mathematics. It can find the Legendre polynomial expression symbolically, which can then be used for evaluation over any number of points. This method is beneficial when precision is a priority or when you wish to manipulate the expression algebraically before numeric evaluation.

Here’s an example:

from sympy import legendre, symbols, lambdify
import numpy as np

x, y = symbols('x y')
coeffs = [1, 2, 3]
X_vals = np.linspace(-1, 1, 5)
Y_vals = np.linspace(-1, 1, 5)

# Get a symbolic Legendre series sum
legendre_series_sum = sum(c * legendre(i, x) * legendre(i, y) for i, c in enumerate(coeffs))

# Lambdify the expression for fast numerical evaluation
func = lambdify((x, y), legendre_series_sum, 'numpy')

# Evaluate the series over a grid
X, Y = np.meshgrid(X_vals, Y_vals)
result = func(X, Y)
print(result)

The output will be the numeric evaluation of the symbolic Legendre series.

In the above code, legendre() generates the symbolic Legendre polynomial, lambdify() converts it to a function that can be used with NumPy arrays, and we evaluate this function over a grid.

Method 4: Using SciPy for Advanced Polynomial Operations

SciPy extends the capabilities of NumPy with additional modules for optimization, special functions, and polynomial manipulation. Similar to NumPy, you can use SciPy to evaluate Legendre polynomials but with additional functionality for handling more complex scenarios or customization.

Here’s an example:

from scipy.special import eval_legendre
import numpy as np

coeffs = [1, 2, 3]
x_vals = np.linspace(-1, 1, 5)
y_vals = np.linspace(-1, 1, 5)

# Compute the cartesian product of x and y
X, Y = np.meshgrid(x_vals, y_vals)
Z = np.zeros_like(X)

# Evaluate the Legendre series for each order on the grid
for n, c in enumerate(coeffs):
    Z += c * np.outer(eval_legendre(n, x_vals), eval_legendre(n, y_vals))

print(Z)

The output is a 5×5 matrix representing the evaluated Legendre series on the grid.

The SciPy function eval_legendre() is used to compute the polynomials and the results are aggregated using NumPy’s outer() function to form the final series on the grid.

Bonus One-Liner Method 5: Compact NumPy Implementation

If you’re looking for a concise, one-liner solution and you’re already familiar with NumPy’s broadcasting rules, you can evaluate the 2D Legendre series by combining several NumPy functions into a single expression.

Here’s an example:

import numpy as np
from numpy.polynomial.legendre import legval

coeffs = [1, 2, 3]
x = np.linspace(-1, 1, 5)
y = np.linspace(-1, 1, 5)

result = legval(x, coeffs)[:, None] * legval(y, coeffs)[None, :]
print(result)

The output will be a matrix with the results of the evaluated Legendre series on the grid.

This one-liner utilizes legval() function for evaluating Legendre polynomials directly, with NumPy’s broadcasting to align the dimensions and perform the Cartesian product implicitly.

Summary/Discussion

  • Method 1: NumPy with Legendre Polynomial Functions. It’s straightforward and leverages NumPy’s optimized code. The downside is the possibility of not matching the coefficient convention you might need.
  • Method 2: Custom Recursive Implementation. Offers maximum control over the evaluation process and does not rely on external libraries. However, it is not as optimized as library-based solutions.
  • Method 3: Using SymPy for Symbolic Computation. Provides high precision and algebraic manipulation before numeric evaluation. The drawback is slower performance due to symbolic computation overhead.
  • Method 4: Using SciPy for Advanced Polynomial Operations. Brings additional features to NumPy’s basic functionality, suitable for complex scenarios. It might be considered overkill for simple tasks.
  • Bonus Method 5: Compact NumPy One-Liner. Is elegantly concise, but relies on understanding of NumPy broadcasting rules and might obscure the code’s intent for less experienced users.