龙贝格求积法

1、4.4 外推原理与外推原理与Romberg求积方法求积方法4.4.1 外推原理外推原理 在科学与工程计算中，很多算法与步长h有关，特别是数值积分、数值微分和微分方程数值解的问题。对于这些算法，我们可以通过外推技巧提高计算精度。例1 计算的近似值。sin x由函数 的Taylor展开式有3524sin3!5!nnnn,( )6sin66hF h若记 ，则有24( ),6120F hhh2411( ),26 412016hFhh由此构造新的表达式414 ( )( )12( )3120 4hFF hF hh 可见计算的近似值的算法F(h)的截断误差是O (h2),而算法F1(h)的截断误差是O (h

2、4)。外推一次，精度就提高了。这就是外推法的基本思想。若重复以上过程，不断外推，即不断折半步长h, 得到计算的算法序列Fk(h)。随着k的增加，算法的截断误差阶越来越高，计算精度越来越好。 可将上述外推思想推广到一般情况。设F(h)是计算F(0)的一种近似算式，带截断误差的表达式为( )(0)(),pspF hFa hO hsppa( )(0)( )(),psphhFFaO hqq1( )( )( )(0)(),1psphq FF hqF hFO hq().sO h其中， 与h无关。如果我们用h和h/q (q1)两种步长分别计算F(h)和F(h /q)，则有削去截断误差的主项，得新的算法我们称

3、这个计算过程为Richardson外推法。 F1(h)逼近F(0)的截断误差是 只要知道F(h)的更加完整的关于h幂的展开式，而无需知道展开式中各个系数的具体数值，就能重复使用Richardson外推法，直到截断误差达到容许误差。用归纳法可以证明下面更一般的定理。定理4.5 假设F(h)逼近F(0)的余项为312123( )(0),pppF hFa ha ha h其中，123,(0,1,2,)kpppa k是与h无关的非零常数。则由01()( )( )( ),( ),0,1,2,1kkpkkkphq FF hpF hF h Fhkq定义的序列Fk(h)有123( )( )( )123( )(0

4、),nnnpppnnnknnnF hFahahah( )(1,2,3,)nn kakh其中与无关，q1. Richardson外推法应用非常广泛且有效，下面介绍应用于数值积分的情形。 变步长梯形求积法变步长梯形求积法算法简单，但精度较差，收敛速度较慢，算法简单，但精度较差，收敛速度较慢，但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙贝格求积公式。贝格求积公式。龙贝格公式龙贝格公式又称又称逐次分半加速法逐次分半加速法。 Romberg算法算法是在积分区间逐次分半的过程中，对用复化梯是在积分区间逐次分半的过程中，对用复化梯形法形法产

5、生的近似值产生的近似值进行进行加权平均加权平均，以获得准确度较高的近似值的，以获得准确度较高的近似值的一种方法，具有公式简练，使用方便，结果较可靠的优点。一种方法，具有公式简练，使用方便，结果较可靠的优点。二、二、Romberg算法算法 根据梯形法的误差公式，可知积分值根据梯形法的误差公式，可知积分值 的截断误差大致与的截断误差大致与 成正比，因此当步长二分后，截断误差将减至原有误差的成正比，因此当步长二分后，截断误差将减至原有误差的1/4，即有即有 将上式移项整理，可得将上式移项整理，可得nT2h（4.4.3） 由此可见，只要二分前后的两个积分值由此可见，只要二分前后的两个积分值 与与 相当

6、接近，就可以保相当接近，就可以保证证 计算结果的误差很小。这样直接用计算结果来估计误差的方法计算结果的误差很小。这样直接用计算结果来估计误差的方法通常称作通常称作误差的事后估计法误差的事后估计法。 nT2nT2nT按式按式(4.4.3),积分近似值与积分近似值与 的误差大致等于的误差大致等于 因此如果因此如果用这个误差值作为用这个误差值作为 的的一种补偿一种补偿，可以期望所得到的，可以期望所得到的 213nnTT2nT2nT（4.4.4）可能是更好的结果。可能是更好的结果。nnnTTT31342(4.4.5) 有可能比有可能比 更好地接近于积分更好地接近于积分 的真值的真值 InT2 dxxf

7、ba)(3122nnnnTTTT即即这就是说，用梯形法二分前后两个积分值这就是说，用梯形法二分前后两个积分值 和和 作线性组合作线性组合nTnT2与与这表明在收敛缓慢的梯形数列这表明在收敛缓慢的梯形数列 的基础上，若将的基础上，若将 kT2nT2nT按按(4.4.6)作线性组合就可产生收敛速度较快的作线性组合就可产生收敛速度较快的Simpson序列：序列：kS21S：、2S、关于关于(4.4.3.6)的证明：由的证明：由(4.4.1)可知：可知：这种新近似值这种新近似值 实质上又是什么呢？可以验证：实质上又是什么呢？可以验证：nTnnST 即：即：144313422nnnnnTTTTS (4.

8、4.6)4.S31211223412nknnnnabkafnabTTT 故（故（4.4.6）式成立）式成立. 1111122322nnkkkhfafxf bhfakh 111102246nnknkkkhfafxfxf bS(由上节梯形公式由上节梯形公式)（由上节（由上节Simpson公式）公式）nnnSSSI22151144151151622222nnnnnSSSSC同理，由上节近似式同理，由上节近似式类似推导可得：类似推导可得：（4.4.7）即将即将Simpson序列序列 按按(4.4.7) 作线性组合就可产生收敛作线性组合就可产生收敛kS2速度更快的新序列速度更快的新序列1242:,kCC

9、 C C -Cotes序列序列即在即在Cotes序列序列 的基础上，产生了一个称为的基础上，产生了一个称为RombergkC2 :2kR1R、2R、4R.序列序列的新序列的新序列（4.4.8）14463163643232nnnnnCCCCR 由近似式 ， 类似推导可得：)(63122nnnCCCI 我们在变步长的过程中运用公式我们在变步长的过程中运用公式(4.4.4-4.4.8)(也称它们为也称它们为加速公式加速公式) ，就能将粗糙的梯形值，就能将粗糙的梯形值Tn逐步加工成精度较高的辛逐步加工成精度较高的辛普森值普森值Sn 、柯特斯值、柯特斯值Cn和龙贝格值和龙贝格值Rn .可以证明：可以证明

10、：当当f（x）满足一定条件时，）满足一定条件时，Romberg序列序列 kR2比比Cotes序列序列 能更快地收敛到积分能更快地收敛到积分 的真值的真值I。kC2 dxxfba越精确的近似值越精确的近似值 也就是将收敛速度缓慢也就是将收敛速度缓慢综上可知：在积分区间逐次分半过程中利用公式综上可知：在积分区间逐次分半过程中利用公式可以将粗糙的近似值可以将粗糙的近似值 逐步地逐步地“加工加工”成越来成越来的梯形序列的梯形序列 逐步地逐步地“加工加工”成收敛速度越来越快的新成收敛速度越来越快的新 kT2新序列新序列 4,3,2,nnnRCS,222kkkRCS(4.4.3-4.4.8) 这种加速的方

11、法称为这种加速的方法称为 Romberg算法。其算法。其“加工加工”过程如下图，过程如下图，其中圆圈中号码表示计算顺序。其中圆圈中号码表示计算顺序。 1T2T4T8T16T1S2S4S8S1C2C4C1R2R22241161641,3315156363 1,2,4,8,.nnnnnnnnnSTT CSSRCCnT1,T2、T4、T6 由梯由梯形公式、复化梯形公式形公式、复化梯形公式的递推化公式求得的递推化公式求得kk2kT212 kS22 kC32 kR0 20=1 T11 21=2 T2 S12 22=4 T4 S2 C13 23=8 T8 S4 C2 R14 24=16 T16 S8 C4

12、 R25 25=32 T32 S16 C8 R4 区间等分数区间等分数 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龙贝格序列龙贝格序列 龙贝格求积算法也可用下表来表示：龙贝格求积算法也可用下表来表示： 例例1 1: 利用利用Romberg算法计算算法计算.,8421TTTT214)(, 1, 0 xxfba3)24(21) 1 ()0(211ffT1 . 351621321)(21212112fTT解解: :由题意由题意12041dxx计算到计算到R1.121413.13333333STT3142441111( )( )3.1(3.7647062.56)3.13117724

13、24TTff242413.14156933STT357184888811( )( )( )( )3.13898928TTffff1211613.1421181515CSS484413.14159333STT121413.13333333STT242413.14156933STT2421613SS1216413.1415862926363RCC 例例2 用用Romberg算法计算算法计算 得到的梯形值，计得到的梯形值，计算结果见算结果见 表表4-5(k代表二分次数代表二分次数)。计算值的误差不超过。计算值的误差不超过0.5 10-6.表表4-5我们看到，这里利用二分我们

14、看到，这里利用二分3次的数据次的数据(它们的精度都很差，只有二三位它们的精度都很差，只有二三位是有效数字是有效数字)，通过三次加速求得，通过三次加速求得 =0.9460831，这个结果的每一位，这个结果的每一位数字都是有效数字，可见加速的效果是十分显著的。三次外推后达到数字都是有效数字，可见加速的效果是十分显著的。三次外推后达到6位有效数字。位有效数字。10sin.xIdxx注意注意 教材中介绍的教材中介绍的Richardson外推法外推法，为便于上机计算，引用记，为便于上机计算，引用记号号 来表示各近似值，其中来表示各近似值，其中k仍代表积分区间的二分次数，而仍代表积分区间的二分次数，而下标

15、下标m则指出了近似值则指出了近似值 所在序列的性质。如所在序列的性质。如m=1在梯形序列中，在梯形序列中，m=2在在Simpson序列中，序列中，m=3在在Cotes序列中，序列中， 引入上面记号引入上面记号后，后，Romberg算法所用到的各个计算公式可统一化为：算法所用到的各个计算公式可统一化为： kmT kmT三、三、Romberg 算法计算公式的简化算法计算公式的简化(0)1( )( )2baTf af b211111121() ,1,2,.222iiiiijbajTTf abai( )(1)(1)14,1,2,.,1,2,.,41mkkkmmmmTTTmki由此可逐行构造由此可逐行构

16、造 出一个三角形数表出一个三角形数表-称为称为T数表数表 ki0123Romberg算法中止准则，一般取同列或同行相邻两数值的误差算法中止准则，一般取同列或同行相邻两数值的误差绝对绝对实际计算中常常只计算到第实际计算中常常只计算到第4列列,只使用只使用3次理查逊外推法。次理查逊外推法。 注：注：值小于事先给定的精度要求。取最后一次的数值积分值作为积值小于事先给定的精度要求。取最后一次的数值积分值作为积分的分的近似值。近似值。 外推次数外推次数分分半半次次数数复化梯复化梯形序列形序列Simpson序列序列Cotes序列序列Romberg序列序列 T1 = T8 = T4 = T2 = ? ? ?

17、 S1 = R1 = S2 = C1 = C2 = S4 =1kT234 kkkTTT01112131TTTT021222TTT0313TT04T1、在上面、在上面“加工加工”过程中的系数过程中的系数 和和 ，当，当 m 4 时，时， 而另一个而另一个 系数则接近于系数则接近于1，也就是，也就是144mm141m2551141m新公式与原公式差别不大，新公式与原公式差别不大， 但工作量却大增。但工作量却大增。 kmkmTT1因此，在实际计算中常规定因此，在实际计算中常规定 m 3，即计算到出现，即计算到出现Romberg序列为止。序列为止。 2、 可用二维数组来存放并参加运算，也可用一维数组。

18、可用二维数组来存放并参加运算，也可用一维数组。 kmT四、几点说明：四、几点说明：3、 对于积分限为无穷的积分对于积分限为无穷的积分 ，可利用变量代换化成有限区，可利用变量代换化成有限区 间的积分然后再进行计算。例如：间的积分然后再进行计算。例如：4、若被积函数有奇异点（间断点）存在于积分区间内，则、若被积函数有奇异点（间断点）存在于积分区间内，则 可将积分可将积分 区间分成小部分，使间断点在子区间的端点处。区间分成小部分，使间断点在子区间的端点处。 也可用变量代换法处理。也可用变量代换法处理。 1321xxdxtx1令令则则代入得：代入得：1030132213211111tttdttttdxxxdx21dxdtt 例例3 用龙贝格方法计算椭圆用龙贝格方法计算椭圆 x2/4 + y2 l 的周长，使结果的周长，使结果具有五位有效数字具有五位有效数字 分析分析 为便于计算，先将椭圆方程采用参数形式表示为便于计算，先将椭圆方程采用参数形式表示, ,再根再根据弧长公式将椭圆周长用积分形式表示由于计算结果要求具据弧长公式将椭圆周长用积分形式表示由于计算结果要求具有五位有效数字，因此需要估计所求积分值有几位整数，从而有五位有效数字，因此需要估计所求积分值有几位整数，从而确定所求积分值的绝对误差限最后再应