引言
定积分 $I[f] = \int_a^b f(x),\mathrm{d}x$ 的数值计算贯穿工程分析的各个环节：有限元刚度矩阵的组装需要对每个单元做积分，概率论中的期望与方差计算是积分，辐射传热的角系数、天线的方向图积分、信号处理中的傅里叶系数也都是积分。然而绝大多数被积函数没有初等原函数（如 $e^{-x^2}$、椭圆积分），数值求积成为唯一出路。本文讨论两条主线：以 Gauss 求积为代表的”固定节点高精度”路线，和以自适应 Simpson 为代表的”误差驱动自动加密”路线，并给出精度与代价的定量比较。
理论基础
插值型求积与代数精度
数值求积的一般形式是用函数在若干节点处的值的加权和逼近积分：
$$I_n[f] = \sum_{i=0}^{n} w_i f(x_i) \approx \int_a^b f(x),\mathrm{d}x$$
若取节点上的 Lagrange 插值多项式代替 $f$ 再积分，即得插值型求积公式，权重 $w_i$ 由基函数的积分唯一确定。衡量精度的核心概念是代数精度：若公式对所有次数不超过 $m$ 的多项式精确成立，而对某个 $m+1$ 次多项式不成立，则称其代数精度为 $m$。代数精度越高，对光滑函数的逼近阶越高。
取等距节点的插值型公式即 Newton–Cotes 公式：$n+1$ 个节点时，代数精度在 $n$ 为奇数时为 $n$，在 $n$ 为偶数时因对称性”白捡”一阶而达到 $n+1$。但当 $n \ge 8$ 时权重出现负值，舍入误差被放大，加之等距高次插值本身的 Runge 振荡，高阶 Newton–Cotes 在工程上并不实用。
Gauss 求积：把节点也变成自由度
Newton–Cotes 浪费了节点位置的 $n+1$ 个自由度。若节点与权重共 $2n+2$ 个参数全部用于提高代数精度，理论上可达到 $2n+1$ 阶——这正是 Gauss 求积，且该上界不可超越。Gauss 节点的构造与正交多项式紧密相连：在区间 $[-1,1]$ 上，Gauss–Legendre 公式的节点恰为 Legendre 多项式 $P_{n+1}$ 的 $n+1$ 个零点，权重
$$w_i = \frac{2}{(1 - x_i^2),[P'_{n+1}(x_i)]^2} > 0$$
权重恒正保证了公式的数值稳定性。其误差项为
$$E_n[f] = \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_{-1}^{1}\prod_{i=0}^{n}(x - x_i)^2,\mathrm{d}x$$
即对 $f \in C^{2n+2}$ 的光滑被积函数，误差随 $n$ 以阶乘速度衰减——实践中十几个节点常可达到双精度机器精度，这是 Gauss 求积在有限元程序中被奉为默认方案的原因。任意有限区间 $[a,b]$ 通过仿射变换 $x = \frac{b-a}{2}t + \frac{a+b}{2}$ 化归标准区间，雅可比因子 $\frac{b-a}{2}$ 乘入权重即可。
自适应 Simpson：让误差估计驱动网格
Gauss 公式虽高效，但对含尖点、边界层或端点奇异性的函数，固定节点的全局公式会把大量节点浪费在平缓区域。自适应策略的基本元件是 Simpson 公式 $S(a,b) = \frac{h}{6}[f(a) + 4f(c) + f(b)]$（$c$ 为中点，$h = b-a$）。设把区间二分后两子区间的 Simpson 和为 $S_2$，由 Richardson 外推可得事后误差估计
$$|I - S_2| \approx \frac{1}{15},|S_2 - S_1|$$
其中 $S_1 = S(a,b)$。于是算法逻辑极为朴素：若当前区间的误差估计超过分摊的容差则递归二分，否则接受并用 $\frac{S_2 - S_1}{15}$ 做外推修正。误差被自动集中到”难积”的区域，对 $\sqrt{x}$ 这类端点导数奇异的函数尤其有效。
数值方法与实现
import numpy as np
from numpy.polynomial.legendre import leggauss
def gauss_quad(f, a, b, n):
"""n 点 Gauss-Legendre 求积（含区间仿射变换）"""
x, w = leggauss(n)
t = 0.5*(b - a)*x + 0.5*(b + a)
return 0.5*(b - a)*float(np.dot(w, f(t)))
def _simpson(f, a, b, fa, fb):
c = 0.5*(a + b)
return (b - a)/6.0*(fa + 4.0*f(c) + fb), c
def adaptive_simpson(f, a, b, tol, depth=25):
"""递归自适应 Simpson 求积"""
fa, fb = f(a), f(b)
S1, c = _simpson(f, a, b, fa, fb)
fc = f(c)
S2l, _ = _simpson(f, a, c, fa, fc)
S2r, _ = _simpson(f, c, b, fc, fb)
S2 = S2l + S2r
if depth <= 0 or abs(S2 - S1) <= 15.0*tol:
return S2 + (S2 - S1)/15.0 # Richardson 外推修正
return (adaptive_simpson(f, a, c, tol/2, depth - 1)
+ adaptive_simpson(f, c, b, tol/2, depth - 1))
f = np.exp
I_exact = np.e - 1.0/np.e # ∫_{-1}^{1} e^x dx
for n in (2, 4, 8, 16):
err = abs(gauss_quad(f, -1.0, 1.0, n) - I_exact)
print(f"n={n:2d} Gauss 误差 {err:.3e}")
g = np.sqrt # 端点导数奇异
print(gauss_quad(g, 0.0, 1.0, 64),
adaptive_simpson(g, 0.0, 1.0, 1e-10), 2.0/3.0)运行可以观察到：对整函数 $e^x$，Gauss 误差从 $10^{-3}$ 量级迅速压到机器精度；而对 $\sqrt{x}$，固定 64 点 Gauss 仍有可见误差，自适应 Simpson 则能贴着 $x=0$ 自动加密达到 $10^{-10}$ 容差。
工程案例分析
有限元程序中的单元积分是 Gauss 求积最大的消费场景。以四节点双线性四边形单元为例，刚度矩阵被积函数在参考单元上是 $2\times 2$ 次多项式（不计雅可比变化），$2\times 2$ 点 Gauss 积分（代数精度 3）恰好精确积分；但故意降阶到 $1$ 点积分会产生著名的”沙漏”零能模式，需要额外稳定化——这提醒我们积分阶数不是越高越贵、越低越省的简单权衡，而与单元的谱性质深度耦合。四面体、三角形单元对应的是 Gauss–Hammer 规则，奇异的近场积分（边界元中的 $1/r$ 核）则需配 Duffy 变换或奇异相消技术。
| 求积公式 | 节点数 | 代数精度 | 权重符号 | 高阶行为 |
|---|---|---|---|---|
| 梯形 | 2 | 1 | 正 | 稳定 |
| Simpson | 3 | 3 | 正 | 稳定 |
| Newton–Cotes（$n \ge 8$） | $\ge 9$ | $n$ 或 $n+1$ | 出现负值 | 不稳定，Runge 振荡 |
| Gauss–Legendre | $n+1$ | $2n+1$ | 恒正 | 谱收敛（光滑函数） |
| Clenshaw–Curtis | $n+1$ | $n$ | 恒正 | 节点嵌套，实际精度接近 Gauss |
振荡积分（如边界积分方程中的 Hankel 核）与无穷区间积分则超出经典多项式框架：前者宜用 Filon/Levin 型方法，把 $e^{i\omega x}$ 吸进权函数；后者用 Gauss–Laguerre、Gauss–Hermite 或双指数变换（tanh-sinh 求积），后者对端点奇异性几乎”免疫”。
结论与展望
Gauss 求积以正交多项式理论换取了插值型公式的最高代数精度，是光滑被积函数与规则区域的首选；自适应 Simpson 以事后误差估计换取了对病态行为的鲁棒性，二者互补而非竞争。工程实践中常见的组合策略是：先用自适应探索定位困难子区域，再在子区域内施加高阶 Gauss 规则。向高维推广时，张量积规则遭遇维数灾难，Smolyak 稀疏网格、拟 Monte Carlo 与深度学习辅助的重要性采样是当前活跃的研究方向。
发表回复