4、数值积分基础：辛普森法则（1/3法则与3/8法则）的原理与推导

聊到数值积分，很多新手第一反应就是梯形法。梯形法简单，但精度嘛……说实话，有时候真不够用。我记得刚入行那会儿，用梯形法算一个流体力学中的压力积分，结果跟实验数据差了将近5%。导师看了一眼就说：「你试试辛普森法则。」

嗯，从那以后，我对辛普森法则就特别上心。今天咱们就来掰扯掰扯这个经典方法——辛普森1/3法则和3/8法则。

4.1 为什么需要更高阶的积分方法？

先想一个问题：梯形法用直线段近似曲线，那如果我用抛物线呢？

说白了，梯形法是用一次多项式（直线）去拟合被积函数。但很多实际问题的被积函数是弯曲的，直线拟合自然会有误差。辛普森法则的核心思想就是：用二次或三次多项式来近似被积函数，从而获得更高的精度。

我在做雷达信号处理仿真时，遇到过一种情况：信号功率谱密度曲线非常陡峭，梯形法需要划分上千个区间才能收敛。换成辛普森1/3法则后，区间数直接砍了一半，精度反而更高。这就是高阶方法的魅力。

4.2 辛普森1/3法则的原理与推导

辛普森1/3法则，名字听着挺玄乎，其实原理很简单：用二次抛物线来近似被积函数。

假设我们要计算：

I = ∫[a, b] f(x) dx

取三个点：左端点 a、中点 m = (a+b)/2、右端点 b。用这三个点构造一条二次抛物线，然后求抛物线下的面积。

推导过程是这样的：

设二次多项式为：

P(x) = Ax² + Bx + C

它必须满足：

P(a) = f(a)
P(m) = f(m)
P(b) = f(b)

解出A、B、C后，对P(x)从a到b积分。嗯，这里我跳过繁琐的代数运算，直接给结果：

∫[a, b] f(x) dx ≈ (b-a)/6 · [f(a) + 4f(m) + f(b)]

你可能会问：为什么叫「1/3」法则？因为误差项里有个1/3的系数。具体来说，这个方法的截断误差是：

E = -(b-a)^5 / 2880 · f⁽⁴⁾(ξ)

注意，这里要求f(x)有四阶连续导数。如果被积函数是三次或更低次的多项式，辛普森1/3法则能给出精确结果。这一点我在做多项式拟合验证时亲自测试过，确实如此。

4.3 复合辛普森1/3法则

实际工程中，积分区间往往很大，直接用一次辛普森1/3法则精度不够。这时候就需要「复合」——把区间分成多个小段，每段用辛普森1/3法则。

具体做法：

1. 将区间[a, b]分成n个等长子区间，n必须是偶数
2. 每两个子区间（即三个点）应用一次辛普森1/3法则
3. 把所有结果加起来

公式长这样：

∫[a, b] f(x) dx ≈ h/3 · [f(x₀) + 4f(x₁) + 2f(x₂) + 4f(x₃) + ... + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)]

其中 h = (b-a)/n，xᵢ = a + i·h。

注意系数规律：首尾是1，奇数下标是4，偶数下标是2。我曾经在写代码时搞反了奇偶顺序，结果算出来的积分值完全不对。排查了半天才发现这个低级错误……嗯，这种细节真的要注意。

4.4 辛普森3/8法则

辛普森1/3法则用三个点（二次多项式），那3/8法则呢？它用四个点（三次多项式）。

公式如下：

∫[a, b] f(x) dx ≈ (b-a)/8 · [f(a) + 3f(x₁) + 3f(x₂) + f(b)]

为什么叫「3/8」？因为系数是3/8。这个方法的截断误差是：

E = -(b-a)^5 / 6480 · f⁽⁴⁾(ξ)

对比一下：1/3法则的误差系数是1/2880 ≈ 0.000347，3/8法则的系数是1/6480 ≈ 0.000154。3/8法则的误差更小，但代价是需要多一个采样点。

从表格能看出来，1/3法则和3/8法则的误差阶数相同，都是O(h⁴)。但3/8法则的系数更小，精度略高。不过在实际工程中，1/3法则更常用，因为它只需要偶数个区间，实现起来更灵活。

4.5 两种法则的对比与选择

我个人习惯这样选：

• 区间数n是偶数：用复合辛普森1/3法则，简单高效
• 区间数n是3的倍数：用复合辛普森3/8法则，精度略高
• 区间数n既不是偶数也不是3的倍数：混合使用。比如n=5时，前两个区间用3/8法则，后三个区间用1/3法则

4.6 代码实现示例

下面是一个简单的Python实现，展示了两种法则的用法：

def simpson_13(f, a, b):
    """辛普森1/3法则，单区间"""
    m = (a + b) / 2
    return (b - a) / 6 * (f(a) + 4*f(m) + f(b))

def composite_simpson_13(f, a, b, n):
    """复合辛普森1/3法则，n必须为偶数"""
    if n % 2 != 0:
        raise ValueError("n必须是偶数")
    h = (b - a) / n
    x = [a + i*h for i in range(n+1)]
    result = f(x[0]) + f(x[-1])
    for i in range(1, n):
        if i % 2 == 0:
            result += 2 * f(x[i])
        else:
            result += 4 * f(x[i])
    return result * h / 3

def simpson_38(f, a, b):
    """辛普森3/8法则，单区间"""
    h = (b - a) / 3
    x1 = a + h
    x2 = a + 2*h
    return (b - a) / 8 * (f(a) + 3*f(x1) + 3*f(x2) + f(b))

你看，代码其实不复杂。关键是要理解背后的原理，这样在遇到边界情况时才知道怎么处理。

4.7 知识体系总览

为了让你更直观地理解本章的知识结构，我画了一张图：

这张图把两种法则的核心要素都串起来了。你想想看，掌握了这些，你在仿真中遇到积分问题，是不是心里就有底了？

好了，辛普森法则的原理和推导就聊到这儿。记住，理论是基础，但真正让你成长的是在实践中不断踩坑、填坑的过程。下次遇到数值积分问题，不妨先试试辛普森法则——它很可能比梯形法更对得起你的仿真时间。

公众号:蓝海资料掘金营，微信deep3321