高斯-勒让德求积公式：原理、Python实现与数值微分实战

做科学计算久了你会发现，微积分课上的求积分技巧到了工程现场基本靠不住——手里的被积函数通常来自实测数据，或者来自另一个数值程序，根本拿不出原函数。这个时候只能老老实实走数值积分。而在只拿到离散采样点、没有解析表达式的地方，数值微分又是绕不开的工具。这两件事里，高斯勒让德求积公式是我个人用得最频繁、也最愿意逢人就安利的一套方案：同样的节点数，它能把代数精度做到传统方法的近两倍，两个点就能精确积出三次多项式，四五个点就能把光滑函数的积分误差压到接近机器精度。

这篇文章不打算只给结论。我会把高斯-勒让德求积公式的原理掰开揉碎讲清楚，然后放出可以直接抄作业的 Python 实现，从手写牛顿迭代到与 SciPy 结果对拍，最后再聊聊数值微分的差商、Richardson 外推，以及一堆我实际踩过的坑。正在学数值分析的本科生、写有限元和谱方法的科研人员、被“求不出原函数的积分”折磨的工程师，应该都能从这里找到有用的东西。

1. 为什么数值积分绕不开高斯-勒让德

1.1 所有求积公式其实都在干同一件事

求积公式的外形高度统一：把积分 ∫_{a}^{b} f(x)dx 近似成一堆权重 w_i 乘以函数值 f(x_i) 的求和。梯形公式是取两个端点，权重各 (b-a)/2；辛普森公式是取三个等距点，权重按 1:4:1 配。两者本质上都是先用插值多项式替代被积函数：在给定节点上构造一个多项式穿过这些点的函数值，再对这个多项式精确求积，结果就是相应的求积公式。这类公式统称插值型求积公式。

这种视角最大的好处是，求积精度完全由“插值多项式在多大程度上逼近 f(x)”决定。等距节点的多项式插值有一个著名隐患：节点数增多时，龙格现象会让高次插值在区间端点附近剧烈振荡，误差不降反升。所以梯形、辛普森这类基于等距节点的公式，工程上通常以“复合”形式使用：把区间切成很多小段，每段用低阶公式。这也是后文精度对比里的主角。

提示：理解“积分 = 对离散点函数值的加权求和”这个框架，后面对高斯公式的理解会顺很多。权重是正的还是负的、节点落在哪里，直接决定数值稳定性。

1.2 代数精度才是衡量求积公式的核心标尺

数值分析里评价求积公式，最常用的指标是代数精度：如果某个公式对任意次数不超过 m 的多项式都精确，但对 m+1 次多项式做不到，就说它的代数精度是 m。为什么盯着多项式看？因为光滑函数在局部可以用泰勒展开近似成多项式，多项式精度上去了，光滑函数的误差自然跟着降。

梯形公式的代数精度是 1，辛普森是 3。n 个节点的插值型求积公式，代数精度一般在 n 左右，极少超过 n+1。而高斯-勒让德给出的答案是：n 个节点，代数精度 2n-1。两个点精确积三次多项式，三个点精确积五次多项式。这个“免费翻倍”是整套公式最值钱的地方，也是它能在有限元里成为默认积分规则的根本原因。

1.3 高斯的选择：把节点从“预定”改成“自由”

传统求积公式里，节点位置是预先定死的，能优化的只有权重。高斯-勒让德的想法很直接：把节点也当成未知数。n 个节点加 n 个权重，一共 2n 个自由度，正好可以用来让公式对 1, x, x², ..., x^{2n-1} 这 2n 个基函数全部精确，于是代数精度冲到 2n-1。

高斯当年证明了这是 n 点求积公式能达到的理论上限。想要达到这个上限，需要解一组非线性方程：让公式对前 2n 个单项式的积分误差全部为零。但实际工程里没人真的去硬解那组方程，因为数学上早就给出了漂亮的封闭解——节点恰好就是勒让德多项式 P_n(x) 的 n 个根。这个结论既是巧合也是必然，下一节拆开看。

2. 原理拆解：节点、权重与误差

2.1 勒让德多项式与三项递推

勒让德多项式是定义在 [-1,1] 上的一族正交多项式，前几项长这样：

• P_0(x) = 1
• P_1(x) = x
• P_2(x) = (3x² - 1) / 2
• P_3(x) = (5x³ - 3x) / 2
• P_4(x) = (35x⁴ - 30x² + 3) / 8

它们的核心性质是正交性：任意两个不同阶的勒让德多项式在 [-1,1] 上的积分为零，P_n 与自己的积分等于 2/(2n+1)。工程上不会真的展开系数去构造 P_n，而是用三项递推：

P_{k+1}(x) = ((2k+1) x P_k(x) - k P_{k-1}(x)) / (k+1)

这个递推数值上非常稳，因为正交多项式族的三项递推是良定的。为什么节点偏偏是勒让德多项式的根？可以粗糙地理解：正交性保证了这些节点在区间上“分布最不浪费”，两端密、中间疏，正好匹配多项式插值误差的分布规律。

2.2 节点是根，权重是闭式公式

只要节点是 P_n 的根，权重可以直接写成闭式：

w_i = 2 / [(1 - x_i²)(P'_n(x_i))²]

推导思路是：用高斯求积公式对拉格朗日插值基函数求积，再利用代数精度对它平方也精确的性质，最后化简成上面的形式。这个式子有个很实用的副产品：所有权重都为正。正权重意味着求积过程不会出现大数相减，舍入误差不会被放大，这是梯形和某些高阶牛顿-柯特斯公式做不到的。

我把 n=2 到 n=5 的节点和权重列出来，方便手工验证：

一个最简单的自查方法：所有权重之和必须等于 2，因为常数函数在 [-1,1] 上的积分就是区间长度 2。

2.3 从 [-1,1] 到 [a,b] 的变量替换

高斯-勒让德所有理论都建立在标准区间 [-1,1] 上，一般区间要先做线性替换：

x = (b-a)/2 · t + (a+b)/2

于是：

∫_{a}^{b} f(x)dx = (b-a)/2 · ∫_{-1}^{1} f((b-a)/2 · t + (a+b)/2) dt

实际代码里最常出的错就是忘了乘 (b-a)/2 这个因子。整个模块积分结果差一个常数倍，排查半天最后发现就死在这行，我见过不止一次。

如果被积函数在整段区间上不太光滑，比如某个子区间变化剧烈，工程上的标准做法是把 [a,b] 切成若干段，每段单独用高斯-勒让德，这就是复合高斯求积。实践中每段用 n=4 到 n=8 个点通常就够，继续加点的收益很有限，节点过度密集反而会引入舍入问题。

2.4 误差项：每加一个点，精度跳一档

高斯-勒让德在 [-1,1] 上的误差有一个精确表达式：

R_n = [2^{2n+1}(n!)⁴ / ((2n+1)((2n)!)³)] · f^{(2n)}(ξ)，ξ ∈ (-1,1)

换到一般区间还要乘 ((b-a)/2)^{2n+1}。这个系数随 n 增大衰减极其快：n=2 时约 1/135，n=3 时约 6.4e-5，n=4 时约 2.9e-7，n=5 时已经到 8e-10 量级。换句话说，只要被积函数足够光滑，每增加一个节点，误差主项大致缩小一到两个数量级。

这解释了数值实验里一个反直觉的现象：高斯求积用几十个点做到的事情，梯形公式往往要几万个点。当然前提是“足够光滑”。一旦被积函数带奇点，f^{(2n)}(ξ) 会暴涨，上面的漂亮结论就不成立，这时候要用第 5 节里的补救手段。

3. 代码实现：从牛顿迭代到 Golub-Welsch

3.1 手写牛顿迭代求根

先给一份“教学级”实现，用牛顿迭代求 P_n 的根。关键在于初始值要选得好。工程里实践证明，用 cos(π(i+0.75)/(n+0.5)) 作为第 i 个根的初始猜测非常稳——这套初始猜测本来是给切比雪夫节点用的，拿来启动勒让德求根也接近最优。

import numpy as np def legendre_poly(n, x): """三项递推计算 P_n(x)""" if n == 0: return 1.0 if n == 1: return x p_prev, p_curr = 1.0, x for k in range(1, n): p_next = ((2 * k + 1) * x * p_curr - k * p_prev) / (k + 1) p_prev, p_curr = p_curr, p_next return p_curr def legendre_deriv(n, x): """P_n'(x) = n (x P_n(x) - P_{n-1}(x)) / (x^2 - 1)""" if n == 0: return 0.0 p = legendre_poly(n, x) p_prev = legendre_poly(n - 1, x) return n * (x * p - p_prev) / (x * x - 1.0) def gauss_legendre_newton(n): """牛顿迭代求节点，闭式公式求权重""" nodes, weights = [], [] for i in range(n): x = np.cos(np.pi * (i + 0.75) / (n + 0.5)) # 初始猜测 for _ in range(50): corr = legendre_poly(n, x) / legendre_deriv(n, x) x -= corr if abs(corr) < 1e-15: break w = 2.0 / ((1.0 - x * x) * legendre_deriv(n, x) ** 2) nodes.append(x) weights.append(w) return np.array(nodes), np.array(weights)

实测下来 n=10 以内，这套手写结果和 NumPy 的leggauss相差在 1e-14 量级，足够验证原理。但 n 超过 30 就开始吃力——高阶勒让德多项式的根对扰动越来越敏感，牛顿迭代的精度会一点点流失。这时候该换生产级方案。

3.2 生产级方案：Golub-Welsch 特征值法

Golub-Welsch 算法是计算高斯求积节点和权重的标准做法。结论很漂亮：构造一个 n×n 对称三对角雅可比矩阵，其非零副对角线元素取

β_k = k / √(4k² - 1), k = 1, 2, ..., n-1

那么矩阵的特征值就是高斯节点，权重由第一个分量的特征向量平方给出：w_i = 2(v_{0,i})²。这个做法把求根问题变成特征值问题，对称三对角矩阵的特征值求解极其稳定，n 到几百都没问题。NumPy 的leggauss内部走的也是类似路线。

def gauss_legendre_golubwelsch(n): """Golub-Welsch：节点是 Jacobi 矩阵特征值，权重来自特征向量""" k = np.arange(1, n, dtype=float) beta = k / np.sqrt(4.0 * k * k - 1.0) J = np.diag(beta, 1) + np.diag(beta, -1) evals, evecs = np.linalg.eigh(J) nodes = evals weights = 2.0 * evecs[0, :] ** 2 return nodes, weights

以 n=2 验证：β = 1/√3 ≈ 0.57735，矩阵 [[0,β],[β,0]] 的特征值正是 ±0.57735，两个特征向量第一分量平方都是 1/2，乘 2 得权重 1。和熟知的二点高斯公式完全一致。

3.3 积分器封装与三种方法精度实测

有了节点和权重，积分器封装就三行：

from numpy.polynomial.legendre import leggauss def gauss_legendre_integral(f, a, b, n): x, w = leggauss(n) t = 0.5 * (b - a) * x + 0.5 * (a + b) return 0.5 * (b - a) * np.dot(w, f(t))

拿 ∫_{0}^{1} e^{-x} dx = 1 - 1/e ≈ 0.6321205588 做实测，结果如下表：

再测一个多项式：∫_{0}^{1} x⁴ dx = 0.2。高斯-勒让德 n=3 的代数精度是 5，结果是 0.2000000000000000，一个点不差；复合梯形想达到这个精度，步长得压到千分之一以下。这就是“代数精度翻倍”最直白的证据。

注意：积分结果对节点顺序不敏感，leggauss返回升序、手写代码返回降序，变换到 [a,b] 后不影响权重和。但如果你在调试，建议先固定用leggauss的结果做基准，别再叠加其他变量。

4. 数值微分：差商公式与误差控制

4.1 从导数定义出发的三套差商

数值微分最朴素的想法就是回到导数定义。常用的三套公式：

中心差商虽然要多算一个点，但精度从 O(h) 直接提到 O(h²)，性价比极高。我几乎不用一阶差商，除非数据真的只有单侧可用，比如边界点。

4.2 Richardson 外推：用两次计算换一个精度量级

中心差商的误差是 c·h² + O(h⁴)，那用 h 和 h/2 各算一次，线性组合消掉 h² 项：

D_extrap = (4·D(h/2) - D(h)) / 3

精度直接从 O(h²) 变成 O(h⁴)。这个技巧叫 Richardson 外推，本质是用两次廉价计算换一个精度量级。实测 f(x) = sin(x) 在 x = π/4 的导数：h=0.1 时中心差商误差约 1.2e-3，h=0.05 时约 3e-4，外推后约 7e-6；把 h 再缩小一个量级，外推结果能到 1e-9 附近。

外推还能继续叠。再用 h/4 算一次，又可以消掉 h⁴ 项，这就是 Romberg 思想的数值微分版本，和 Romberg 积分同出一套血统。收敛速度确实快，但每次外推都要额外算函数值，实际用到二阶外推就差不多了。

4.3 二阶导数、步长选择和舍入误差

二阶导数的中心差商是 (f(x+h) - 2f(x) + f(x-h)) / h²。这里有个反直觉的坑：和积分不同，数值微分天然病态。步长 h 太大，截断误差压不住；h 太小，f(x+h) 和 f(x-h) 非常接近，两个大数相减把有效数字洗掉了，舍入误差开始统治结果。

浮点下有个经验最优步长：中心一阶差商的最优 h 大约在 eps^(1/3) 量级，float64 下约 6e-6；二阶导数的最优 h 大约在 eps^(1/4) 量级，约 1e-4。所以千万别一味把 h 往小调。实测 h 从 1e-5 一路压到 1e-10，误差会先降后涨，第一次遇到这种现象的人多半都会被吓一跳。

如果你能拿到函数源码而不是黑盒，现代做法是用自动微分（JAX、PyTorch 的 grad），链式法则能把导数算到接近机器精度，完全不踩差商的坑。数值微分真正的用武之地是黑盒函数、实验数据，以及嵌入在大型数值流程里的局部求导。

提示：高斯-勒让德的节点在高阶谱方法里还有第二重身份——在勒让德-高斯-洛巴托节点上构造的谱微分矩阵，本质是“全局版数值微分”，流体力学的谱元法、地震波场模拟都在用。同一个“高斯点即最优离散点”的思想贯穿了积分和微分两条线。

5. 常见问题与避坑经验

5.1 我见过的三处经典翻车点

第一，忘了区间缩放因子 (b-a)/2。这是在代码审查里抓到过最多的错误，整个模块的积分结果差一个常数倍，写的时候少乘一个因子，排查能浪费一整天。

第二，节点权重自查对不上。权重和不是 2，或者 ∫_{-1}^{1} x² dx 用 n=2 对不上 2/3，都是实现出了问题的明显信号。我习惯在写完求积函数后立刻跑这几个自测，而不是直接扔进大工程里。

第三，拿高斯公式硬套离散采样数据。高斯节点是固定的，和你的采样网格几乎不可能重合，硬套只能先插值，插值又引入自己的误差。这种场景老实跑复合梯形或复合辛普森，或者先做样条插值再积分，别迷信高斯的名头。

5.2 端点奇异和振荡积分怎么办

端点奇异是高斯求积最大的克星。比如 ∫_{0}^{1} ln(x) dx，对数函数在 0 处发散，高斯-勒让德硬算收敛速度从指数掉回代数级，非常浪费。常用三个办法：奇异点减去法，把 ln(x) 拆成可积部分加光滑部分；变量替换，比如 x = t² 把奇点拉开；或者换带权的高斯-雅可比公式。我优先用变量替换，代码改动最小，效果也稳定。

振荡积分是另一个场景。∫ sin(ωx)·g(x) dx 在 ω 很大时，高斯求积需要的点数大致正比于 ω，点数一多，计算成本直线上去。这时 Filon 或 Levin 方法更划算。一个低成本技巧是分段：在每个振荡半周期内布点，保证函数在小区间里单调，高斯求积的精度会明显恢复。

5.3 工程里的实战建议

在有限元里，高斯积分是单元刚度矩阵的默认积分规则。一维单元用 1 到 4 个积分点很常见；对高阶单元或者形函数复杂的情况，积分点跟着阶数走。一个长期有效的经验：被积函数越光滑，越值得用高斯公式；越不光滑，越不值得。工程上九成积分问题用每段 4 到 8 个积分点就能解决，盲上几十个点往往只是 CPU 白烧。

高维积分是另一个大坑。张量积形式的 d 维高斯求积要 n^d 个点，d=6、n=4 就是 4096 个点，每个点还要算一次被积函数，成本瞬间爆炸。我的经验是 d≤3 用张量积高斯，d≥4 直接切蒙特卡洛或者稀疏网格。别迷信“高斯精度高”——维度上去了，高斯也救不了你。

我在实际项目中养成的习惯是：把常用阶数的勒让德节点和权重以常量表写死在库里，而不是每次调用时重新算一遍。省掉的那点时间是小事，关键是排除了初始化代码出错的隐患。还有个小技巧值得分享：不确定该用多少点时，就跑 n 点和 n+1 点的结果，两个结果差小于容差就收手，这就是自适应高斯求积的雏形，也是我调试数值代码时最常用的一个判断手段。