Romberg 算法

Romberg 算法是一种自适应的利用自适应的复化梯形公式进行数值积分的算法。它是 Richardson 外推法的推广。

步骤

假设我们所考虑的积分是$\int_0^1 f(x) \mathrm{d}x =: I$，不在$[0,1]$区间上的积分可以做相应的变量代换转化到该区间上去。

我们先使用两点梯形公式计算初始迭代值$T_1 = \dfrac{1}{2} (f(0) + f(1))$，将区间不断做二等分，得到一系列复化梯形公式

$$T_{2^k} = \dfrac{1}{2} T_{2^{k-1}} + \dfrac{1}{2^k} \sum_{j=1}^{2^{k-1}} f\!\left( \dfrac{2j-1}{2^k} \right), k = 1,2,\cdots.$$

这些都是后续迭代的初始值，上式给出了初始值的通式。下面我们开始第一次迭代，令

$$S_{2^k} = \dfrac{4T_{2^{k+1}}-T_{2^k}}{3}, k = 0,1,2,\cdots.$$

上式的$S_{2^k}$是 Simpson 公式的复化自适应求积公式，进一步迭代

$$C_{2^k} = \dfrac{16S_{2^{k+1}}-S_{2^k}}{15}, k = 0,1,2,\cdots.$$

上式的$C_{2^k}$是 Cotes 公式的复化自适应求积公式，进一步迭代

$$R_{2^k} = \dfrac{64C_{2^{k+1}}-C_{2^k}}{63}, k = 0,1,2,\cdots.$$

上式的$R_{2^k}$是 Romberg 公式的复化自适应求积公式。

这样一直进行下去，精度不断提高，一般意义下的递推公式为

$$T_{j,k} = \dfrac{4^j T_{j-1,k+1} - T_{j-1,k}}{4^j-1}, j,k = 0,1,2,\cdots.$$

在上面的计算中，$T_{2^k} = T_{0, k}, S_{2^k} = T_{1, k}, C_{2^k} = T_{2, k}, R_{2^k} = T_{3, k}.$ 上述计算的精度可以控制为

$$T_{j,k} - I = O(h^{2j+2}).$$

判断迭代停止的条件可以是先验的，即预先给出绝对误差$\varepsilon > 0$后判断停止的条件可以是

$$|T_{N, 2} - T_{N, 1}| < \varepsilon, \exists N \in \N^+.$$

这时可以以$T_{N, 1}$作为积分$I$的数值解。

参考资料

1. 黄云清, 《数值计算方法》, 科学出版社, 北京, 2012-06, ISBN 978-7-0302-3428-5.