《现代数值计算》课件4.4 Romberg求积法

1、4.4 外推原理与Romberg求积方法4.4.1 外推原理 在科学与工程计算中，很多算法与步长h有关，特别是数值积分、数值微分和微分方程数值解的问题。对于这些算法，我们可以通过外推技巧提高计算精度。例1 计算的近似值。由函数 的Taylor展开式有若记 ，则有由此构造新的表达式 可见计算的近似值的算法F(h)的截断误差是O (h2),而算法F1(h)的截断误差是O (h4)。外推一次，精度就提高了。这就是外推法的基本思想。若重复以上过程，不断外推，即不断折半步长h, 得到计算的算法序列Fk(h)。随着k的增加，算法的截断误差阶越来越高，计算精度越来越好。 可将上述外推思想推广到一般情况。设F

2、(h)是计算F(0)的一种近似算式，带截断误差的表达式为其中， 与h无关。如果我们用h和h/q (q1)两种步长分别计算F(h)和F(h /q)，则有削去截断误差的主项，得新的算法我们称这个计算过程为Richardson外推法。 F1(h)逼近F(0)的截断误差是 只要知道F(h)的更加完整的关于h幂的展开式，而无需知道展开式中各个系数的具体数值，就能重复使用Richardson外推法，直到截断误差达到容许误差。用归纳法可以证明下面更一般的定理。定理4.5 假设F(h)逼近F(0)的余项为其中，是与h无关的非零常数。则由定义的序列Fk(h)有无关，q1. Richardson外推法应用非常广泛

3、且有效，下面介绍应用于数值积分的情形。 变步长梯形求积法算法简单，但精度较差，收敛速度较慢，但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格求积公式。龙贝格公式又称逐次分半加速法。 Romberg算法是在积分区间逐次分半的过程中，对用复化梯形法产生的近似值进行加权平均，以获得准确度较高的近似值的一种方法，具有公式简练，使用方便，结果较可靠的优点。二、Romberg算法 根据梯形法的误差公式，可知积分值 的截断误差大致与 成正比，因此当步长二分后，截断误差将减至原有误差的1/4，即有 将上式移项整理，可得（4.4.3） 由此可见，只要二分前后的两个积分值 与 相当接近，就可以保证 计

4、算结果的误差很小。这样直接用计算结果来估计误差的方法通常称作误差的事后估计法。 按式(4.4.3),积分近似值与 的误差大致等于 因此如果用这个误差值作为 的一种补偿，可以期望所得到的 （4.4.4）可能是更好的结果。(4.4.5) 有可能比 更好地接近于积分 的真值 I即这就是说，用梯形法二分前后两个积分值 和 作线性组合与这表明在收敛缓慢的梯形数列 的基础上，若将按(4.4.6)作线性组合就可产生收敛速度较快的Simpson序列：、关于(4.4.3.6)的证明：由(4.4.1)可知：这种新近似值 实质上又是什么呢？可以验证：即： (4.4.6) 故（4.4.6）式成立.(由上节梯形公式)（

5、由上节Simpson公式）同理，由上节近似式类似推导可得：（4.4.7）即将Simpson序列 按(4.4.7) 作线性组合就可产生收敛速度更快的新序列-Cotes序列即在Cotes序列 的基础上，产生了一个称为Romberg、.序列的新序列（4.4.8） 由近似式 ， 类似推导可得： 我们在变步长的过程中运用公式(4.4.4-4.4.8)(也称它们为加速公式) ，就能将粗糙的梯形值Tn逐步加工成精度较高的辛普森值Sn 、柯特斯值Cn和龙贝格值Rn .可以证明：当f（x）满足一定条件时，Romberg序列 比Cotes序列 能更快地收敛到积分 的真值I。越精确的近似值 也就是将收敛速度缓慢综上

6、可知：在积分区间逐次分半过程中利用公式可以将粗糙的近似值 逐步地“加工”成越来的梯形序列 逐步地“加工”成收敛速度越来越快的新新序列(4.4.3-4.4.8) 这种加速的方法称为Romberg算法。其“加工”过程如下图，其中圆圈中号码表示计算顺序。 T1,T2、T4、T6 由梯形公式、复化梯形公式的递推化公式求得0 20=1 T11 21=2 T2 S12 22=4 T4 S2 C13 23=8 T8 S4 C2 R14 24=16 T16 S8 C4 R25 25=32 T32 S16 C8 R4 区间等分数 梯形序列 辛普森序列 柯特斯序列 龙贝格序列 龙贝格求积算法也可用下表来表示： 例

7、1: 利用Romberg算法计算解:由题意计算到R1. 例2 用Romberg算法计算 得到的梯形值，计算结果见 表4-5(k代表二分次数)。计算值的误差不超过0.510-6.表4-5我们看到，这里利用二分3次的数据(它们的精度都很差，只有二三位是有效数字)，通过三次加速求得 =0.9460831，这个结果的每一位数字都是有效数字，可见加速的效果是十分显著的。三次外推后达到6位有效数字。注意 教材中介绍的Richardson外推法，为便于上机计算，引用记号 来表示各近似值，其中k仍代表积分区间的二分次数，而下标m则指出了近似值 所在序列的性质。如m=1在梯形序列中，m=2在Simpson序列中

8、，m=3在Cotes序列中， 引入上面记号后，Romberg算法所用到的各个计算公式可统一化为：三、Romberg 算法计算公式的简化由此可逐行构造 出一个三角形数表-称为T数表 ki0123Romberg算法中止准则，一般取同列或同行相邻两数值的误差绝对实际计算中常常只计算到第4列,只使用3次理查逊外推法。 注：值小于事先给定的精度要求。取最后一次的数值积分值作为积分的近似值。 外推次数分半次数复化梯形序列Simpson序列Cotes序列Romberg序列 T1 = T8 = T4 = T2 = ? ? ? S1 = R1 = S2 = C1 = C2 = S4 =1、在上面“加工”过程中的

9、系数 和 ，当 m 4 时， 而另一个 系数则接近于1，也就是新公式与原公式差别不大， 但工作量却大增。因此，在实际计算中常规定 m 3，即计算到出现Romberg序列为止。 2、 可用二维数组来存放并参加运算，也可用一维数组。四、几点说明：3、 对于积分限为无穷的积分 ，可利用变量代换化成有限区 间的积分然后再进行计算。例如：4、若被积函数有奇异点（间断点）存在于积分区间内，则 可将积分 区间分成小部分，使间断点在子区间的端点处。 也可用变量代换法处理。 令则代入得： 例3 用龙贝格方法计算椭圆 x2/4 + y2 l 的周长，使结果具有五位有效数字 分析 为便于计算，先将椭圆方程采用参数形式表示,再根据弧长公式将椭圆周长用积分形式表示由于计算结果要求具有五位有效数字，因此需要估计所求积分值有几位整数，从而确定所求积分值的绝对误差限最后再应用龙贝格方法计算积分 解 令 x 2cosq，y sinq ， 则椭圆的周长为下表给出了用龙贝格方法计算积分I= 1+3sin2 dx 的过程./200 1 2.356 1941 2 2.419 921 2.441 1632 4 2.422 103 2.422 830 2.421 608 3 8 2.422 112 2.422 115 2.422 067 2.422 074 4