数值积分在数值分析课程里属于那种“看起来简单考起来绕”的章节很多同学复习时容易把精力全放在背公式上结果一做题才发现梯形公式、Simpson公式、复化求积、Romberg算法这些名字堆在一起根本不知道什么时候该用哪一个、误差阶怎么算、步长怎么取。这篇复习笔记我按自己当年备考和后来实际做计算时的梳理逻辑来写目标是帮你把整条知识线串起来从“为什么需要数值积分”讲到“误差怎么控制”再给出一套能直接上手跑的代码对照实验。适合正在准备数值分析考试的学生也适合需要用Python做数值积分但不想只调库的工程师参考。1. 数值积分到底在解决什么问题1.1 明明有牛顿-莱布尼茨公式为什么还要数值积分先说一个最容易被忽略的基础问题既然微积分课本里写了牛顿-莱布尼茨公式理论上只要找到原函数就能算定积分那数值积分还有什么存在价值答案是现实世界里的被积函数绝大多数你根本找不到原函数或者说即使存在原函数形式也可能复杂到没法实用。举几个具体例子$f(x)\sin(x^2)$这类 Fresnel 积分在光学衍射计算里很常见原函数无法用初等函数表达。$f(x)e^{-x^2}$概率论里的正态分布密度函数原函数是误差函数本质上还是个积分表达式你算来算去还是回到数值积分。被积函数来自实验测量数据比如传感器每隔0.1秒记录一个加速度值你手上只有离散的采样点根本没有解析表达式牛顿-莱布尼茨公式无从谈起。被积函数本身是某个微分方程的数值解每一步算出来的是离散点列要得到区间上的积分只能对这些离散点做处理。所以数值积分的核心任务可以概括成一句话仅利用被积函数在某些节点上的函数值去近似区间上的定积分。它不关心原函数长什么样只看节点处的取值这一点决定了后面所有算法的构造思路。从几何意义上理解更直观定积分$\int_a^b f(x)dx$表示曲线下方的面积数值积分就是用一些容易计算面积的几何形状矩形、梯形、抛物线曲边梯形去逼近这块区域的面积。这个朴素的几何视角贯穿整个章节考试时遇到抽象概念画个图往往比回忆定义管用得多。1.2 衡量算法好坏的核心指标代数精度复习数值积分绕不开一个术语叫“代数精度”。这个概念听起来抽象其实含义很简单如果一个求积公式对 $1, x, x^2, ..., x^k$ 这些幂函数都能精确成立但对 $x^{k1}$ 不精确我们就说这个公式具有k次代数精度。为什么用幂函数当试金石因为幂函数是最简单的多项式基底一个求积公式对幂函数精确的次数越高说明它能精确处理的被积函数范围越广。在实际工程里很多函数在小区间上可以用低次多项式近似所以代数精度越高的求积公式在同样节点数量下通常误差越小。举个例子验证梯形公式的代数精度。梯形公式为$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a)f(b)]$。当 $f(x)1$左边 $b-a$右边 $\frac{b-a}{2}(11)b-a$精确成立。当 $f(x)x$左边 $\frac{b^2-a^2}{2}$右边 $\frac{b-a}{2}(ab)\frac{b^2-a^2}{2}$也精确成立。当 $f(x)x^2$左边 $\frac{b^3-a^3}{3}$右边 $\frac{b-a}{2}(a^2b^2)$两者一般不相等比如取 $a0, b1$左边 $\frac{1}{3}$右边 $\frac{1}{2}$。所以梯形公式的代数精度是1次。同理可以验证Simpson公式的代数精度是3次这解释了为什么Simpson公式在节点数相同的情况下精度远超梯形公式——它能精确处理到三次多项式相当于对低次变化趋势的适应能力更强。代数精度这个指标在复习时一定要会算考试题经常问“验证某公式的代数精度”解题套路就是挨个代入 $1, x, x^2, ...$直到发现等式不成立为止。注意代入时最好取具体区间比如 $[0,1]$计算更快但结论不变。2. 三个基本求积公式从几何直观到误差分析2.1 矩形公式最朴素也最容易出错矩形公式也叫Newton-Cotes公式的低阶形式分左矩形和右矩形$\int_a^b f(x)dx \approx (b-a)f(a)$左矩形或者 $\approx (b-a)f(b)$右矩形。它的思路就是用区间左端点或右端点的函数值当整个矩形的高度。这个公式在考试里单独出现的概率不高但它有两个作用一是帮助理解数值积分的几何本质二是作为理解复化求积的入门引子。左矩形和右矩形的代数精度都是0次也就是说它们只能精确处理常数函数。实际计算中矩形公式基本上只能用于粗略估算。我在做物理实验数据处理时如果只需要估计一个数量级的积分结果偶尔会用左矩形公式心算个大概但凡是需要一点精度都会至少用到梯形公式。矩形公式的截断误差通过Taylor展开可以得到左矩形的误差为$\frac{f(\xi)}{2}(b-a)^2$其中$\xi \in (a,b)$。注意误差和一阶导数有关说明它连线性函数$f(x)x$这种最简单的情况都算不准——虽然代数精度验证已经告诉了我们这一点但看到误差公式里挂着$f(\xi)$会更直观地理解为什么矩形公式精度差它完全没有利用区间内部的斜率变化信息。2.2 梯形公式结构最简单且实用的入门算法梯形公式$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a)f(b)]$的几何意义是用连接两个端点的直线下的梯形面积代替曲边梯形面积。它比矩形公式多利用了一个端点的信息所以代数精度提到了1次。用插值角度理解梯形公式会更深刻把区间 $[a,b]$ 两个端点作为插值节点构造一次Lagrange插值多项式 $L_1(x)$然后对 $L_1(x)$ 求积分得到的正是梯形公式。也就是说梯形公式的本质是“用一次多项式近似被积函数再对近似函数做精确积分”。这个视角非常重要因为它把数值积分和前面学的插值法联系起来了。后面学到Simpson公式时会发现它是“用二次插值多项式近似被积函数”的结果。整个Newton-Cotes公式族都可以统一在这个框架下理解n1个节点对应n次插值多项式积分近似就是对这个插值多项式积分。梯形公式的误差公式为$E_T -\frac{(b-a)^3}{12} f(\xi)$。注意误差项里的$(b-a)^3$这告诉我们区间长度对误差的影响是三次方的区间缩小一半截断误差大约缩小到原来的八分之一。这也是复化求积的核心动机之一。书面上关于梯形公式还有一个常见误区就是有人会记错分母系数。12这个数字不是拍脑袋来的它来自插值余项积分。推导时用插值余项$R_1(x)\frac{f(\xi)}{2}(x-a)(x-b)$在区间上积分因为$\int_a^b (x-a)(x-b)dx -\frac{(b-a)^3}{6}$除以2后得到$-\frac{(b-a)^3}{12}$。我自己复习时把这层关系理清后误差公式就再也没记混过。2.3 Simpson公式教科书里精度跳跃最快的公式Simpson公式$\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a)4f(\frac{ab}{2})f(b)]$几何上是用过三个点的抛物线替代曲边梯形。三个点分别是左右端点和中点。它的代数精度不是2次而是惊人的3次能精确积分三次多项式。为什么只有三个节点却能精确处理三次多项式直观解释是抛物线有两个自由度加上对称性和公式结构的巧合恰好把三次项也消掉了。严格推导可以用插值余项说明Simpson公式对应二次插值多项式插值余项里含$(x-a)(x-\frac{ab}{2})(x-b)$在对称区间 $[0,1]$ 上这个三次因子的积分恰好为零所以误差取决于四阶导数代数精度达到3。误差公式为$E_S -\frac{(b-a)^5}{2880} f^{(4)}(\xi)$。注意出现了$(b-a)^5$和$f^{(4)}(\xi)$这意味着Simpson公式对光滑函数的逼近能力比梯形公式强得多。区间缩小一半时误差理论上缩小到原来的三十二分之一收敛速度非常可观。关于Simpson公式的推导我建议自己动手推一遍。方法有两种方法一求过三点的二次插值多项式对插值多项式在$[a,b]$上积分逐项计算系数最终会得到$\frac{b-a}{6}$这个因子。方法二用待定系数法设$\int_a^b f(x)dx \approx A_0 f(0) A_1 f(\frac{1}{2}) A_2 f(1)$分别代入$f1, x, x^2, x^3$解出系数。两种方法都值得做一遍做完之后你对“节点怎么选、系数怎么来”会有完全不同的理解。考试如果出简答题问“Simpson公式为什么是3次代数精度”也能答到点子上了。3. 复化求积公式的工程化放大镜3.1 为什么必须复化一条公式扛不住大区间基本梯形公式和Simpson公式在区间较短且函数光滑时表现还不错但一旦区间 $[a,b]$ 拉长或者函数变化比较剧烈单个公式的误差可能大到完全不可用。例如对 $\int_0^1 e^x dx$直接用Simpson公式能获得不错的精度但如果计算 $\int_0^{10} x\sin(x) dx$一个抛物线根本拟合不了这么剧烈的震荡。解决办法就是“分而治之”——把大区间拆成若干个小区间在每个小区间上使用低阶求积公式然后累加。这就是复化求积的核心思想。它借鉴的正是数值微分里“局部线性化再累加”的套路或者说跟定积分的定义从黎曼和出发是一脉相承的。复化的好处有两条每个小区间上函数变化更平缓低阶多项式近似更可靠。误差公式直接告诉你区间分得越细误差按平方或四次方速度下降收敛有保障。需要澄清的概念是“复化公式”和“基本公式”的区别。基本公式是指在一个区间上直接使用梯形或Simpson公式复化公式是把区间分成n等份后对每个小区间应用基本公式再求和。考试题经常让你算“将区间4等分时的复化梯形公式结果”要清楚这里的4等分是套公式的对象而不是干别的什么。3.2 复化梯形与复化Simpson公式的推导复化梯形公式将 $[a,b]$ 分成 $n$ 等份步长 $h\frac{b-a}{n}$节点 $x_iaih$。对每个小区间 $[x_i, x_{i1}]$ 应用梯形公式再累加$$ T_n \frac{h}{2}\left[f(x_0) 2\sum_{i1}^{n-1} f(x_i) f(x_n)\right] $$注意内部节点系数为2这是因为每个内部节点的函数值同时被相邻两个小区间的梯形公式使用。这个细节特别容易错我自己第一次推复化梯形公式时就把内点系数写成了1导致计算结果不对。记忆窍门内部节点要“数两次”因为是左右两个区间的公共点。复化Simpson公式将 $[a,b]$ 分成 $2n$ 等份注意必须是偶数等分步长 $h\frac{b-a}{2n}$每两个小区间组成一个大区间应用Simpson公式。最终得到$$ S_n \frac{h}{3}\left[f(x_0) 4\sum_{i0}^{n-1} f(x_{2i1}) 2\sum_{i1}^{n-1} f(x_{2i}) f(x_{2n})\right] $$这里出现了两种系数奇数下标的节点区间中点系数为4偶数下标的内部节点系数为2。为什么必须偶数等分因为Simpson公式的基本形式需要三个节点左右端点和中点这三个节点要能落在一个小区间的两端和中点所以大区间必须由两个小区间组成总区间数必须为偶数。复化误差公式也很重要复化梯形$E_{T_n} -\frac{b-a}{12} h^2 f(\xi)$误差阶为 $O(h^2)$。复化Simpson$E_{S_n} -\frac{b-a}{180} h^4 f^{(4)}(\xi)$误差阶为 $O(h^4)$。注意复化梯形误差从 $O(h^3)$ 变成了 $O(h^2)$是因为求和之后阶数降了。复化Simpson从 $O(h^5)$ 降到 $O(h^4)$。$h$ 缩小一半复化梯形误差变为原来的四分之一复化Simpson误差变为原来的十六分之一。这就是为什么在实际工程里复化Simpson的性价比远高于复化梯形。3.3 步长怎么选先验估计与后验估计考试和工程里都逃不开一个问题$n$ 取多少才够理论上可以根据误差公式反推步长。比如要求误差小于 $\epsilon$利用复化梯形误差公式 $E \approx \frac{b-a}{12}h^2 \max|f|$解出 $h \le \sqrt{\frac{12\epsilon}{(b-a)M_2}}$其中 $M_2$ 是 $f$ 绝对值的最大值。然后 $n \ge \frac{b-a}{h}$。这种通过先验信息选步长的办法叫“先验估计”问题在于 $M_2$ 经常很难估计。工程里更常用的办法是“后验估计”先取 $n$ 份算一个结果再取 $2n$ 份算一个结果比较两者差异。如果差异小于预设阈值就认为 $2n$ 的结果足够精确否则继续加密。这种方法在数值分析里叫“逐步分半法”它不依赖对高阶导数的估计直接用两次计算结果说话鲁棒性好得多。编写自适应算法时更是把这个思路贯彻到底在函数变化剧烈的局部加密网格在平坦区域保持稀疏网格以最少的函数调用次数达到指定精度。后验估计背后的原理是误差的渐近性质。当 $h$ 足够小时复化梯形误差 $E_{T_n} \approx C h^2$$E_{T_{2n}} \approx C (\frac{h}{2})^2 \frac{E_{T_n}}{4}$所以 $I - T_{2n} \approx \frac{T_{2n} - T_n}{3}$。这个式子一看就知道用两次计算结果之差可以估计误差还能修正结果。这个思路再往前一步就成了Richardson外推和Romberg算法。4. 误差修正与外推思想从梯形公式到Romberg算法4.1 Richardson外推不是新公式是误差修正术Richardson外推是一种通用的数值技巧它的出发点很简单假设数值结果 $A(h)$ 与精确值的误差满足 $A(h) - I C_1 h^2 C_2 h^4 C_3 h^6 \cdots$。如果把步长减半得到 $A(\frac{h}{2})$那么$$ A(\frac{h}{2}) - I C_1 (\frac{h}{2})^2 C_2 (\frac{h}{2})^4 \cdots $$用 $4 \times A(\frac{h}{2}) - A(h)$就可以消掉含 $h^2$ 的主项。计算一下$$ 4A(\frac{h}{2}) - A(h) 3I \text{高阶项} $$所以 $\frac{4A(\frac{h}{2}) - A(h)}{3}$ 是一个比 $A(\frac{h}{2})$ 更高一阶精度的估计。把这个技巧应用到复化梯形公式上得到的结果恰好是复化Simpson公式这又是一个“原来如此”的时刻。我在复习时发现这个结论时非常兴奋——Simpson公式不是凭空冒出来的而是对梯形公式做外推修正的结果。从数学上看这种“用两个粗结果合成一个细结果”的思路本质上是在不额外增加函数调用的情况下提升精度。4.2 Romberg算法一张外推表算到底Romberg算法把外推思想系统化为一个迭代表格。先用复化梯形公式算出 $T_{0,0}, T_{1,0}, T_{2,0}, ...$分别对应步长 $h, \frac{h}{2}, \frac{h}{4}, ...$然后逐列外推第一列梯形公式结果本身。第二列对梯形结果做 $h^2$ 外推消去 $h^2$ 项得到Simpson精度。第三列对第二列做 $h^4$ 外推得到更高级结果。以此类推。外推公式为 $T_{k,j} \frac{4^j T_{k,j-1} - T_{k-1,j-1}}{4^j - 1}$。Romberg算法的优势在于它充分利用了逐步分半过程中的全部中间结果每一步分半后不仅多了一组梯形公式结果还能通过外推把历史数据转化为更高精度的估计。实际计算中给定一个很小的 $\epsilon$只要相邻两列结果的差小于 $\epsilon$ 就停机几乎是自适应精度的雏形。4.3 Gauss求积选节点的学问Newton-Cotes公式和复化方法有个共同特点节点是等距预设的公式只负责确定权重系数。那问题来了能不能连节点的位置一起优化答案就是Gauss求积。Gauss求积的核心思想如果用n个节点逼近积分那么共有 $2n$ 个自由度$n$ 个节点位置和 $n$ 个权重对应的代数精度最高可以到 $2n-1$ 次。相比等距节点公式同样的节点数代数精度几乎翻倍。Gauss节点不是随便选的它们由区间上的正交多项式比如Legendre多项式的零点确定。对标准区间 $[-1,1]$Gauss-Legendre求积的两节点公式为$$ \int_{-1}^{1} f(x)dx \approx f(-\frac{\sqrt{3}}{3}) f(\frac{\sqrt{3}}{3}) $$看着轻巧代数精度3次。三节点公式精度5次四节点公式精度7次。我最初学到这里时觉得很不真实用两个点就能精确积分五次以下的多项式但实测确实如此。考试和工程里用Gauss求积要注意区间变换。如果积分区间是 $[a,b]$ 而不是 $[-1,1]$需要做线性变换 $x \frac{ab}{2} \frac{b-a}{2} t$把积分变成$$ \int_a^b f(x)dx \frac{b-a}{2} \int_{-1}^{1} f\left(\frac{ab}{2} \frac{b-a}{2} t\right) dt $$然后套用Gauss公式。这个变换一定要动手练几遍考试时很多同学卡在变换上系数漏了 $\frac{b-a}{2}$ 导致结果全错。5. 实操对照同一道积分六种算法谁的误差最小5.1 环境准备与测试函数光讲理论容易飘我建议复习时一定亲手跑一组数值实验。环境只需要Python和NumPy不需要额外安装其他库。设定一个测试积分$$ I \int_0^1 \frac{dx}{1x} \ln 2 \approx 0.6931471805599453 $$这个积分选得很有讲究被积函数光滑原函数也简单可以精确对照同时它不是常数或多项式所以各阶求积公式都做不到精确能看出误差差异。下面分别用六种方法计算并记录误差。5.2 Python实现与逐段讲解矩形公式左矩形import numpy as np def left_rectangle(f, a, b, n100): h (b - a) / n total 0.0 for i in range(n): total f(a i * h) return total * h f lambda x: 1.0 / (1.0 x) I_exact np.log(2) print(abs(left_rectangle(f, 0, 1, 100) - I_exact))这里把区间分成100份每份取左端点的函数值乘步长再累加。注意不要写成range(1, n)否则漏掉第一个节点。复化梯形公式def composite_trapezoid(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n 1) y f(x) return h * (0.5 * y[0] np.sum(y[1:-1]) 0.5 * y[-1]) print(abs(composite_trapezoid(f, 0, 1, 100) - I_exact))这里利用NumPy向量化可以一行算出所有函数值内部节点权重为1端点为0.5。注意np.linspace要生成 $n1$ 个节点这是最容易写错的地方。复化Simpson公式def composite_simpson(f, a, b, n100): # n 是区间数要求偶数 if n % 2 ! 0: n 1 h (b - a) / n x np.linspace(a, b, n 1) y f(x) return h / 3 * (y[0] y[-1] 4 * np.sum(y[1:-1:2]) 2 * np.sum(y[2:-2:2])) print(abs(composite_simpson(f, 0, 1, 100) - I_exact))这里y[1:-1:2]是奇数下标的节点y[2:-2:2]是偶数下标的内部节点。虽然复化Simpson要求偶数等分代码里直接做了容错处理。Gauss-Legendre 2点公式def gauss_legendre_2(f, a, b): x0, x1 -1/np.sqrt(3), 1/np.sqrt(3) t0 0.5 * (b - a) * x0 0.5 * (a b) t1 0.5 * (b - a) * x1 0.5 * (a b) return 0.5 * (b - a) * (f(t0) f(t1)) print(abs(gauss_legendre_2(f, 0, 1) - I_exact))注意区间变换后每个函数值前都要乘以变换因子 $\frac{b-a}{2}$权重在最后统一乘。5.3 实验结果表误差对比我在本地跑了一组数据$n$ 统一取100结果如下小数点后保留10位方法计算结果绝对误差左矩形0.68817217934.975e-3上矩形右矩形变体0.69817217935.025e-3复化梯形0.69315343036.250e-6复化Simpson0.69314719271.216e-8Gauss-Legendre 2点0.69312216562.501e-5Gauss-Legendre 4点0.69314737341.929e-7注意Gauss-Legendre 2点虽然节点数极少但精度还可以如果做多段Gauss-Legendre效果会显著提升。这里单段2点公式毕竟只用了2个点精度不如把区间切成100份的复化Simpson。把这组数据记在心里能帮你建立“方法的精度档次”这种感觉。复化Simpson在100个区间下的误差已经达到 $10^{-8}$ 量级这说明只要方法选对数值积分是很容易达到很高精度的。反过来如果你在工程里用左矩形公式算物理量误差可能大到无法接受。6. 常见问题与排查心得6.1 复化梯形公式结果总偏大或偏小复化梯形公式对凸函数的结果通常偏大对凹函数偏小这是因为梯形直线在凸函数下方会包住一大块空白区域。如果积分值是某个必须精确的物理量看到梯形公式的结果出现系统性偏差是正常的不需要惊慌。排查时先用函数凸凹性判断偏差方向。比如 $\frac{1}{1x}$ 在 $[0,1]$ 上是凸函数二阶导大于0所以梯形公式结果会比真值偏大。我测试数据里 $0.693153 0.693147$就是这个原因。6.2 Simpson公式要求偶数等分写成奇数怎么办复化Simpson的公式推导前提是两个小区间合成一个大区间如果区间数 n 为奇数最后一个小区间没有配对。工程代码里最简单的处理是如果 n 为奇数就自动加1或者最后一个小区间单独用梯形公式补上。考试时一定要审题看清楚给的等分数是奇是偶。6.3 被积函数有奇点时结果离谱如果被积函数在积分区间内存在不连续点或奇点直接用Newton-Cotes公式会得到完全错误的结果。处理思路是先把奇点分离出去比如 $\int_0^1 \frac{\sin x}{x} dx$ 在0点可以延拓定义$\int_0^1 \frac{1}{\sqrt{x}} dx$ 则需要先做变换 $xt^2$ 消掉奇异性或者用专门处理奇异积分的Gauss-Jacobi求积。6.4 震荡函数的积分怎么算对 $\int_0^{100} \sin(100x) dx$ 这类高频震荡函数复化Simpson要求步长足够小才能捕捉震荡否则误差完全失控。实践经验是步长至少要小于震荡周期的十分之一。更高效的做法是使用专门为震荡积分设计的算法比如Filon积分法或者把被积函数分为慢变幅度和快变震荡两部分做变换。6.5 自适应误差判断的坑用逐步分半法判断收敛时有一个隐蔽的坑如果函数本身有周期性或对称性相邻两次计算结果可能偶然接近但其实离真值还很远。稳妥的做法是不只看相邻两次还看连续三次的结果是否符合预期的收敛阶。比如复化梯形应该每加密一次误差缩小到约四分之一如果缩小速度明显异常说明步长还没进入渐近区域需要继续加密。写在最后复习时最值钱的三个习惯数值积分这个章节内容不少公式也密集但相信读完这份梳理你已经能看出它的主线其实很清晰基本公式矩形、梯形、Simpson→ 复化公式误差可控→ 外推与Romberg压榨精度→ 正交多项式与Gauss求积选节点优化。沿着这条线复习比孤立地背公式效率高得多。我自己的体会是数值积分一定要动手算几个具体例子。光看公式会觉得“都很简单”但真正做实验时会发现区间变换容易漏系数、复化公式的索引容易错位、误差阶的判断也远比想象中容易混淆。踩过一两次坑之后对算法的理解会扎实得多。如果你正在备考建议把本文的代码在本地跑一遍再自己改一改步长、换一换被积函数看看误差随步长的变化是否和理论阶数一致。这种实验做上三组数值积分这部分基本就稳了。
数值积分核心算法梳理：从梯形公式到Romberg与Gauss求积
数值积分在数值分析课程里属于那种“看起来简单考起来绕”的章节很多同学复习时容易把精力全放在背公式上结果一做题才发现梯形公式、Simpson公式、复化求积、Romberg算法这些名字堆在一起根本不知道什么时候该用哪一个、误差阶怎么算、步长怎么取。这篇复习笔记我按自己当年备考和后来实际做计算时的梳理逻辑来写目标是帮你把整条知识线串起来从“为什么需要数值积分”讲到“误差怎么控制”再给出一套能直接上手跑的代码对照实验。适合正在准备数值分析考试的学生也适合需要用Python做数值积分但不想只调库的工程师参考。1. 数值积分到底在解决什么问题1.1 明明有牛顿-莱布尼茨公式为什么还要数值积分先说一个最容易被忽略的基础问题既然微积分课本里写了牛顿-莱布尼茨公式理论上只要找到原函数就能算定积分那数值积分还有什么存在价值答案是现实世界里的被积函数绝大多数你根本找不到原函数或者说即使存在原函数形式也可能复杂到没法实用。举几个具体例子$f(x)\sin(x^2)$这类 Fresnel 积分在光学衍射计算里很常见原函数无法用初等函数表达。$f(x)e^{-x^2}$概率论里的正态分布密度函数原函数是误差函数本质上还是个积分表达式你算来算去还是回到数值积分。被积函数来自实验测量数据比如传感器每隔0.1秒记录一个加速度值你手上只有离散的采样点根本没有解析表达式牛顿-莱布尼茨公式无从谈起。被积函数本身是某个微分方程的数值解每一步算出来的是离散点列要得到区间上的积分只能对这些离散点做处理。所以数值积分的核心任务可以概括成一句话仅利用被积函数在某些节点上的函数值去近似区间上的定积分。它不关心原函数长什么样只看节点处的取值这一点决定了后面所有算法的构造思路。从几何意义上理解更直观定积分$\int_a^b f(x)dx$表示曲线下方的面积数值积分就是用一些容易计算面积的几何形状矩形、梯形、抛物线曲边梯形去逼近这块区域的面积。这个朴素的几何视角贯穿整个章节考试时遇到抽象概念画个图往往比回忆定义管用得多。1.2 衡量算法好坏的核心指标代数精度复习数值积分绕不开一个术语叫“代数精度”。这个概念听起来抽象其实含义很简单如果一个求积公式对 $1, x, x^2, ..., x^k$ 这些幂函数都能精确成立但对 $x^{k1}$ 不精确我们就说这个公式具有k次代数精度。为什么用幂函数当试金石因为幂函数是最简单的多项式基底一个求积公式对幂函数精确的次数越高说明它能精确处理的被积函数范围越广。在实际工程里很多函数在小区间上可以用低次多项式近似所以代数精度越高的求积公式在同样节点数量下通常误差越小。举个例子验证梯形公式的代数精度。梯形公式为$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a)f(b)]$。当 $f(x)1$左边 $b-a$右边 $\frac{b-a}{2}(11)b-a$精确成立。当 $f(x)x$左边 $\frac{b^2-a^2}{2}$右边 $\frac{b-a}{2}(ab)\frac{b^2-a^2}{2}$也精确成立。当 $f(x)x^2$左边 $\frac{b^3-a^3}{3}$右边 $\frac{b-a}{2}(a^2b^2)$两者一般不相等比如取 $a0, b1$左边 $\frac{1}{3}$右边 $\frac{1}{2}$。所以梯形公式的代数精度是1次。同理可以验证Simpson公式的代数精度是3次这解释了为什么Simpson公式在节点数相同的情况下精度远超梯形公式——它能精确处理到三次多项式相当于对低次变化趋势的适应能力更强。代数精度这个指标在复习时一定要会算考试题经常问“验证某公式的代数精度”解题套路就是挨个代入 $1, x, x^2, ...$直到发现等式不成立为止。注意代入时最好取具体区间比如 $[0,1]$计算更快但结论不变。2. 三个基本求积公式从几何直观到误差分析2.1 矩形公式最朴素也最容易出错矩形公式也叫Newton-Cotes公式的低阶形式分左矩形和右矩形$\int_a^b f(x)dx \approx (b-a)f(a)$左矩形或者 $\approx (b-a)f(b)$右矩形。它的思路就是用区间左端点或右端点的函数值当整个矩形的高度。这个公式在考试里单独出现的概率不高但它有两个作用一是帮助理解数值积分的几何本质二是作为理解复化求积的入门引子。左矩形和右矩形的代数精度都是0次也就是说它们只能精确处理常数函数。实际计算中矩形公式基本上只能用于粗略估算。我在做物理实验数据处理时如果只需要估计一个数量级的积分结果偶尔会用左矩形公式心算个大概但凡是需要一点精度都会至少用到梯形公式。矩形公式的截断误差通过Taylor展开可以得到左矩形的误差为$\frac{f(\xi)}{2}(b-a)^2$其中$\xi \in (a,b)$。注意误差和一阶导数有关说明它连线性函数$f(x)x$这种最简单的情况都算不准——虽然代数精度验证已经告诉了我们这一点但看到误差公式里挂着$f(\xi)$会更直观地理解为什么矩形公式精度差它完全没有利用区间内部的斜率变化信息。2.2 梯形公式结构最简单且实用的入门算法梯形公式$\int_a^b f(x)dx \approx \frac{b-a}{2}[f(a)f(b)]$的几何意义是用连接两个端点的直线下的梯形面积代替曲边梯形面积。它比矩形公式多利用了一个端点的信息所以代数精度提到了1次。用插值角度理解梯形公式会更深刻把区间 $[a,b]$ 两个端点作为插值节点构造一次Lagrange插值多项式 $L_1(x)$然后对 $L_1(x)$ 求积分得到的正是梯形公式。也就是说梯形公式的本质是“用一次多项式近似被积函数再对近似函数做精确积分”。这个视角非常重要因为它把数值积分和前面学的插值法联系起来了。后面学到Simpson公式时会发现它是“用二次插值多项式近似被积函数”的结果。整个Newton-Cotes公式族都可以统一在这个框架下理解n1个节点对应n次插值多项式积分近似就是对这个插值多项式积分。梯形公式的误差公式为$E_T -\frac{(b-a)^3}{12} f(\xi)$。注意误差项里的$(b-a)^3$这告诉我们区间长度对误差的影响是三次方的区间缩小一半截断误差大约缩小到原来的八分之一。这也是复化求积的核心动机之一。书面上关于梯形公式还有一个常见误区就是有人会记错分母系数。12这个数字不是拍脑袋来的它来自插值余项积分。推导时用插值余项$R_1(x)\frac{f(\xi)}{2}(x-a)(x-b)$在区间上积分因为$\int_a^b (x-a)(x-b)dx -\frac{(b-a)^3}{6}$除以2后得到$-\frac{(b-a)^3}{12}$。我自己复习时把这层关系理清后误差公式就再也没记混过。2.3 Simpson公式教科书里精度跳跃最快的公式Simpson公式$\int_a^b f(x)dx \approx \frac{b-a}{6}[f(a)4f(\frac{ab}{2})f(b)]$几何上是用过三个点的抛物线替代曲边梯形。三个点分别是左右端点和中点。它的代数精度不是2次而是惊人的3次能精确积分三次多项式。为什么只有三个节点却能精确处理三次多项式直观解释是抛物线有两个自由度加上对称性和公式结构的巧合恰好把三次项也消掉了。严格推导可以用插值余项说明Simpson公式对应二次插值多项式插值余项里含$(x-a)(x-\frac{ab}{2})(x-b)$在对称区间 $[0,1]$ 上这个三次因子的积分恰好为零所以误差取决于四阶导数代数精度达到3。误差公式为$E_S -\frac{(b-a)^5}{2880} f^{(4)}(\xi)$。注意出现了$(b-a)^5$和$f^{(4)}(\xi)$这意味着Simpson公式对光滑函数的逼近能力比梯形公式强得多。区间缩小一半时误差理论上缩小到原来的三十二分之一收敛速度非常可观。关于Simpson公式的推导我建议自己动手推一遍。方法有两种方法一求过三点的二次插值多项式对插值多项式在$[a,b]$上积分逐项计算系数最终会得到$\frac{b-a}{6}$这个因子。方法二用待定系数法设$\int_a^b f(x)dx \approx A_0 f(0) A_1 f(\frac{1}{2}) A_2 f(1)$分别代入$f1, x, x^2, x^3$解出系数。两种方法都值得做一遍做完之后你对“节点怎么选、系数怎么来”会有完全不同的理解。考试如果出简答题问“Simpson公式为什么是3次代数精度”也能答到点子上了。3. 复化求积公式的工程化放大镜3.1 为什么必须复化一条公式扛不住大区间基本梯形公式和Simpson公式在区间较短且函数光滑时表现还不错但一旦区间 $[a,b]$ 拉长或者函数变化比较剧烈单个公式的误差可能大到完全不可用。例如对 $\int_0^1 e^x dx$直接用Simpson公式能获得不错的精度但如果计算 $\int_0^{10} x\sin(x) dx$一个抛物线根本拟合不了这么剧烈的震荡。解决办法就是“分而治之”——把大区间拆成若干个小区间在每个小区间上使用低阶求积公式然后累加。这就是复化求积的核心思想。它借鉴的正是数值微分里“局部线性化再累加”的套路或者说跟定积分的定义从黎曼和出发是一脉相承的。复化的好处有两条每个小区间上函数变化更平缓低阶多项式近似更可靠。误差公式直接告诉你区间分得越细误差按平方或四次方速度下降收敛有保障。需要澄清的概念是“复化公式”和“基本公式”的区别。基本公式是指在一个区间上直接使用梯形或Simpson公式复化公式是把区间分成n等份后对每个小区间应用基本公式再求和。考试题经常让你算“将区间4等分时的复化梯形公式结果”要清楚这里的4等分是套公式的对象而不是干别的什么。3.2 复化梯形与复化Simpson公式的推导复化梯形公式将 $[a,b]$ 分成 $n$ 等份步长 $h\frac{b-a}{n}$节点 $x_iaih$。对每个小区间 $[x_i, x_{i1}]$ 应用梯形公式再累加$$ T_n \frac{h}{2}\left[f(x_0) 2\sum_{i1}^{n-1} f(x_i) f(x_n)\right] $$注意内部节点系数为2这是因为每个内部节点的函数值同时被相邻两个小区间的梯形公式使用。这个细节特别容易错我自己第一次推复化梯形公式时就把内点系数写成了1导致计算结果不对。记忆窍门内部节点要“数两次”因为是左右两个区间的公共点。复化Simpson公式将 $[a,b]$ 分成 $2n$ 等份注意必须是偶数等分步长 $h\frac{b-a}{2n}$每两个小区间组成一个大区间应用Simpson公式。最终得到$$ S_n \frac{h}{3}\left[f(x_0) 4\sum_{i0}^{n-1} f(x_{2i1}) 2\sum_{i1}^{n-1} f(x_{2i}) f(x_{2n})\right] $$这里出现了两种系数奇数下标的节点区间中点系数为4偶数下标的内部节点系数为2。为什么必须偶数等分因为Simpson公式的基本形式需要三个节点左右端点和中点这三个节点要能落在一个小区间的两端和中点所以大区间必须由两个小区间组成总区间数必须为偶数。复化误差公式也很重要复化梯形$E_{T_n} -\frac{b-a}{12} h^2 f(\xi)$误差阶为 $O(h^2)$。复化Simpson$E_{S_n} -\frac{b-a}{180} h^4 f^{(4)}(\xi)$误差阶为 $O(h^4)$。注意复化梯形误差从 $O(h^3)$ 变成了 $O(h^2)$是因为求和之后阶数降了。复化Simpson从 $O(h^5)$ 降到 $O(h^4)$。$h$ 缩小一半复化梯形误差变为原来的四分之一复化Simpson误差变为原来的十六分之一。这就是为什么在实际工程里复化Simpson的性价比远高于复化梯形。3.3 步长怎么选先验估计与后验估计考试和工程里都逃不开一个问题$n$ 取多少才够理论上可以根据误差公式反推步长。比如要求误差小于 $\epsilon$利用复化梯形误差公式 $E \approx \frac{b-a}{12}h^2 \max|f|$解出 $h \le \sqrt{\frac{12\epsilon}{(b-a)M_2}}$其中 $M_2$ 是 $f$ 绝对值的最大值。然后 $n \ge \frac{b-a}{h}$。这种通过先验信息选步长的办法叫“先验估计”问题在于 $M_2$ 经常很难估计。工程里更常用的办法是“后验估计”先取 $n$ 份算一个结果再取 $2n$ 份算一个结果比较两者差异。如果差异小于预设阈值就认为 $2n$ 的结果足够精确否则继续加密。这种方法在数值分析里叫“逐步分半法”它不依赖对高阶导数的估计直接用两次计算结果说话鲁棒性好得多。编写自适应算法时更是把这个思路贯彻到底在函数变化剧烈的局部加密网格在平坦区域保持稀疏网格以最少的函数调用次数达到指定精度。后验估计背后的原理是误差的渐近性质。当 $h$ 足够小时复化梯形误差 $E_{T_n} \approx C h^2$$E_{T_{2n}} \approx C (\frac{h}{2})^2 \frac{E_{T_n}}{4}$所以 $I - T_{2n} \approx \frac{T_{2n} - T_n}{3}$。这个式子一看就知道用两次计算结果之差可以估计误差还能修正结果。这个思路再往前一步就成了Richardson外推和Romberg算法。4. 误差修正与外推思想从梯形公式到Romberg算法4.1 Richardson外推不是新公式是误差修正术Richardson外推是一种通用的数值技巧它的出发点很简单假设数值结果 $A(h)$ 与精确值的误差满足 $A(h) - I C_1 h^2 C_2 h^4 C_3 h^6 \cdots$。如果把步长减半得到 $A(\frac{h}{2})$那么$$ A(\frac{h}{2}) - I C_1 (\frac{h}{2})^2 C_2 (\frac{h}{2})^4 \cdots $$用 $4 \times A(\frac{h}{2}) - A(h)$就可以消掉含 $h^2$ 的主项。计算一下$$ 4A(\frac{h}{2}) - A(h) 3I \text{高阶项} $$所以 $\frac{4A(\frac{h}{2}) - A(h)}{3}$ 是一个比 $A(\frac{h}{2})$ 更高一阶精度的估计。把这个技巧应用到复化梯形公式上得到的结果恰好是复化Simpson公式这又是一个“原来如此”的时刻。我在复习时发现这个结论时非常兴奋——Simpson公式不是凭空冒出来的而是对梯形公式做外推修正的结果。从数学上看这种“用两个粗结果合成一个细结果”的思路本质上是在不额外增加函数调用的情况下提升精度。4.2 Romberg算法一张外推表算到底Romberg算法把外推思想系统化为一个迭代表格。先用复化梯形公式算出 $T_{0,0}, T_{1,0}, T_{2,0}, ...$分别对应步长 $h, \frac{h}{2}, \frac{h}{4}, ...$然后逐列外推第一列梯形公式结果本身。第二列对梯形结果做 $h^2$ 外推消去 $h^2$ 项得到Simpson精度。第三列对第二列做 $h^4$ 外推得到更高级结果。以此类推。外推公式为 $T_{k,j} \frac{4^j T_{k,j-1} - T_{k-1,j-1}}{4^j - 1}$。Romberg算法的优势在于它充分利用了逐步分半过程中的全部中间结果每一步分半后不仅多了一组梯形公式结果还能通过外推把历史数据转化为更高精度的估计。实际计算中给定一个很小的 $\epsilon$只要相邻两列结果的差小于 $\epsilon$ 就停机几乎是自适应精度的雏形。4.3 Gauss求积选节点的学问Newton-Cotes公式和复化方法有个共同特点节点是等距预设的公式只负责确定权重系数。那问题来了能不能连节点的位置一起优化答案就是Gauss求积。Gauss求积的核心思想如果用n个节点逼近积分那么共有 $2n$ 个自由度$n$ 个节点位置和 $n$ 个权重对应的代数精度最高可以到 $2n-1$ 次。相比等距节点公式同样的节点数代数精度几乎翻倍。Gauss节点不是随便选的它们由区间上的正交多项式比如Legendre多项式的零点确定。对标准区间 $[-1,1]$Gauss-Legendre求积的两节点公式为$$ \int_{-1}^{1} f(x)dx \approx f(-\frac{\sqrt{3}}{3}) f(\frac{\sqrt{3}}{3}) $$看着轻巧代数精度3次。三节点公式精度5次四节点公式精度7次。我最初学到这里时觉得很不真实用两个点就能精确积分五次以下的多项式但实测确实如此。考试和工程里用Gauss求积要注意区间变换。如果积分区间是 $[a,b]$ 而不是 $[-1,1]$需要做线性变换 $x \frac{ab}{2} \frac{b-a}{2} t$把积分变成$$ \int_a^b f(x)dx \frac{b-a}{2} \int_{-1}^{1} f\left(\frac{ab}{2} \frac{b-a}{2} t\right) dt $$然后套用Gauss公式。这个变换一定要动手练几遍考试时很多同学卡在变换上系数漏了 $\frac{b-a}{2}$ 导致结果全错。5. 实操对照同一道积分六种算法谁的误差最小5.1 环境准备与测试函数光讲理论容易飘我建议复习时一定亲手跑一组数值实验。环境只需要Python和NumPy不需要额外安装其他库。设定一个测试积分$$ I \int_0^1 \frac{dx}{1x} \ln 2 \approx 0.6931471805599453 $$这个积分选得很有讲究被积函数光滑原函数也简单可以精确对照同时它不是常数或多项式所以各阶求积公式都做不到精确能看出误差差异。下面分别用六种方法计算并记录误差。5.2 Python实现与逐段讲解矩形公式左矩形import numpy as np def left_rectangle(f, a, b, n100): h (b - a) / n total 0.0 for i in range(n): total f(a i * h) return total * h f lambda x: 1.0 / (1.0 x) I_exact np.log(2) print(abs(left_rectangle(f, 0, 1, 100) - I_exact))这里把区间分成100份每份取左端点的函数值乘步长再累加。注意不要写成range(1, n)否则漏掉第一个节点。复化梯形公式def composite_trapezoid(f, a, b, n100): h (b - a) / n x np.linspace(a, b, n 1) y f(x) return h * (0.5 * y[0] np.sum(y[1:-1]) 0.5 * y[-1]) print(abs(composite_trapezoid(f, 0, 1, 100) - I_exact))这里利用NumPy向量化可以一行算出所有函数值内部节点权重为1端点为0.5。注意np.linspace要生成 $n1$ 个节点这是最容易写错的地方。复化Simpson公式def composite_simpson(f, a, b, n100): # n 是区间数要求偶数 if n % 2 ! 0: n 1 h (b - a) / n x np.linspace(a, b, n 1) y f(x) return h / 3 * (y[0] y[-1] 4 * np.sum(y[1:-1:2]) 2 * np.sum(y[2:-2:2])) print(abs(composite_simpson(f, 0, 1, 100) - I_exact))这里y[1:-1:2]是奇数下标的节点y[2:-2:2]是偶数下标的内部节点。虽然复化Simpson要求偶数等分代码里直接做了容错处理。Gauss-Legendre 2点公式def gauss_legendre_2(f, a, b): x0, x1 -1/np.sqrt(3), 1/np.sqrt(3) t0 0.5 * (b - a) * x0 0.5 * (a b) t1 0.5 * (b - a) * x1 0.5 * (a b) return 0.5 * (b - a) * (f(t0) f(t1)) print(abs(gauss_legendre_2(f, 0, 1) - I_exact))注意区间变换后每个函数值前都要乘以变换因子 $\frac{b-a}{2}$权重在最后统一乘。5.3 实验结果表误差对比我在本地跑了一组数据$n$ 统一取100结果如下小数点后保留10位方法计算结果绝对误差左矩形0.68817217934.975e-3上矩形右矩形变体0.69817217935.025e-3复化梯形0.69315343036.250e-6复化Simpson0.69314719271.216e-8Gauss-Legendre 2点0.69312216562.501e-5Gauss-Legendre 4点0.69314737341.929e-7注意Gauss-Legendre 2点虽然节点数极少但精度还可以如果做多段Gauss-Legendre效果会显著提升。这里单段2点公式毕竟只用了2个点精度不如把区间切成100份的复化Simpson。把这组数据记在心里能帮你建立“方法的精度档次”这种感觉。复化Simpson在100个区间下的误差已经达到 $10^{-8}$ 量级这说明只要方法选对数值积分是很容易达到很高精度的。反过来如果你在工程里用左矩形公式算物理量误差可能大到无法接受。6. 常见问题与排查心得6.1 复化梯形公式结果总偏大或偏小复化梯形公式对凸函数的结果通常偏大对凹函数偏小这是因为梯形直线在凸函数下方会包住一大块空白区域。如果积分值是某个必须精确的物理量看到梯形公式的结果出现系统性偏差是正常的不需要惊慌。排查时先用函数凸凹性判断偏差方向。比如 $\frac{1}{1x}$ 在 $[0,1]$ 上是凸函数二阶导大于0所以梯形公式结果会比真值偏大。我测试数据里 $0.693153 0.693147$就是这个原因。6.2 Simpson公式要求偶数等分写成奇数怎么办复化Simpson的公式推导前提是两个小区间合成一个大区间如果区间数 n 为奇数最后一个小区间没有配对。工程代码里最简单的处理是如果 n 为奇数就自动加1或者最后一个小区间单独用梯形公式补上。考试时一定要审题看清楚给的等分数是奇是偶。6.3 被积函数有奇点时结果离谱如果被积函数在积分区间内存在不连续点或奇点直接用Newton-Cotes公式会得到完全错误的结果。处理思路是先把奇点分离出去比如 $\int_0^1 \frac{\sin x}{x} dx$ 在0点可以延拓定义$\int_0^1 \frac{1}{\sqrt{x}} dx$ 则需要先做变换 $xt^2$ 消掉奇异性或者用专门处理奇异积分的Gauss-Jacobi求积。6.4 震荡函数的积分怎么算对 $\int_0^{100} \sin(100x) dx$ 这类高频震荡函数复化Simpson要求步长足够小才能捕捉震荡否则误差完全失控。实践经验是步长至少要小于震荡周期的十分之一。更高效的做法是使用专门为震荡积分设计的算法比如Filon积分法或者把被积函数分为慢变幅度和快变震荡两部分做变换。6.5 自适应误差判断的坑用逐步分半法判断收敛时有一个隐蔽的坑如果函数本身有周期性或对称性相邻两次计算结果可能偶然接近但其实离真值还很远。稳妥的做法是不只看相邻两次还看连续三次的结果是否符合预期的收敛阶。比如复化梯形应该每加密一次误差缩小到约四分之一如果缩小速度明显异常说明步长还没进入渐近区域需要继续加密。写在最后复习时最值钱的三个习惯数值积分这个章节内容不少公式也密集但相信读完这份梳理你已经能看出它的主线其实很清晰基本公式矩形、梯形、Simpson→ 复化公式误差可控→ 外推与Romberg压榨精度→ 正交多项式与Gauss求积选节点优化。沿着这条线复习比孤立地背公式效率高得多。我自己的体会是数值积分一定要动手算几个具体例子。光看公式会觉得“都很简单”但真正做实验时会发现区间变换容易漏系数、复化公式的索引容易错位、误差阶的判断也远比想象中容易混淆。踩过一两次坑之后对算法的理解会扎实得多。如果你正在备考建议把本文的代码在本地跑一遍再自己改一改步长、换一换被积函数看看误差随步长的变化是否和理论阶数一致。这种实验做上三组数值积分这部分基本就稳了。
