第四章 数值积分与数值微分 4.3 Romberg算法§ 4.3 Romberg算 法§ 综合前几节的内容 ,我们知道 梯形公式 ,Simpson公式 ,Cotes公式的代数精度分别为 1次 ,3次和 5次 复合梯形 、 复合 Simpson、 复合 Cotes公式的收敛阶分别为 2阶 、 4阶和 6阶 无论从代数精度还是收敛速度 ,复合梯形公式都是较差的 有没有办法改善梯形公式呢 ? 一 、 复合梯形公式的递推化 等份分割为的积分区间将定积分 nbadxxfI b a ],[)(∫= njjhaxk ,,1,0, L=+= n abh ?=各节点为 ])()(2)([2 1 1 ∑? = ++?= n j jn bfxfafn abT 复合梯形 (Trapz)公式为 则不变等份 ， 而分割为如果将 ,/)(2],[ nabhnba ?= )]()(2)(2)([4 1 0 2 1 1 1 2 bfxfxfafn abT n j j n j jn +++ ?= ∑∑ ? = + ? = --------(1) --------(2) )]()(2)(2)([4 1 0 2 1 1 1 2 bfxfxfafn abT n j j n j jn +++ ?= ∑∑ ? = + ? = hjahxx j j )21(21 2 1 ++=+=+其中 ∑∑ ? = + ? = ?+++?= 1 0 2 1 1 1 )(24)]()(2)([4 n j j n j j xfn abbfxfaf n ab ∑? = + ?+= 1 0 2 1 )(22 1 n j j n xfn abT ∑? = ++?+= 1 0 ))21((221 n j n hjafn abT --------(3)∑ ? = ?++?+= 1 0 )2)12((221 n j n n abjaf n abT abhn ?== 时 ，1 则由 (1)(2)(3)式 ,有 )]()([21 bfafabT +?= )21(221 12 hafabTT +?+= )0(0T= )1(0T= )1(0 ?= kTTn记12 ?= kn若 L,2,1=k jhaxj +=1 2 ? ?= k abh hxx j j 2 1 2 1 +=+ 12)2 1( ? ?++= k abja 12 ? ?+= k abja k abja 2)12( ?++= 因此 (1)(2)(3)式可化为如下递推公式 )]()([2 bfafab +?=)0(0T ∑? = ? ? ++?+?= 12 0 00 1 )2)12((2)1(21)( k j kk abjafabkTkT L,2,1=k ?? ???(4)------- 上式称为 递推的梯形公式 递推梯形公式加上一个控制精度 ,即 可成为自动选取步长的复合梯形公式 具体的方法请同学们完成 思考 二 、 外推加速公式 由复合梯形公式的余项公式 )(31 22 nnn TTTI ?≈? nn TTI 3 1 3 4 2 ?≈可得 n n j j n Txfn abTI 3 1))( 22 1( 3 4 1 0 2 1 ? ?+≈ ∑? = + 由 (3)式 ∑? = + ?+≈ 1 0 2 1 )(6 )(4 3 1 n j j n xfn abT 12 ?= kn设 )1( 3 1)( 3 4 00 ??= kTkT ∑∑ ? = + ? = ?+++?≈ 1 0 2 1 1 1 )(6 )(4])(2)()((2[31 n j j n j j xfn abxfbfaf n abI ])(4)(2)()((6 1 0 2 1 1 1 ∑∑ ? = + ? = +++?= n j j n j j xfxfbfafn ab nS= 复合 Simpson公式 nn TTI 3 1 3 4 2 ?≈ 令引入 ),1(1 ?kT )1(1 ?kT )1(31)(34 00 ??= kTkT nS= 12 ?= kS --------(5) --------(6) )(151 22 nnn SSSI ?≈? 因此由复合 Simpson公式的余项 可得 nn SSI 1511516 2 ?≈ )1(1 ?= kT12 ?= kS即 )1(151)(1516 11 ??= kTkT nS )(1 kT=nS2当然 )1(2 ?kT )1(151)(1516 11 ??= kTkT令 nC= 自己 证明 --------(6) nC= --------(7) )1(2 ?= kT12 ?= kCnC --------(8)即 )(2 kT=nC2当然 同样由复合 Cotes公式的余项 )(631 22 nnn CCCI ?≈? nn CCI 63 1 63 64 2 ?≈ )1(63 1)( 63 64 22 ??= kTkT得 )1(3 ?kT令 )1(631)(6364 22 ??= kTkT --------(9) )1(1 ?kT )1(31)(34 00 ??= kTkT )1(151)(1516 11 ??= kTkT)1(2 ?kT )1(3 ?kT )1(631)(6364 22 ??= kTkT )]()([2 bfafab +?=)0(0T ∑? = ? ? ++?+?= 12 0 00 1 )2)12((2)1(21)( k j kk abjafabkTkT L,2,1=k ? ? ? ? ? ? ? 外推 加速 公式 以上 整个过程称为 Romberg算法 将 上 述 结 果 综 合 后 )]1()(4[14 1)1( 11 ???=? ?? kTkTkT mmmmm 其中外推加速公式可简化为 --------(9) )0(0T )1(0T )0(1T )2(0T )1(1T )0(2T )3(0T )2(1T )1(2T )0(3T L,2,1=mm可以推广到并且 L,2,1=k Romberg算法的收敛 阶高达 m+1的两倍 Romberg算法求解步骤 Romberg算法的 代 数 精度为 m的两倍 romberg.m例 : ∫20 sin p xdx计算积分 前侧矩形公式 z1 = 0.99212545660563 z11 = 0.99212545660563 后侧矩形公式 z2 = 1.00783341987358 z22 = 1.00783341987358 梯形公式 z3 = 0.99997943823961 Simpson公式 z4 = 1.00000000201613 8阶 simpson公式 z5 =1.00000000000000 自选步长梯形公式 z6 = 0.99999921563419 自选步长 Simpson公式 z7 =1.00000051668471 Romberg公式 z8 = 0.99999999999802 Mote-Carlo算法 z9 = 0.99821071589516 -0.00787454339437 -0.00787454339437 0.00783341987358 0.00783341987358 -0.00002056176039 0.00000000201613 -0.00000000000000 -0.00000078436581 0.00000051668471 -0.00000000000198 -0.00178928410484 Jifenbijiao.m 积分法 积分值 绝对误差 如何构造 Romberg算法