19.3.4 龙贝格方法

为提高数值积分的精度, 推荐从梯形求和序列出发, 重复对分积分步长的龙贝格 (Romberg) 方法.

19.3.4.1 龙贝格方法的算法

该方法包含以下步骤.

1. 确定梯形和

根据 19.3.2.2 的 (19.76) 确定积分 ${\int }_a^{b}f\left(x\right)\mathrm{d}x$ 关于步长

$$\begin{array}{}\text{(19.85)}& h_i=\frac{b-a}{2^i}\left(i=0,1,2,\cdots ,m\right)\end{array}$$

的近似梯形和 $T\left(h_i\right)$ ,这里考虑递归关系

$$T\left(h_i\right)=T\left(\frac{h_{i-1}}{2}\right)$$$$=\frac{h_{i-1}}{2}[\frac{1}{2}f\left(a\right)+f\left(a+\frac{h_{i-1}}{2}\right)+f\left(a+h_{i-1}\right)+f\left(a+\frac{3}{2}{h}_{i-1}\right)$$$$+f\left(a+2h_{i-1}\right)+\cdots +f\left(a+\frac{2n-1}{2}{h}_{i-1}\right)+\frac{1}{2}f\left(b\right)]$$$$=\frac{1}{2}T\left(h_{i-1}\right)+\frac{h_{i-1}}{2}\sum _{j=0}^{n-1}f\left(a+\frac{h_{i-1}}{2}+j{h}_{i-1}\right)\left(i=1,2,\cdots ,m;n=2^{i-1}\right).$$

(19.86)

递归公式 (19.86) 告诉我们,由 $T\left(h_{i-1}\right)$ 计算 $T\left(h_i\right)$ 仅需在新增插值节点计算函数值.

2. 三角格式

令 $T_{0i}=T\left(h_i\right)\left(i=0,1,2,\cdots \right)$ ,进行递归计算

$$\begin{array}{}\text{(19.87)}& T_{ki}=T_{k-1,i}+\frac{T_{k-1,i}-T_{k-1,i-1}}{4^k-1}\left(k=1,2,\cdots ,m;i=k,k+1,\cdots \right).\end{array}$$

根据 (19.87) 计算得到的数值经常排列成三角格式, 其元素按逐列方式计算:

$$T\left(h_0\right)=T_{00}$$$$T\left(h_1\right)=T_{01}{T}_{11}$$$$\begin{array}{}\text{(19.88)}& T\left(h_2\right)=T_{02}{T}_{12}{T}_{22}\end{array}$$$$T\left(h_3\right)=T_{03}{T}_{13}{T}_{23}{T}_{33}$$

... ...

该格式会一直持续下去 (对于给定的列数), 直到在右边最下面的数值几乎相同. 第二列的值 $T_{1i}\left(i=1,2,\cdots \right)$ 相应于由辛普森公式计算得到的值.

19.3.4.2 外推原理

龙贝格方法应用了所谓外推原理. 当 $k=1$ 时通过推导 (19.86) 证明之. 需要计算的积分表示为 $I$ ,相应的梯形和 (19.76) 表示为 $T\left(h\right)$ . 若 $I$ 的被积函数在积分区间是 $2m+2$ 次连续可微的,则能证明关于 $h$ 的渐进展开对求积公式的误差 $R$ 是成立的, 且有形式

$$\begin{array}{}\text{(19.89a)}& R\left(h\right)=I-T\left(h\right)=a_1{h}^2+a_2{h}^4+\cdots +a_m{h}^{2m}+O\left(h^{2m+2}\right)\end{array}$$

或

$$\begin{array}{}\text{(19.89b)}& T\left(h\right)=I-a_1{h}^2-a_2{h}^4-\cdots -a_m{h}^{2m}+O\left(h^{2m+2}\right),\end{array}$$

系数 $a_1,a_2,\cdots ,a_m$ 为与 $h$ 无关的常数.

根据(19.89b),考虑 $T\left(h\right)$ 与 $T\left(\frac{h}{2}\right)$ 及其线性组合

$$T_1\left(h\right)={\alpha }_{1}T\left(h\right)+{\alpha }_{2}T\left(\frac{h}{2}\right)=\left({\alpha }_1+{\alpha }_2\right)I-a_1\left({\alpha }_1+\frac{{\alpha }_2}{4}\right)h^2$$$$\begin{array}{}\text{(19.90)}& -a_2\left({\alpha }_1+\frac{{\alpha }_2}{16}\right)h^4-\cdots \end{array}$$

若以 ${\alpha }_1+{\alpha }_2=1$ 及 ${\alpha }_1+\frac{{\alpha }_2}{4}=0$ 代入,则 $T_1\left(h\right)$ 有 4 阶误差,而 $T\left(h\right)$ 与 $T\left(\frac{h}{2}\right)$ 仅有 2 阶误差. 公式为

$$\begin{array}{}\text{(19.91)}& T_1\left(h\right)=-\frac{1}{3}T\left(h\right)+\frac{4}{3}T\left(\frac{h}{2}\right)=T\left(\frac{h}{2}\right)+\frac{T\left(\frac{h}{2}\right)-T\left(h\right)}{3}.\end{array}$$

这是 $k=1$ 时的公式 (19.87),根据 (19.87) 重复应用上述过程就得到近似 $T_{ik}$ ,且

$$\begin{array}{}\text{(19.92)}& T_{ki}=I+O\left(h_i^{2k+2}\right).\end{array}$$

$◼$ 定积分 $I={\int }_0^1\frac{\sin x}{x}\mathrm{d}x$ (参见第 681 页 8.2.5,1.) 不能用基本方法得到. 计算该积分的近似值 (保留 8 位数字).

1. 龙贝格方法

$\mathrm{k}=0$	$\mathrm{k}=1$	$\mathrm{k}=2$	$\mathrm{k}=3$
0.92073549
0.93979328	0.94614588
0.94451352	0.94608693	0.94608300
0.94569086	0.94608331	0.94608307	0.94608307

由龙贝格方法得到近似值 0.94608307. 按 10 位数字得到的是 0.9460830704. 根据(19.92)验证了误差阶为 $O\left({\left(1/8\right)}^8\right)\approx 6\cdot {10}^{-8}$ .

2. 梯形与辛普森公式

根据龙贝格方法的格式,对于 $h_3=1/8$ 能直接得到梯形公式的近似值 0.94569086, 而辛普森公式得到的近似值为 0.94608331 .

根据 (19.79) 修正埃尔米特梯形公式得到的结果为

$$\$◼\$\approx 0.94569086+\frac{0.30116868}{64\cdot 12}=0.94608301.$$

3. 高斯公式

由公式 (19.83) 得

$$n=1:I\approx \frac{1}{2}\left[c_{0}f\left(\frac{1}{2}{x}_0+\frac{1}{2}\right)+c_{1}f\left(\frac{1}{2}{x}_1+\frac{1}{2}\right)\right]$$$$=0.94604113$$$$n=2:I\approx \frac{1}{2}\left[c_{0}f\left(\frac{1}{2}{x}_0+\frac{1}{2}\right)+c_{1}f\left(\frac{1}{2}{x}_1+\frac{1}{2}\right)+c_{2}f\left(\frac{1}{2}{x}_2+\frac{1}{2}\right)\right]$$$$=0.94608313$$$$n=3:I\approx \frac{1}{2}\left[c_{0}f\left(\frac{1}{2}{x}_0+\frac{1}{2}\right)+\cdots +c_{3}f\left(\frac{1}{2}{x}_3+\frac{1}{2}\right)\right]$$$$=0.94608307\text{.}$$

可见对于 $n=3$ 的高斯积分公式,仅用四个函数值就可得到 8 位数字的准确值. 而用梯形公式要达到这个精度需要非常多 (大于 1000) 个函数值.

注 (1) 傅里叶分析在周期函数的积分中起重要作用 (参见第 633 页 7.4.1.1,1.). 其数值实现的细节可在调和分析的目录下找到 (参见第 1287 页 19.6.4). 实际计算基于所谓快速傅里叶变换 FFT (参见第 1288 页 19.6.4.2).

(2) 在许多应用中, 考虑积分的特殊性质是有用的. 对这些特殊情况发展了进一步的积分法则. 关于收敛性、误差分析、最优积分公式的大量讨论可见文献 (例如见 $\left[19,7\right])$ .

(3) 文献中讨论了求多重积分值的数值方法 (例如见 [19.34]).