高斯-勒让德求积公式：原理、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 的多项式都精确但对 m1 次多项式做不到就说它的代数精度是 m。为什么盯着多项式看因为光滑函数在局部可以用泰勒展开近似成多项式多项式精度上去了光滑函数的误差自然跟着降。梯形公式的代数精度是 1辛普森是 3。n 个节点的插值型求积公式代数精度一般在 n 左右极少超过 n1。而高斯-勒让德给出的答案是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) 1P_1(x) xP_2(x) (3x² - 1) / 2P_3(x) (5x³ - 3x) / 2P_4(x) (35x⁴ - 30x² 3) / 8它们的核心性质是正交性任意两个不同阶的勒让德多项式在 [-1,1] 上的积分为零P_n 与自己的积分等于 2/(2n1)。工程上不会真的展开系数去构造 P_n而是用三项递推P_{k1}(x) ((2k1) x P_k(x) - k P_{k-1}(x)) / (k1)这个递推数值上非常稳因为正交多项式族的三项递推是良定的。为什么节点偏偏是勒让德多项式的根可以粗糙地理解正交性保证了这些节点在区间上“分布最不浪费”两端密、中间疏正好匹配多项式插值误差的分布规律。2.2 节点是根权重是闭式公式只要节点是 P_n 的根权重可以直接写成闭式w_i 2 / [(1 - x_i²)(P_n(x_i))²]推导思路是用高斯求积公式对拉格朗日插值基函数求积再利用代数精度对它平方也精确的性质最后化简成上面的形式。这个式子有个很实用的副产品所有权重都为正。正权重意味着求积过程不会出现大数相减舍入误差不会被放大这是梯形和某些高阶牛顿-柯特斯公式做不到的。我把 n2 到 n5 的节点和权重列出来方便手工验证n节点 x_i权重 w_i2±0.57735026921.000000000030, ±0.77459666920.8888888889, 0.55555555564±0.3399810436, ±0.86113631160.6521451549, 0.347854845150, ±0.5384693101, ±0.90617984590.5688888889, 0.4786286705, 0.2369268851一个最简单的自查方法所有权重之和必须等于 2因为常数函数在 [-1,1] 上的积分就是区间长度 2。2.3 从 [-1,1] 到 [a,b] 的变量替换高斯-勒让德所有理论都建立在标准区间 [-1,1] 上一般区间要先做线性替换x (b-a)/2 · t (ab)/2于是∫_{a}^{b} f(x)dx (b-a)/2 · ∫_{-1}^{1} f((b-a)/2 · t (ab)/2) dt实际代码里最常出的错就是忘了乘 (b-a)/2 这个因子。整个模块积分结果差一个常数倍排查半天最后发现就死在这行我见过不止一次。如果被积函数在整段区间上不太光滑比如某个子区间变化剧烈工程上的标准做法是把 [a,b] 切成若干段每段单独用高斯-勒让德这就是复合高斯求积。实践中每段用 n4 到 n8 个点通常就够继续加点的收益很有限节点过度密集反而会引入舍入问题。2.4 误差项每加一个点精度跳一档高斯-勒让德在 [-1,1] 上的误差有一个精确表达式R_n [2^{2n1}(n!)⁴ / ((2n1)((2n)!)³)] · f^{(2n)}(ξ)ξ ∈ (-1,1)换到一般区间还要乘 ((b-a)/2)^{2n1}。这个系数随 n 增大衰减极其快n2 时约 1/135n3 时约 6.4e-5n4 时约 2.9e-7n5 时已经到 8e-10 量级。换句话说只要被积函数足够光滑每增加一个节点误差主项大致缩小一到两个数量级。这解释了数值实验里一个反直觉的现象高斯求积用几十个点做到的事情梯形公式往往要几万个点。当然前提是“足够光滑”。一旦被积函数带奇点f^{(2n)}(ξ) 会暴涨上面的漂亮结论就不成立这时候要用第 5 节里的补救手段。3. 代码实现从牛顿迭代到 Golub-Welsch3.1 手写牛顿迭代求根先给一份“教学级”实现用牛顿迭代求 P_n 的根。关键在于初始值要选得好。工程里实践证明用 cos(π(i0.75)/(n0.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)实测下来 n10 以内这套手写结果和 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, dtypefloat) 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以 n2 验证β 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 做实测结果如下表方法采样点数计算结果绝对误差复合梯形 h1/450.6354143.3e-3复合辛普森 h1/450.6321341.4e-5高斯-勒让德 n330.6321202033.6e-7高斯-勒让德 n550.6321205591e-12 量级再测一个多项式∫_{0}^{1} x⁴ dx 0.2。高斯-勒让德 n3 的代数精度是 5结果是 0.2000000000000000一个点不差复合梯形想达到这个精度步长得压到千分之一以下。这就是“代数精度翻倍”最直白的证据。注意积分结果对节点顺序不敏感leggauss返回升序、手写代码返回降序变换到 [a,b] 后不影响权重和。但如果你在调试建议先固定用leggauss的结果做基准别再叠加其他变量。4. 数值微分差商公式与误差控制4.1 从导数定义出发的三套差商数值微分最朴素的想法就是回到导数定义。常用的三套公式公式名称表达式截断误差主项前向差商(f(xh) - f(x)) / h(h/2) f(ξ)后向差商(f(x) - f(x-h)) / h-(h/2) f(ξ)中心差商(f(xh) - f(x-h)) / (2h)(h²/6) f(ξ)二阶中心差商(f(xh) - 2f(x) f(x-h)) / h²(h²/12) f⁗(ξ)中心差商虽然要多算一个点但精度从 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 的导数h0.1 时中心差商误差约 1.2e-3h0.05 时约 3e-4外推后约 7e-6把 h 再缩小一个量级外推结果能到 1e-9 附近。外推还能继续叠。再用 h/4 算一次又可以消掉 h⁴ 项这就是 Romberg 思想的数值微分版本和 Romberg 积分同出一套血统。收敛速度确实快但每次外推都要额外算函数值实际用到二阶外推就差不多了。4.3 二阶导数、步长选择和舍入误差二阶导数的中心差商是 (f(xh) - 2f(x) f(x-h)) / h²。这里有个反直觉的坑和积分不同数值微分天然病态。步长 h 太大截断误差压不住h 太小f(xh) 和 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 用 n2 对不上 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 个点d6、n4 就是 4096 个点每个点还要算一次被积函数成本瞬间爆炸。我的经验是 d≤3 用张量积高斯d≥4 直接切蒙特卡洛或者稀疏网格。别迷信“高斯精度高”——维度上去了高斯也救不了你。我在实际项目中养成的习惯是把常用阶数的勒让德节点和权重以常量表写死在库里而不是每次调用时重新算一遍。省掉的那点时间是小事关键是排除了初始化代码出错的隐患。还有个小技巧值得分享不确定该用多少点时就跑 n 点和 n1 点的结果两个结果差小于容差就收手这就是自适应高斯求积的雏形也是我调试数值代码时最常用的一个判断手段。