正交多项式实战：用Python手写Legendre多项式拟合曲线（附代码）

最新推荐文章于 2026-06-27 20:00:00 发布

原创

最新推荐文章于 2026-06-27 20:00:00 发布 · 1k 阅读

正交多项式实战：用Python手写Legendre多项式拟合曲线（附代码）

在工程仿真、科学计算和数据分析的日常工作中，我们常常需要用一个简洁的数学模型去描述一组看似杂乱无章的观测数据，或者去逼近一个计算成本高昂的复杂函数。这个过程，我们称之为“函数逼近”或“曲线拟合”。传统的最小二乘法大家都不陌生，但当你尝试拟合一个高阶多项式时，可能会发现那个需要求解的“法方程”矩阵变得异常庞大且病态，计算不仅缓慢，结果还可能极不稳定。这时，一个在理论数学中熠熠生辉的工具——正交多项式——就能从幕后走到台前，成为解决这类问题的利器。

特别是Legendre多项式，它凭借在区间 [-1, 1] 上权函数为1（即“等权”）的优良特性，成为了离散数据最小二乘拟合中的“明星”基底。它最大的魔力在于，能将原本稠密、可能病态的法方程矩阵，转化为一个纯粹的对角阵。这意味着，拟合系数的求解从解一个复杂的线性方程组，退化成了一个个独立的、简单的内积计算。对于需要快速迭代、验证不同模型阶数的开发者而言，这无疑是效率的飞跃。

本文将从工程实践的角度出发，避开繁琐的理论推导，直接聚焦于如何用Python从零开始实现基于Legendre多项式的曲线拟合。我们将深入代码细节，探讨数据规范化、递推公式的稳定实现、以及如何处理非标准区间等实际问题，让你能立刻将这套方法应用到自己的项目中。

1. 正交多项式：为何是拟合的“降维打击”？

在深入代码之前，我们有必要先理解正交多项式为何能简化最小二乘拟合。假设我们有一组数据点 (x_i, y_i), i=1,...,m，想用一个n阶多项式 S(x) = a_0 + a_1*x + ... + a_n*x^n 去拟合。经典最小二乘法的目标是让残差平方和最小，这导出了一个关于系数 a_k 的线性方程组（法方程）：

(φ_0, φ_0) a_0 + (φ_0, φ_1) a_1 + ... + (φ_0, φ_n) a_n = (y, φ_0)
(φ_1, φ_0) a_0 + (φ_1, φ_1) a_1 + ... + (φ_1, φ_n) a_n = (y, φ_1)
...
(φ_n, φ_0) a_0 + (φ_n, φ_1) a_1 + ... + (φ_n, φ_n) a_n = (y, φ_n)

其中，φ_k(x) = x^k，而内积 (f, g) = Σ_i f(x_i) g(x_i)（对于离散数据）。这个方程组的系数矩阵（Gram矩阵）通常是一个稠密的、随着n增大条件数急剧恶化的矩阵，求解它既耗时又不稳定。

现在，如果我们选择一组在给定数据点集上离散正交的多项式序列 {P_0(x), P_1(x), ..., P_n(x)} 作为新的基底，即它们满足： (P_j, P_k) = 0, 当 j ≠ k。

那么，上述法方程矩阵中所有的非对角线元素都变成了0！方程组瞬间退化为n+1个独立的方程：

(P_k, P_k) * a_k = (y, P_k), for k = 0, 1, ..., n

每个系数 a_k 可以直接、独立地计算出来：a_k = (y, P_k) / (P_k, P_k)。计算复杂度从 O(n³) 量级（求解线性方程组）降低到了 O(n²) 量级（计算内积），并且完全避免了矩阵求逆带来的数值不稳定问题。这就是正交多项式在拟合中的核心优势。

提示：这里的“离散正交”内积定义依赖于具体的数据点集。而像Legendre多项式这样的经典正交多项式族，是在连续区间 [-1,1] 上、关于权函数 ρ(x)=1 连续正交的。当我们把数据点规范映射到 [-1,1] 区间后，可以近似利用其正交性，或者通过Gram-Schmidt过程在数据点上构造一组严格离散正交的多项式。

2. Legendre多项式：定义、性质与递推计算

Legendre多项式 P_n(x) 是定义在区间 [-1, 1] 上，关于权函数 ρ(x) = 1 正交的一组多项式序列。其正交性定义为：

∫_{-1}^{1} P_m(x) P_n(x) dx = 0, 当 m ≠ n
∫_{-1}^{1} [P_n(x)]^2 dx = 2 / (2n + 1)

最常用的生成方法是三项递推公式，它允许我们高效、稳定地计算任意高阶的Legendre多项式值：

P_0(x) = 1
P_1(x) = x
P_{n+1}(x) = [(2n+1) * x * P_n(x) - n * P_{n-1}(x)] / (n+1), 对于 n ≥ 1

这个递推关系是数值计算中的基石。下面，我们用Python函数来实现它，该函数能同时计算从0阶到n阶的所有Legendre多项式在给定x点处的值。

import numpy as np

def legendre_polynomials(x, n):
    """
    计算点x处从0阶到n阶的Legendre多项式值。

    参数
    ----------
    x : float 或 np.ndarray
        自变量值，可以是一个标量或数组。
    n : int
        所需计算的多项式最高阶数。

    返回
    -------
    P : np.ndarray
        形状为 (n+1, ...) 的数组。P[k] 对应 k 阶Legendre多项式在x处的值。
        如果x是标量，则返回形状为 (n+1,) 的一维数组。
        如果x是数组，则返回形状为 (n+1, x.shape) 的数组。
    """
    # 将输入转换为至少1维的数组以便于广播
    x = np.asarray(x)
    original_shape = x.shape
    x = x.ravel()  # 展平以便操作
    m = len(x)

    # 初始化结果数组
    P = np.zeros((n + 1, m))

    # 0阶和1阶多项式
    P[0, :] = 1.0
    if n >= 1:
        P[1, :] = x

    # 使用三项递推公式计算高阶项
    for k in range(1, n):
        # P_{k+1} = [(2k+1)*x*P_k - k*P_{k-1}] / (k+1)
        P[k + 1, :] = ((2.0 * k + 1.0) * x * P[k, :] - k * P[k - 1, :]) / (k + 1.0)

    # 恢复原始形状（如果x是多维的）
    if len(original_shape) > 1:
        new_shape = (n + 1,) + original_shape
        P = P.reshape(new_shape)
    elif original_shape == ():  # x是标量
        P = P[:, 0]  # 压扁为 (n+1,) 的一维数组

    return P

这个函数有几个值得注意的工程细节：