在Python中使用3D系数数组对x和y的笛卡尔积上评估2D Legendre级数

简介

笛卡尔坐标系的坐标点可以用一对有序数对 $(x,y)$ 表示。对于每个坐标点 $(x,y)$ ，我们可以使用一个函数 $f(x,y)$ 对其进行描述。在多项式拟合中，常常使用的是二元多项式（即二次方程）描述 $f(x,y)$ ，而Legendre系数由于其良好的性质，也常被用于二元多项式拟合。

Legendre级数

对于任意的一元函数 $f(x)$ ，在 $[-1,1]$ 区间上，我们可以使用Legendre多项式进行拟合，即
$f(x)=\sum_{n=0}^{\infty}a_nP_n(x)$
其中 $a_n$ 表示 $f(x)$ 在Legendre多项式 $P_n(x)$ 下的展开系数，通常我们只需选取一部分系数，即可对 $f(x)$ 进行较好的拟合。

对于二元函数 $f(x,y)$ ，我们也可以进行类似的展开。我们用Legendre多项式 $P_n(x)$ 和 $P_m(y)$ 的笛卡尔积得到 $L_{nm}(x,y)=P_n(x)P_m(y)$ ，因此有
$f(x,y)=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}a_{nm}L_{nm}(x,y)$
其中 $a_{nm}$ 仍然表示在 $L_{nm}(x,y)$ 下 $f(x,y)$ 的展开系数。同样地，我们通常只需选取一部分系数，即可对 $f(x,y)$ 进行较好的拟合。

对于Legendre多项式 $P_n(x)$ 和 $P_m(y)$ ，可以使用scipy.special库中的函数eval_legendre进行计算。

我们先定义一个求取Legendre级数的函数eval_legendre_2d（示例代码如下，为Python语言）。

from scipy.special import eval_legendre

def eval_legendre_2d(x, y, a_nm, n_max, epsilon=1e-5):
    """
    评估2D Legendre级数

    Parameters
    ----------
    x : float or ndarray
        x坐标
    y : float or ndarray
        y坐标
    a_nm : ndarray
        Legendre系数数组
    n_max : int
        最高展开阶数
    epsilon : float
        计算精度

    Returns
    -------
    float or ndarray
        笛卡尔积(x, y)上的2D Legendre级数值
    """
    ret = 0.0
    for n in range(n_max + 1):
        for m in range(n + 1):
            L_nm = eval_legendre(n, x) * eval_legendre(m, y)
            if abs(L_nm) < epsilon:
                continue
            a_nm_nm = a_nm[n][m] if (n < len(a_nm) and m < len(a_nm[n])) else 0.0
            ret += a_nm_nm * L_nm
    return ret

其中，n_max为最高展开阶数，a_nm为Legendre系数数组，epsilon为计算精度。该函数使用双重循环求解2D Legendre级数，第一重循环枚举 $n$ ，第二重循环枚举 $m$ 。对于每个阶数对应的Legendre多项式 $L_{nm}(x,y)$ ，计算当前坐标 $(x,y)$ 的函数值。

此外，程序中还使用了适当的优化，对于某些趋近于0的 $L_{nm}(x,y)$ 项，不进行计算，直接跳过。

下面，我们给出一个例子，使用eval_legendre_2d函数对二元函数 $f(x,y)=(x+2y+1)^2$ 进行拟合。

import numpy as np# 程序示例

首先，我们需要生成数据集。此处选取x, y的笛卡尔积，作为我们的数据点。

```python
x = np.linspace(-1, 1, 11)
y = np.linspace(-1, 1, 11)
xx, yy = np.meshgrid(x, y)
data = np.stack([xx.ravel(), yy.ravel()], axis=1)

接着，我们计算数据点上的函数值 $f(x,y)$ 。

f = (data[:, 0] + 2 * data[:, 1] + 1) ** 2

f中存储了数据点的函数值。接下来，我们需要求解Legendre系数 $a_{nm}$ 。

n_max = 5

# 计算Legendre系数
a_nm = np.zeros([n_max + 1, n_max + 1])
for n in range(n_max + 1):
    for m in range(n + 1):
        x_2n1 = 2 * np.arange(n + 1) + 1
        y_2m1 = 2 * np.arange(m + 1) + 1
        p = x_2n1.reshape(-1, 1) * yy.reshape(1, -1)
        q = xx.reshape(-1, 1) * y_2m1.reshape(1, -1)
        a_nm[n][m] = np.vdot(p, q * f) * (n + 0.5) * (m + 0.5) * ((-1) ** (n - m))
        if m != n:
            a_nm[m][n] = a_nm[n][m]

可以看到，我们使用了另一种方法，直接求取Legendre系数。相对于其他常见的方法（如基于广义逆的方法、最小二乘法等），这种方法在计算效率上更为优秀。

最后，我们使用eval_legendre_2d求解Legendre级数，并绘制出结果。

# x, y网格
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
xx, yy = np.meshgrid(x, y)

# 计算Legendre级数
z = eval_legendre_2d(xx, yy, a_nm, n_max)

# 绘制结果
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xx, yy, z.reshape(xx.shape))
plt.show()

另外，由于Jupyter Notebook支持多种语言内核，我们在实现中也可以使用其他语言来加速计算。例如，以下是使用numba库加速计算的代码示例（示例代码为Python语言）。

from numba import njit

@njit
def eval_legendre_2d_njit(x, y, a_nm, n_max, epsilon=1e-5):
    ret = 0.0
    for n in range(n_max + 1):
        for m in range(n + 1):
            L_nm = 1.0
            pn = 1.0
            for i in range(n):
                pn *= x
                L_nm *= (2 * i + 1) / (i + 1) * (1 - x**2)
            pn *= (2 * n + 1)
            L_nm *= (-1)**m * (pn**0.5) * (1 - y**2)**0.5
            if abs(L_nm) < epsilon:
                continue
            a_nm_nm = a_nm[n][m] if (n < len(a_nm) and m < len(a_nm[n])) else 0.0
            ret += a_nm_nm * L_nm
    return ret

由于numba库可以自动将Python代码编译成机器码，因此使用njit修饰函数eval_legendre_2d_njit即可获得一定的计算速度提升。对于一些较大规模的计算，这种方法非常有效。