辛普森法则的数值积分

John D. Cook

4.97/5 (19投票s)

2012 年 2 月 22 日

BSD

4分钟阅读

70163

3036

如何实现和使用辛普森法则

引言

辛普森法则是一种简单而有效的数值积分技术。然而，实际应用需要比方法介绍中通常呈现的更多。本文将展示如何在 C# 中实现该算法，并将讨论如何准备积分进行评估。

背景

微积分书籍通常会快速介绍梯形法则和辛普森法则，作为近似计算无法手工精确计算的积分的方法。

给定一个区间 [a, b] 和一个偶数 n，辛普森法则近似于积分

$\int_a^b f(x) \, dx$

$\frac{h}{3} \left[f(x_0) + 4 f(x_1) + 2 f(x_2) + 4 f(x_3) + 2f(x_4) + \cdots + 4 f(x_{n-1}) + f(x_n) \right]$

其中 h = (b - a)/n 和 x_i = a + ih。

虽然这是对算法的正确描述，但它留下了两个相关的问题未解答

近似有多好？
我应该如何选择 n？

在实践中，您在开始之前很少知道 n 需要多大才能满足您的精度要求。如果您可以在运行时估计误差，则可以使用该估计来确定 n。

一种常见的方法是重复应用辛普森法则，在每个阶段将积分点的数量加倍。两个连续阶段之间的差异给出了误差的估计值，因此您继续加倍阶段，直到误差估计低于您的容差。

实现这种方法有几个实际问题。首先，如果天真地实现，该方法将多次在相同的点重新评估被积函数。在每个新阶段，一半的积分点与前一阶段相同。一些代数运算使得可以重用前一阶段的积分估计，并且仅在新积分点评估被积函数。

第二个实际实现问题是不要过早地检查停止条件。如果您在每个阶段之后检查停止规则，则在早期阶段，您的积分点可能会错过被积函数中的可变区域，并且您会被误导认为该方法已经收敛，而实际上并没有。此处给出的代码仅在辛普森法则的五个阶段后检查收敛性。

Using the Code

类 SimpsonIntegrator 只有两个方法，Integrate 的重载实现。 Integrate 的更简单版本只返回积分的值。更详细的版本返回使用的函数评估次数和误差估计。更简单的版本通过调用更详细的方法并丢弃额外信息来实现。此外，更详细的版本允许用户指定最大积分阶段数；更简单的版本为此参数设置默认值。

这是一个在区间 [0, 2π] 上积分函数 f(x) = cos(x) + 1 的示例

    double twopi = 2.0*Math.PI;
    double epsilon = 1e-8;
    double integral = SimpsonIntegrator.Integrate(x => Math.Cos(x) + 1.0, 0, 
        twopi, epsilon);
    double expected = twopi;
    Assert.IsTrue(Math.Abs(integral - expected) < epsilon);

请注意，C# 对隐式函数的支持对于简单的被积函数非常方便。

根据辛普森法则的理论推导，该方法应该仅用一个积分阶段就能精确地积分三次多项式。这是一个调用更详细接口的示例，以指定该方法应仅使用一个阶段来积分区间 [a, b] 上的 f(x) = x³。我们还通过检查被积函数仅被评估三次来验证代码是否尊重了我们的请求，即仅使用一个阶段。（辛普森法则的第一阶段评估积分在区间的两端和中点。）

    double a = 2, b = 5;
    double epsilon = 1e-10;
    double integral;
    int functionEvalsUsed;
    double estimatedError;
    double expected;
    int maxStages = 1;

    integral = SimpsonIntegrator.Integrate(x => x*x*x, a, b, epsilon, 
        maxStages, out functionEvalsUsed, out estimatedError);
    expected = (b*b*b*b - a*a*a*a)/4.0;
    Assert.AreEqual(functionEvalsUsed, 3);
    Assert.IsTrue(Math.Abs(integral - expected) < epsilon);

辛普森法则仅适用于闭区间上的平滑函数。如果这不符合您的积分描述，那么在稍作准备后仍然可以应用该方法。另一方面，直接在不适当的地方应用该方法可能会产生不良结果。

中间有扭结的函数不平滑，因此辛普森法则不直接适用。但是，可以将这样的函数简单地分为两个积分，其中被积函数在两个积分中分别平滑。例如，要积分 exp(-|x|) 在 [-3, 5] 上，我们需要将问题分解为 [-3, 0] 上的一个积分和 [0, 5] 上的另一个积分，因为被积函数在 0 处有一个急剧的弯曲。

辛普森法则不能直接计算无限范围上的积分。但是，通常可以使用变量替换将积分转换为有限范围内表现良好的积分。

辛普森法则也不适用于变为无穷大的函数。有关如何使用辛普森法则进行此类积分的示例，请参阅奇点的注意事项。