快速数值积分
4.89/5 (64投票s)
使用最优算法对有限区间内的光滑函数进行数值积分。
引言
在实际应用中遇到的绝大多数积分都必须通过数值方法进行求值。不幸的是，没有一种通用代码能够高效准确地计算任何一个给定的积分。然而，我们可以编写针对特定问题类型的积分例程，使其能很好地工作。
本文介绍的方法是对有限区间内的光滑函数进行积分。从某种意义上说，该方法具有最优效率，在给定函数评估次数的情况下，能达到最高的精度。
更准确地说，本文介绍的方法假定被积函数是解析的。这意味着函数要求非常光滑，在函数或其任何导数中都没有不连续点。这排除了使用“if”语句拼接起来的函数，例如“如果 x 小于 0，则返回 x³，如果 x 大于 0，则返回 x²”。这也排除了在积分区间中间取无穷大的函数。另一方面，它允许在积分区间的端点取无穷大的函数。
本文介绍的方法是双指数变换。有关数学细节的解释，请参阅 Masatake Mori 和 Masaaki Sugihara 的文章《数值分析中的双指数变换》，发表于《计算与应用数学杂志》，第 127 卷（2001 年），页码 287-296。该方法背后的数学相当复杂，此处不再赘述。然而，实现该方法的代码简洁易用。
Using the Code
积分器实现为一个模板类 DEIntegrator。（名称中的“DE”代表代码核心的“双指数”算法。）模板参数是被积函数的类。代码使用函数类而不是函数指针，因为前者更灵活。函数通常带有需要附加的参数，而在 C++ 中，函数类是处理这种情况的方式。（在函数式语言中，您可能会使用闭包。）
要在代码中使用 DEIntegrator，您需要将两个文件添加到项目中：DEIntegrator.h 和 DEIntegratorConstants.h，并在您想要调用积分器的文件中添加 #include "DEIntegrator.h"。
DEIntegrator 的 Integrate 方法有两个重载版本。第一个是最容易使用的。
static double Integrate
(
const TFunctionObject& f, // [in] integrand
double a, // [in] left limit of integration
double b, // [in] right limit of integration
double targetAbsoluteError // [in] desired bound on error
)
您只需传入要积分的函数、积分区间以及您希望计算结果的精度。返回值就是积分的值。
另一个重载版本允许您在积分计算完成后获取更多信息。
static double Integrate
(
const TFunctionObject& f, // [in] integrand
double a, // [in] left limit of integration
double b, // [in] right limit of integration
double targetAbsoluteError, // [in] desired bound on error
int& numFunctionEvaluations, // [out] number of function evaluations used
double& errorEstimate // [out] estimated error in integration
)
另外两个参数让您知道计算积分所需的函数评估次数以及算法的误差估计。请注意，targetAbsoluteError 指定您希望误差有多小，而 errorEstimate 估计实际误差有多小。通常，后者会远小于前者，即该方法比您要求的做得更好。然而，如果该方法在处理您的被积函数时遇到困难（例如，由于端点处存在强奇异点），那么 errorEstimate 可能会大于 targetAbsoluteError，从而告知您存在问题。
示例 1
假设我们要对函数 f(x) = exp(-x/5) (2 + sin(2x)) 在 0 到 10 的区间上进行积分，如下图所示：
事实证明，该积分可以精确计算。其值为 9.1082396073230。这是一个有些人工的例子：通常，我们使用数值积分来处理那些无法精确计算的积分。但是，精确积分对于说明和测试很有用。假设我们希望将积分的数值结果保留到小数点后六位。
首先，我们编写函数对象：
class DemoFunction
{
public:
double operator()(double x) const
{
return exp(-x/5.0)*(2.0 + sin(2.0*x));
}
};
然后，我们如下计算其积分：
DemoFunction f;
int evaluations;
double errorEstimate;
double integral = DEIntegrator<DemoFunction>::Integrate
(f, 0, 10, 1e-6, evaluations, errorEstimate);
std::cout << integral << ", " << errorEstimate
<< ", " << evaluations << "\n";
输出是：
9.10823960732284, 5.21017118337852e-009, 97
因此，在这个例子中，我们要求误差小于 10-6，而该方法估计误差小于 5 × 10-9，实际上误差更接近 10-13。该方法是“少说多做”。在此过程中，它使用了 97 次函数评估。
示例 2
作为另一个例子，我们考虑函数 f(x) = x-1/3(1-x)5，如下图所示。我们要对该函数在 0 到 1 的区间上进行积分。这个例子是为了说明被积函数不必在所有地方都是解析的，只需在积分区间的内部是解析的即可。换句话说，如果函数在至少一个端点处发散，也是可以接受的。
我们如下实现被积函数的函数对象：
class DemoFunction2
{
public:
double operator()(double x) const
{
return pow(1.0 - x, 5.0)*pow(x, -1.0/3.0);
}
};
现在，假设我们只想要估计的积分值。我们不关心需要多少函数评估次数或估计误差。在这种情况下，计算积分的代码更加简洁。
DemoFunction2 f2;
integral = DEIntegrator<DemoFunction2>::Integrate(f2, 0, 1, 1e-6);
std::cout << integral << "\n";
在这种情况下，该积分碰巧也是一个可以精确计算的积分。实际上，积分值为 2187/5236。软件计算出的积分值为 0.41768525570112。误差约为 2 × 10-10。
误差估计
基本的数值积分中的误差会随着积分点数量的增加而以某个常数幂的速率下降。例如，辛普森法则的误差以 N-4 的速率减小，其中 N 是积分点数量。但是，本文介绍的方法中的误差呈指数级下降。具体来说，误差的量级为 exp(-cN/log(N))。有关数学细节，请参阅此页面。
在内部，该代码通过比较连续迭代之间的差异与理论预测的差异来监控自身的误差率。这使得用户可以指定误差容差，并且代码事后可以估计其表现如何。
讨论
本文讨论的方法能够准确高效地对有限区间内的解析函数进行积分。从技术意义上讲，该方法对于它所集成的函数类别来说是尽可能高效的，具体描述在此。
与大学课程中介绍的大多数方法不同，本文介绍的方法不假设被积函数的行为大致像多项式。被积函数可能在端点处是奇异的，如示例 2 所示，而不会损害方法的准确性。此外，与更基本的方法不同，收敛速率是指数级的而不是多项式的。这意味着本文介绍的方法收敛速度非常快。
双指数规则有用于无界区域积分的变体。不幸的是，这些通常难以实现；软件需要了解有关被积函数的更多详细信息。一种更简单的方法是将积分问题转化为可以使用本文软件计算的问题。最明显的方法是将无界积分截断为有界积分。另一种方法是使用变量替换将积分转换为有界区间上的新积分。
不建议截断无界区域上的积分。决定截断点可能很困难。选择过小的区间可能导致答案不准确。选择过大的区间可能导致浪费时间对被积函数几乎为零的区域进行积分。
使用变量替换将积分转换可能更有效和准确。对于 (0, ∞) 上的积分，一种常见的变量替换是 x = t/(1-t)。这会将 f(x) 从 0 到 ∞ 的积分转换为 f(t/(1-t)) /(1-t)2 从 0 到 1 的积分。对于 (-∞, ∞) 上的积分，变量 x = tan(t) 将 f(x) 从 -∞ 到 ∞ 的积分转换为 f(tan(t)) (1 + sec2(t)) 从 -1 到 1 的积分。这些变量替换通常效果很好。通过为您的积分问题量身定制转换，有可能获得更好的性能，但这通常不是必需的。
历史
- 2008 年 12 月 5 日：初始版本。
- 2008 年 12 月 8 日：添加了对误差估计的讨论。
- 2009 年 6 月 29 日：添加了对无界区域积分的讨论。
- 2015 年 5 月 4 日：修复了处理目标绝对误差中的错误。