企业官网建设流程全解析

数值积分是很多工程计算绕不开的坎。有些函数看似简单，原函数却根本写不出来，比如exp(-x^2)、带振荡项的 Bessel 函数，或者只有离散采样数据而没有解析表达式的情况。这时候就得走近似路线，而牛顿-柯特斯求积公式，就是数值积分里最经典、最容易被作为第一课的那类方法。

这篇文章我会从零手写一个 Python 版的牛顿-柯特斯求积实现，把 Cotes 系数的计算、单段公式、复合公式、误差表现、典型坑位全部过一遍。适合刚学数值分析、或者要在项目里快速实现数值积分的同学参考。读完你能直接抄代码去跑自己的积分问题，也能搞清楚为什么梯形公式精度差、辛普森公式是“性价比之王”、高阶公式为什么容易翻车。

1. 牛顿-柯特斯公式到底在解决什么问题

1.1 数值积分的基本困境

求定积分∫_a^b f(x)dx的原始思路是找原函数，然后用牛顿-莱布尼茨公式。但真实世界里，大量被积函数没有初等原函数，哪怕有，算起来也很痛苦。数值积分的思想就一句话：用一个容易积分的函数去近似被积函数，然后对这个近似函数做精确积分。

牛顿-柯特斯方法选的“容易积分的函数”是插值多项式。把区间[a, b]分成n等份，节点取x_i = a + i*h（h = (b-a)/n），用这些节点构造n次拉格朗日插值多项式P_n(x)去逼近f(x)，再用P_n(x)的积分替代原积分。这就是牛顿-柯特斯求积公式的基本框架。

公式写出来是：

∫_a^b f(x)dx ≈ (b-a) * Σ_{i=0}^{n} C_i * f(x_i)

其中C_i就是 Cotes 系数，只和n、i有关，和具体的a、b、f(x)都无关。这个性质很关键：系数可以先算好存起来，换不同积分区间和不同函数时直接复用。

1.2 用插值换积分：公式怎么来的

很多人背公式时只记住了(b-a)/6 * (f(a) + 4f(c) + f(b))这种辛普森形式，却不知道它怎么冒出来的。其实过程并不复杂。

先做一个坐标变换，把[a, b]映射到[0, n]。令x = a + t*h，那么dx = h*dt。于是：

∫_a^b f(x)dx = h * ∫_0^n f(a + t*h) dt

对每个节点j = 0,1,...,n，构造拉格朗日插值基函数：

L_j(t) = Π_{k≠j} (t - k) / (j - k)

然后用Σ f(a + j*h) * L_j(t)近似f(a + t*h)，积分后就得到：

C_j = (1/n) * ∫_0^n L_j(t) dt

所以 Cotes 系数本质上是拉格朗日基函数在[0, n]上的平均积分。这个推导看起来抽象，但一旦写成 Python 代码，就是几个循环的事。我建议你亲手推一遍 n=1 和 n=2 的情况，梯形公式和辛普森公式的系数会非常直观地浮现出来。

注意：上面说的是闭型公式，也就是积分节点包含端点a和b。还有一种开型公式不取端点，用于处理端点处被积函数奇异的情况。本文默认讨论闭型。

2. 常见形式与选型思路

2.1 梯形、辛普森与高阶公式

牛顿-柯特斯公式按n区分，最常用的几个是：

这里有个重要概念：代数精度。如果一个求积公式对所有次数不超过m的多项式都能精确成立，就说它有m次代数精度。

梯形公式(b-a)/2 * (f(a)+f(b))是一次插值多项式的精确积分，所以至少 1 次代数精度。辛普森公式用的是二次插值多项式，按理说只有 2 次代数精度，但因为它对三次多项式也精确——这个“白赚”的一阶精度让辛普森公式在工程中特别受欢迎。

偶数阶牛顿-柯特斯公式有个普遍规律：n=2k时，公式实际代数精度是2k+1，比插值多项式本身的次数还高一阶。所以n=2的辛普森和n=3的辛普森 3/8 都是 3 次代数精度，但 3/8 公式多一个节点，精度没有本质提升，一般只在区间长度不能整除成偶数份的时候用它补齐。

2.2 为什么工程中很少用高阶节点

一个很自然的想法是：既然n=2比n=1准，n=3比n=2准，那我一直加大n不就行了吗？理论上确实如此，但实际操作中高阶牛顿-柯特斯公式有两个致命问题。

第一个问题是 Cotes 系数会出现负值。从n=8开始，部分系数变成负数。负系数意味着求和时会出现大数相减，浮点舍入误差会被显著放大。数值分析里有句话叫“绝对误差下界与系数绝对值之和成正比”，当系数绝对值之和很大时，即使被积函数计算得很准，积分结果的误差也可能失控。

第二个问题是等距插值本身的局限性。插值节点等距分布时，多项式在高次情况下会在边界附近产生龙格振荡，也就是插值多项式剧烈摆动，导致牛顿-柯特斯公式对某些函数不仅不收敛，误差反而越来越大。

实际工程里的经验法则是：闭型牛顿-柯特斯公式做到n=4或n=5就行，再高就得不偿失。真正需要高精度时，优先用复合求积（把区间切细、每段用低阶公式），而不是盲目提高单段公式阶数。

3. Python实现：从系数到完整代码

3.1 最核心的一步：Cotes系数计算

写牛顿-柯特斯公式的 Python 实现，最关键的函数就是计算 Cotes 系数。我用的方法是按定义直接算：构造拉格朗日基函数，符号积分，除以n。

import numpy as np import sympy as sp def newton_cotes_coeffs(n): """ 计算闭型牛顿-柯特斯公式的 Cotes 系数。 原理： 把区间映射到 [0, n]，节点取 j = 0,1,...,n， 对每个节点构造拉格朗日基函数 L_j(t)， C_j = (1/n) * ∫_0^n L_j(t) dt """ t = sp.symbols('t') coeffs = [] for j in range(n + 1): basis = 1 for k in range(n + 1): if k != j: basis *= (t - k) / (j - k) val = sp.integrate(basis, (t, 0, n)) / n coeffs.append(float(val)) return np.array(coeffs)

这个实现很短，但已经把“插值-积分”的推导过程完整翻译成了代码。你可以先跑一下n=1和n=2，会看到输出分别是[0.5, 0.5]和[0.16666667, 0.66666667, 0.16666667]，和教材里完全一致。

用符号积分的好处是清晰、不容易错。缺点是在n较大时sympy.integrate会变得很慢。如果你只是写工具函数，平时最多用到n=6，这个速度完全能接受。如果对性能有强迫症，可以一次性把常用系数硬编码成查表格式，不过那就失去了推导过程的学习意义。

3.2 完整积分函数与复合形式

有了系数，求积函数就很简单了：

def newton_cotes(f, a, b, n): """ 单段 n 阶牛顿-柯特斯求积。 f: 被积函数 a, b: 积分区间 n: 等分数目 """ coeffs = newton_cotes_coeffs(n) xs = np.linspace(a, b, n + 1) ys = [f(x) for x in xs] return (b - a) * np.dot(coeffs, ys)

这里我用列表推导式逐点调用f(x)，而不是直接把整个数组传给f。原因是很多被积函数是普通 Python 函数，只接受标量输入；列表推导式对兼容两边的情况最稳。如果确认被积函数支持向量化运算，例如numpy函数，也可以直接写f(xs)，速度会更快。

接下来是复合求积。思路很简单：把[a, b]切成m段，每段用n阶公式，再把结果累加。

def composite_newton_cotes(f, a, b, n, m): """ 复合牛顿-柯特斯求积。 区间 [a, b] 分成 m 等份，每个子区间使用 n 阶公式。 """ sub = (b - a) / m total = 0.0 for i in range(m): left = a + i * sub right = left + sub total += newton_cotes(f, left, right, n) return total

这个函数虽然简单，但实际使用价值很高。比如辛普森公式做复合，就是工程里最经典的“复合辛普森求积”。也可以用更高效的方式写，把公共节点合并，避免重复计算端点值：

def composite_simpson_fast(f, a, b, m): """ 高效的复合辛普森公式，m 个等距子区间。 注意：m 个区间需要 2m+1 个节点。 """ h = (b - a) / (2 * m) xs = np.linspace(a, b, 2 * m + 1) ys = [f(x) for x in xs] return h / 3 * (ys[0] + ys[-1] + 4 * sum(ys[1:-1:2]) + 2 * sum(ys[2:-2:2]))

复合公式才是解决实际问题的正路。单段公式哪怕阶数选得再合适，面对变化剧烈的函数也容易歇菜；把区间切细之后，每段上的函数变化平缓，低阶公式也能保持很好的精度。

3.3 用分数核对系数的小技巧

用sympy计算系数时，默认会把它转成浮点数。但学习阶段，我更建议你顺手打印出分数形式，方便和教材对照：

def cotes_coeffs_fraction(n): t, j, k = sp.symbols('t j k') coeffs = [] for j in range(n + 1): basis = 1 for k in range(n + 1): if k != j: basis *= (t - k) / (j - k) coeffs.append(sp.simplify(sp.integrate(basis, (t, 0, n)) / n)) return coeffs

跑一下n=4，你会看到输出是[7/90, 16/45, 2/15, 16/45, 7/90]这样的精确分数。一旦在某个n看到系数出现负数，就要提起警觉：这个阶数不适合上机使用。

4. 数值实验：公式到底准不准

4.1 三个典型测试函数

光说不练假把式。我选了三个有解析真值的函数做测试：

1. f(x) = e^x，[0, 1]，真值e - 1 ≈ 1.718281828459045
2. f(x) = sin(x)，[0, π]，真值2
3. f(x) = 1/(1+x^2)，[0, 1]，真值π/4 ≈ 0.7853981633974483

测试代码就一行一行跑：

import math f1 = lambda x: math.exp(x) true1 = math.e - 1 for n in range(1, 7): approx = newton_cotes(f1, 0, 1, n) err = true1 - approx print(f"n={n}, approx={approx:.15f}, err={err:.3e}")

我在本机跑出来的部分结果如下，不同环境下最后几位可能有细微差异，但趋势是稳定的。

对e^x这种平滑函数，单段牛顿-柯特斯的收敛速度非常恐怖。n=5直接逼近了机器精度。这也是为什么教材总拿e^x举例——它实在是太“配合”了。

但换到sin(x)在[0, π]上的积分，画风突变：

n=1的梯形公式直接算出 0，因为sin(0)和sin(π)都是 0，公式只用了两端点，完全丢失了内部的“山峰”。这暴露了单段公式的盲区：节点少的时候，积分结果对函数内部的剧烈变化毫无感知。解决办法就是复合——把区间切成若干段，让每个子区间内部的峰谷都被照顾到。

第三个函数1/(1+x^2)是龙格现象的标准受害者。它在[0,1]上其实挺平滑，但用等距节点逼近它时，高次插值的边界振荡会悄悄给积分结果带来误差。单段n=4的布尔公式误差大约在1e-7量级，看起来还行；但一旦把区间拉大到[-5,5]，高次公式的误差就会明显上升。这也提醒我：测试数值积分算法，别只用单调光滑函数，得多换几种形态。

4.2 误差与收敛性分析

牛顿-柯特斯公式的误差和n、被积函数的光滑性密切相关。梯形公式的误差项是O(h^3)量级，辛普森公式是O(h^5)。复合时，如果把区间等分成m段，每段长度h=(b-a)/m，那么复合梯形误差是O(h^2)，复合辛普森误差是O(h^4)。

这个“每复合一次，误差阶数提升”的规律，可以很直观地验证。我用复合梯形和复合辛普森算同一个积分∫_0^1 e^x dx，把m从 4 增加到 8，看误差缩小的比例。复合梯形误差大约缩小 4 倍，复合辛普森误差大约缩小 16 倍。这个 4 倍和 16 倍不是巧合，正是h^2和h^4收敛阶的体现。

写代码验证收敛阶的时候，有个细节值得注意：用“误差比”而不是“误差绝对值”来判断收敛速度更直观。比如上一步误差是E_m，下一步把区间加密一倍，如果E_{2m} / E_m ≈ 1/16，就说明当前算法是 4 阶收敛。这个现象比单纯报一个误差数字更有说服力，也更容易在项目汇报里展示。

4.3 单点积分的“舒适区”

从实验里能总结出一个大致的选型准则：

• 被积函数平滑、区间不大、对精度要求一般：n=2复合辛普森，无脑好用。
• 被积函数有振荡但周期不明显，精度要求高：增加复合段数，而不是增加 n。
• 被积函数在端点有奇异：换成开型牛顿-柯特斯或高斯求积。
• 积分区间非常大：先考虑变换或分段，别指望一个单段公式硬扛。

我个人在实际工程里，90% 的数值积分需求是复合辛普森解决的。不是因为其他方法不好，而是它代码量小、行为稳定、误差可预测。

5. 实战中的坑与排查技巧

5.1 遇到的主要问题与对策

再优雅的公式，到了实际代码里都会碰到各种意料之外的情况。下面是我实践中整理的一份问题速查表。

最典型的高阶翻车场景我再说一下。用n=8的牛顿-柯特斯公式算某个看起来正常的光滑函数积分，结果可能比n=4还差，甚至出现负值。问题不在算法本身，而在于系数里的负号。比如n=8时部分 Cotes 系数绝对值很大且正负交替，浮点计算时细微的舍入误差被放大了好几个数量级。所以看到高阶公式，第一反应不应该是“更高阶等于更精确”，而是“这个系数表里有没有负号”。

5.2 容易被忽略的精度细节

第一个容易被忽略的细节是：端点处的函数值计算。牛顿-柯特斯公式的权重对端点非常敏感。如果f(a)或f(b)计算本身就有较大误差，这个误差会直接进入积分结果。极端情况下，f(a)是两个接近的大数相减得到的，那积分结果的误差可能比公式的理论误差还要大。这一点在处理来自实验数据的分段函数时尤其要小心。

第二个细节是：复合公式不要重复计算公共节点。很多新手写复合辛普森，会用循环对每个子区间单独调用辛普森公式，这样相邻区间共享的端点被重复算了多次，函数调用次数多了一倍。虽然结果数值上不出错，但性能白白损失。我自己写程序时有一个习惯：能用数组切片就不要循环逐个加，能用一次f(xs)就不要反复调用几十次 Python 函数。

第三个细节是：区间尺度。公式里的(b-a)是整体乘上去的。如果b-a非常大，比如1000，而系数和函数值的量级都很小，最终结果可能出现有效位数丢失。遇到这种情况，最好做变量替换把区间缩放到一个合理尺度，比如归一化到[0, 1]，算完再乘回缩放因子。这个操作对浮点精度的影响，经常比换一个高阶公式还明显。

第四个细节是：误差估计。工程上不要只看一个结果就交付，至少要比较两个不同加密程度的计算结果。比如用m=8和m=16的复合辛普森各算一次，如果两者在小数点后 6 位一致，基本可以认为结果已经收敛到位。如果差别还大，就继续加密，直到结果稳定。这个方法虽然土，但比任何误差公式都可靠。

第五个细节是：周期性函数的特殊福利。如果被积函数是周期函数，且积分区间刚好是周期的整数倍，梯形公式的收敛速度会异常快，快到和辛普森公式不相上下。这是 Euler-Maclaurin 公式背后的一个结论。所以遇到周期函数，可以直接用复合梯形，不要觉得它“太低阶”就嫌弃。我试过用梯形公式算某个周期函数的积分，加密几次后精度惊人地高，完全颠覆了对梯形公式的刻板印象。

我个人在实际操作中的体会是，牛顿-柯特斯公式是理解数值积分的最佳入门模型，但它不是万能的。真正干活时，我会观察被积函数的形态，再决定用单段公式还是复合公式。日常需求里，复合辛普森是默认选项；遇到振荡特别厉害的函数，直接切到高斯求积或者自适应积分。

最后再分享一个小技巧。写完牛顿-柯特斯代码后，别急着拿去算复杂函数，先用一组已知真值的多项式测试。比如f(x)=x^3+x^2+x+1在[0,1]积分，如果辛普森公式没有精确返回1.75，说明代码里系数算错了。这个自测习惯能省下大量排查时间。做数值计算的程序，第一件事永远是验证可信度，而不是追求高精度。