数值积分

引言

若函数在 $f(x)$ 在区间 $[a,b]$ 上连续，且原函数为 $F(x)$，则有 Newton-Leibniz 公式 $\int_a^b f(x)\mathrm d x=F(b)-F(a)$。在实际情况中，$f(x)$ 的原函数往往过于复杂甚至难以求解，我们不得不用数值积分的方式得到定积分的近似解。

数值积分就是对于 $f(x)$，选择合适的离散点 $a\le x_0<x_1<\cdots<x_n\le b$，并给离散点的函数值 $f(x_k)$ 合适的权值 $A_k$，从而构造出近似解： $\int_a^b f(x)\mathrm dx\approx\sum\limits_{k=0}^n A_kf(x_k)$。记 $I(f)=\int_a^b f(x)\mathrm dx$，$Q(f)=\sum\limits_{k=0}^n A_kf(x_k)$，则截断误差（余项）$R(f)=I(f)-Q(f)$。

数值积分方法

Newton-Cotes 求积公式

选择 $n+1$ 个离散点 $a\le x_0<x_1<\cdots<x_n\le b$，用 Lagrange 插值多项式近似表示 $f(x)$，则 $Q(f)=\displaystyle\int_a^b\left[\sum\limits_{k=0}^nl_k(x)f(x_k)\right]\mathrm dx=\sum\limits_{k=0}^n\left[\int_a^b l_k(x)\mathrm dx\right]f(x_k)$，其中 $A_k=\int_a^bl_k(x)\mathrm dx$。截断误差 $E(f)=I(f)-Q(f)=\displaystyle\int_a^b\left[f(x)-\sum\limits_{k=0}^nl_k(x)f(x_k)\right]\mathrm dx=\int_a^b\frac{f^{(n+1)}(\xi)}{(n+1)!}\omega_{n+1}\mathrm dx$。

令 $x_0\sim x_n$ 等距分布，根据采样点数量 $n+1$ 的不同，依次对应不同的求积公式。

代数精度

若对任一不超过 $m$ 次的多项式 $p_m(x)$，Newton-Cotes 求积公式总有 $I(p_m)=Q(p_m)$，而对某一个 $m+1$ 次多项式 $p_{m+1}(x)$，$I(p_{m+1})\neq Q(p_{m+1})$，则称此求积公式的代数精度为 $m$。

要验证求积公式的代数精度为 $m$，只需要验证对于 $f(x)=1,x,\cdots,x^m$，$I(f)=Q(f)$ 成立，对于 $f(x)=x^{m+1}$，$I(f)\neq Q(f)$。

有结论：$n$ 次插值型求积公式的代数精度至少为 $n$ 次；反之，如果求积公式 $\int_a^b f(x)\mathrm dx\approx\sum\limits_{k=0}^n A_kf(x_k)$ 的代数精度至少为 $n$ 次，则它必定是插值型的。

梯形求积公式

取 $n=1$，则 $x_0=a$，$x_1=b$，求积公式为 $Q(f)=\frac{b-a}{2}[f(a)+f(b)]$，实际上就是用梯形面积代替积分。余项 $E_1=\frac{f''(\xi)}{2!}\int_a^b(x-a)(x-b)\text dx=-\frac{1}{12}f''(\eta)(b-a)^3$，$\eta\in(a,b)$。

Simpson 求积公式

取 $n=2$，求积公式 $Q(f)=\frac{b-a}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]$。

Simpson $\frac{3}{8}$ 求积公式

取 $n=3$，求积公式 $Q(f)=\frac{b-a}{8}\left[f(a)+3f\left(\frac{2a+b}{3}\right)+3f\left(\frac{a+2b}{3}\right)+f(b)\right]$。

Cotes 求积公式

取 $n=4$，求积公式 $Q(f)=\frac{b-a}{90}\left[7f(a)+32f\left(\frac{3a+b}{4}\right)+12f\left(\frac{a+b}{2}\right)+32f\left(\frac{a+3b}{4}\right)+7f(b)\right]$。

复化求积公式

随着选取的离散点的个数增多，Newton-Cotes 公式的精度也会越高。但是对于 $n\ge 7$ 时的 Newton-Cotes 公式，它的各项系数开始出现负值。根据误差理论的研究，当积分公式开始出现负值时，可能导致舍入误差较大，因此通过增加求积节点数来提高计算精度并不是一种好方法。

不妨将积分区间划分成若干个小区间，在每个小区间上应用低次 Newton-Cotes 公式，然后把所有小区间上的计算结果累加得到整个区间上的积分。

复化梯形公式

将 $[a,b]$ 分成 $n$ 等分，其中节点 $x_i=a+ih$，$h=\frac{b-a}{n}\ (i=0,1,\cdots,n)$。

$\displaystyle I=\int_a^b f(x)\text dx=\sum\limits_{k=0}^{n-1}\int_{x_{k}}^{x_{k+1}}f(x)\text dx\approx\sum\limits_{k=0}^{n-1}\frac{h}{2}\left[f(x_k)+f(x_{k+1})\right]=\frac{h}{2}\left[f(a)+\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]$，即 $\displaystyle T_n=\frac{h}{2}\left[f(a)+\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]$。余项 $R_n(f)=-\frac{b-a}{12}h^2f''(\eta)$，$\eta\in(a,b)$。

复化Simpson公式

推导过程与复化梯形公式类似，$\displaystyle S_n=\frac{h}{6}\left[f(a)+4\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})+2\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]$。

步长 $h$ 的选取

通常采取将区间不断对分的方法，即 $n=2^k$；反复使用复化求积公式，直到相邻两次计算结果之差的绝对值小于指定精度 $\varepsilon$。

以复化梯形公式举例，首先将 $[a,b]$ 分成 $n$ 等分，步长为 $h=\frac{b-a}{n}$。$\displaystyle T_n=\sum\limits_{k=0}^{n-1} \frac{h}{2}\left[f(x_k)+f(x_{k+1})\right]=\frac{h}{2}\left[f(a)+\sum\limits_{k=1}^{n-1}f(x_k)+f(b)\right]$，$\displaystyle T_{2n}=\sum\limits_{k=0}^{n-1}\frac{h}{4}\left\{\left[f(x_k)+f(x_{k+\frac{1}{2}})\right]+\left[f(x_{k+\frac{1}{2}})+f(x_{k+1})\right]\right\}=\frac{h}{4}\sum\limits_{k=0}^{n-1}\left[f(x_k)+f(x_{k+\frac{1}{2}})\right]+\frac{h}{2}\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})=\frac{1}{2}T_n+\frac{h}{2}\sum\limits_{k=0}^{n-1}f(x_{k+\frac{1}{2}})$。

由复化梯形公式的余项知，$I-T_n=-\frac{b-a}{12}h^2f''(\eta_1)$，$I-T_{2n}=-\frac{b-a}{12}\left(\frac{h}{2}\right)^2f''(\eta_2)$；假设 $f''(x)$ 变化不大，即 $f''(\eta_1)\approx f''(\eta_2)$，得到 $\frac{I-T_n}{I-T_{2n}}\approx4$，即 $I\approx T_{2n}+\frac{1}{3}(T_{2n}-T_n)$。当 $T_{2n}-T_n<\varepsilon$，有 $I-T_{2n}<\frac{\varepsilon}{3}$，因此不断将区间对分的方法是可行的。

Romberg 求积公式

在上面的推导中，$I\approx T_{2n}+\frac{1}{3}(T_{2n}-T_n)$，可认为 $\overline{T}=T_{2n}+\frac{1}{3}(T_{2n}-T_n)$ 是一个较 $T_{2n}$ 更精确的值。实际上，将 $T_{2n}$ 用 $T_n$ 表示，再将 $T_n$ 展开，可以得到 $S_n=\overline{T}$；即 $S_n=\frac{4}{3}T_{2n}-\frac{1}{3}T_n$，类似的有 $C_n=\frac{16}{15}S_{2n}-\frac{1}{16}S_n$，$R_n=\frac{64}{63}C_{2n}-\frac{1}{63}C_n$。

称 $R_n$ 为 Romberg 公式，用 Romberg 方法求数值积分类似填表的过程。

梯形序列	Simpson序列	Cotes序列	Romberg序列
$T_1$
$T_2$	$S_1$
$T_4$	$S_2$	$C_1$
$T_8$	$S_4$	$C_2$	$R_1$
$T_{16}$	$S_8$	$C_4$	$R_2$

用 $\displaystyle I=\int_a^b \frac{4}{1+x^2}\text dx$ 且 $\varepsilon=10^{-4}$ 举例，Romberg 方法的代码实现如下。

#include<iostream>
#include<cmath>
using namespace std;
double ans[20][4];
inline double f(double x)
{
	return 4.0 / (1 + x * x);
}
int main()
{
	double a, b;
	// 输入积分上下限
	cin >> a >> b;
	double h = b - a;
	int n = 1;
	// 计算 T[1]
	ans[0][0] = (f(a) + f(b)) * h / 2.0;
	printf("%.15lf\n", ans[0][0]);
	// 填表
	for (int i = 1;;i++)
	{
		// 计算 T_{2n}
		ans[i][0] = 0.5 * ans[i - 1][0];
		for (int k = 0;k <= n - 1;k++)
		{
			ans[i][0] += h * 0.5 * f((a + k * h + a + (k + 1) * h) / 2.0);
		}
		// 更新步长
		n <<= 1;
		h = (b - a) / (double)n;
		printf("%.15lf ", ans[i][0]);
		
		if (i >= 1)
		{
			// Simpson
			ans[i][1] = ans[i][0] * 4 / 3.0 - 1.0 / 3.0 * ans[i - 1][0];
			printf("%.15lf ", ans[i][1]);
		}
		if (i >= 2)
		{
			// Cotes
			ans[i][2] = ans[i][1] * 16 / 15.0 - 1.0 / 15.0 * ans[i - 1][1];
			printf("%.15lf ", ans[i][2]);
		}
		if (i >= 3)
		{
			// Romberg
			ans[i][3] = ans[i][2] * 64 / 63.0 - 1.0 / 63.0 * ans[i - 1][2];
			printf("%.15lf", ans[i][3]);
		}
		if (i >= 4)
		{
			// 检查误差
			if (fabs(ans[i][3] - ans[i - 1][3]) < 0.0001)
			{
				return 0;
			}
		}
		putchar('\n');
	}
	return 0;
}