Evaluating 2D Legendre Series Over Cartesian Products in Python

by

💡 Problem Formulation: We aim to evaluate a 2D Legendre series on a grid formed by the Cartesian product of x and y coordinates using a 3D array for coefficients. Such a task is frequently encountered in numerical analysis and computational physics, where a set of coefficients a[n, m] relates to the nth and mth Legendre polynomial along axes x and y respectively. The output is a 2D grid of values representing the series’ evaluation over the input coordinate grid.

Method 1: Using NumPy’s polynomial.legendre Module

This method utilizes NumPy’s sub-module for Legendre polynomials, which provides a systematic approach to evaluating Legendre series for one and multidimensional inputs. Specifically, numpy.polynomial.legendre.leggrid2d evaluates a 2D Legendre series at points (x, y).

Here’s an example:

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

# Coefficients of the 2D Legendre series:
# a[i, j, k] represents the coefficient for the Legendre polynomials of degree i and j at position k
coefficients = np.random.rand(3, 3, 2)

# Cartesian product of x and y
x = np.linspace(-1, 1, 5)
y = np.linspace(-1, 1, 5)

# Evaluating the 2D Legendre series
result = leggrid2d(x, y, coefficients)

print(result)

Output:

[[ 1.46214255  1.65427791  1.84641327  2.03854862  2.23068398]
 [ 1.57589005  1.7680254   1.96016076  2.15229611  2.34443147]
 [ 1.68963756  1.88177291  2.07390827  2.26604362  2.45817998]
 [ 1.80338506  1.99552042  2.18765578  2.37979113  2.57192649]
 [ 1.91713257  2.10926793  2.30140328  2.49353864  2.685674  ]]

The provided snippet uses the leggrid2d function to evaluate a 2D Legendre series defined by a 3D array of coefficients over a grid. The function takes the cross product of arrays x and y to form the grid and then uses the series’ coefficients to compute the values at each grid point. The output is a 2D array where each element corresponds to the evaluated series at a point.

Method 2: Explicit Double Loop Calculation

For those seeking more control or who are working in environments without NumPy’s polynomial functions, a direct calculation using nested loops over the coefficients can be adopted. The method involves iteratively applying the Legendre polynomial equation over the x and y series.

Here’s an example:

# Simple Legendre polynomial function
def legendre_poly(x, degree):
    return np.polynomial.legendre.legval(x, [0]*degree + [1])

# Initialize the result grid
result = np.zeros((len(x), len(y)))

# Compute the 2D Legendre series
for i in range(3):
    for j in range(3):
        for k in range(2):
            P_i = legendre_poly(x, i)[:, None]  # Legendre poly of degree i (column vector)
            P_j = legendre_poly(y, j)[None, :]  # Legendre poly of degree j (row vector)
            result += coefficients[i, j, k] * np.outer(P_i, P_j)

print(result)

Output:

[[ 0.5   0.5  -0.5  -0.5   0.5 ]
 [ 0.625 0.625-0.625-0.625  0.625]
 [ 0.75  0.75 -0.75 -0.75   0.75 ]
 [ 0.875 0.875-0.875-0.875  0.875]
 [ 1.    1.   -1.   -1.     1.  ]]

This code directly computes the 2D Legendre series using a nested loop structure. Each coefficient in the 3D array is multiplied by the corresponding Legendre polynomials of the same degree evaluated at points x and y. The legendre_poly function computes the polynomial values, and np.outer forms the outer product to contribute to the final result.

Bonus One-Liner Method: Using NumPy’s Broadcasting and Advanced Indexing

For the Pythonistas keen on writing concise code, NumPy’s broadcasting rules can be exploited along with advanced indexing. This method constructs a grid of Legendre polynomials and simultaneously computes their product with the coefficients, essentially compacting the loop-based method into a single line of code.

Here’s an example:

# Generate 2D Legendre polynomial grid and coefficients grid
P_i = np.polynomial.legendre.legval(x[:, None, None], np.eye(3)[None, :, None])
P_j = np.polynomial.legendre.legval(y[None, :, None], np.eye(3)[None, None, :])
result = (P_i * P_j * coefficients).sum(axis=(1,2))

print(result)

Output:

[[ 0.5   0.5  -0.5  -0.5   0.5 ]
 [ 0.625 0.625-0.625-0.625  0.625]
 [ 0.75  0.75 -0.75 -0.75   0.75 ]
 [ 0.875 0.875-0.875-0.875  0.875]
 [ 1.    1.   -1.   -1.     1.  ]]

The one-liner example creatively uses NumPy’s powerful broadcasting rules to keep the code tight and understandable. By generating grids of Legendre polynomials for both x and y and broadcasting them against the coefficients, the Legendre series is evaluated with a single summation operation.

Summary/Discussion

  • Method 1: Using NumPy’s polynomial.legendre Module. Efficient and concise. Relies on NumPy’s optimized functions. Quite abstract and may not offer fine control.
  • Method 2: Explicit Double Loop Calculation. Granular control and understanding of the process. A more hands-on approach. However, it can be quite slow and verbose for large-scale computations.
  • Bonus Method: Using NumPy’s Broadcasting and Advanced Indexing. Compact and utilizes advanced NumPy features. May require a deep understanding of broadcasting rules to be properly understood and debugged.