Romberg算法

Cha4.2.m
程序演示
Cha41.m: 复化N-C公式 Cha42.m: 梯形法的递推
龙贝格 (Romberg) 加速
梯形法递推公式算法简单，编程方便
但收敛速度较 慢 解决方法：龙贝格 (Romberg) 加速
复化积分公式的渐近状态
思想：利用余项公式ห้องสมุดไป่ตู้积分的定义
梯形公式：
Cotes公式
16 1 16 1 S 2 n S n 而 Cn S 2 n S n 15 15 15 15
Cotes公式的加速
Romberg公式
64 1 64 1 I C2 n Cn 令 Rn C2 n Cn 63 63 63 63
Richardson 外推加速法
上述加速过程不是出于偶然，而且可以继续进 行下去。其理论保证是梯形法的余项展开式 定理：设 f(x)C [a, b], 记 Tn = T (h), 则有
n1
梯形法递推公式
1 h n1 1 h n1 T2 n Tn f ( xi 1 2 ) Tn f ( a ih 0.5h) 2 2 i 0 2 2 i 0
ba T1 f (a) f (b) 2 h0 0 1 T2 T1 f (a ih0 0.5h0 ) 2 2 i 0
定义
如果一种复化求积公式 I n 满足下列关系
I In lim p C , C 0 h 0 h
则称该求积公式是 p 阶收敛的. 有如下误差估计式
h2 I Tn [ f ' (b) f ' (a)] 12 h4 I Sn [ f ' ' ' (b) f ' ' ' (a)] 4 180 2
I [ f ]=0.4

1
0
x 3 dx , 要求计算精度
Cha44.m
解：逐步计算可得
k
0 1 2
T0(k)
0.50000000
T1(k)
T2(k)
T3(k)
T4(k)
T5(k)
0.42677670 0.40236893 0.40701811 0.40043192 0.40030278
3 3 3
0.40181246 0.40007725 0.40005361 0.40004965 0.40046340 0.40001371 0.40000948 0.40000878 0.40000862 0.40011767 0.40000243 0.40000168 0.40000155 0.40000152 0.40000152
梯形法的余项展开式的推导
T (h) I[ f ] 1h2 2h4 i h2i
必要性: Richardson外推法的基础
推导方法：Taylor展开方法
作业
1.证明：梯形公式的Romberg加速为辛普森公式 2. 教材 P 104/11
ba h n
h0 b a
h1 ba 2
n=1
n=2 n=4
1 h1 1 T4 T2 f (a ih1 0.5h1 ) 2 2 i 0 1 h2 3 T8 T4 f ( a ih2 0.5h2 ) 2 2 i 0
ba h2 4
梯形法递推公式
hk 1 1 T2k T2k 1 2 2
2k 1 1

i 0
f ( a ihk 1 0.5hk 1 )
hk 1
ba k 1 2
k 1
记 T
(k)
T2k
2k 1 1
n2
T
(k)
1 ( k 1) hk 1 T 2 2

i 0
f (a ihk 1 0.5hk 1 )
hk 1 ba 2 k 1
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分 解：
T (0)

1
0
sin( x ) dx , 要求 x
ba f ( a) f ( b) 0.920735492 2 1 (0) h0 0 T (1) T f (a ih0 0.5h0 ) 0.939793285 2 2 i 0 1 (1) h1 1 ( 2) T T f (a ih1 0.5h1 ) 0.944513522 2 2 i 0 1 ( 2) h2 1 ( 3) T T f ( a ih2 0.5h2 ) 0.945690864 2 2 i 0
步长折半：[xi , xi+1/2] ， [xi +1/2 , xi+1]
h T2 n f ( xi ) f ( xi 1 2 ) f ( xi 1 2 ) f ( xi 1 ) i 0 4 n1 h f ( xi ) 2 f ( xi 1 2 ) f ( xi 1 ) i 0 4 h n1 h n1 1 h n1 f ( xi ) f ( xi 1 ) f ( xi 1 2 ) Tn f ( xi 1 2 ) 4 i 0 2 i 0 2 2 i 0
1 S1 4T2 T1 3 0.946145882
1 S2 4T4 T2 3 0.946086934
1 C1 16 S2 S1 0.946083004 15
Cha43.m
举例
例：用 Romberg 算法计算定积分
满足
( ( ) | Tmk ) Tmk 1 | 107
T (h) I[ f ] 1h2 2h4 i h2i
ba h n
其中系数
k , k 1,2,3,...
与h无关
梯形法的加速
T (h) I[ f ] 1h2 2h4 i h2i I[ f ] O( h2 )
h h h h T I [ f ] 1 2 i 2 2 2 2
2 4 2i
4T ( h/2) T (h) 3I[ f ] ( 3 / 4)2h4 ( 15 / 16)3h6

Richardson 外推算法
Romberg 算法
记: T
(k) 0
T2k , T
(k) 1
S2k , T
(k) 2
C2k , T
(k) 3
R2k
T0( k ) : k 次等分后梯形公式计算所得的近似值 ( Tmk ) : m 次加速后所得的近似值
(0)
(1) (2) (3)
( Tmk ) ( ( 4 mTmk 1) Tmk )1 1 4m 1
I Tn 1 [ f ' (b) f ' (a)] 2 h 12
辛普森公式：
I Sn 1 [ f ' ' ' (b) f ' ' ' (a )] 4 4 h 180 2
Cotes公式：
I Cn 2 [ f (5) (b) f ( 5) (a)] h6 945 46
I [ f ]=0.946083070367

1
0
sin( x ) dx x
T4 0.944513522
k
0 1 2 3 4 5 6 7 8 9 10
T2 0.939793285
T0(k)
0.920735492 0.939793285 0.944513522 0.945690864 0.945985030 0.946058561 0.946076943 0.946081539 0.946082687 0.946082975 0.946083046
数值分析及计算软件
Chap 4 数值积分与数值微分
4.1 Romberg 算法
Romberg 算法
利用复化梯形公式、复化Simpson公式、复化 Cotes公式等计算定积分时，如何选取步长 h 太 大: 太 小: 计算精度难以保证
？
计算量增大
解决办法：变步长算法 通常采取将区间不断对分的方法，即取 n = 2k ， 反复采用复化求积公式，直到所得到的计算结果 满足指定的精度为止。
1 h S ( h) 4T T ( h) I [ f ] 1h4 2h6 I [ f ] O( h4 ) 3 2 1 6 h C ( h) 16 S S ( h) I [ f ] 1h6 2h8 I [ f ] O( h ) 15 2 1 h R( h) 64C C ( h) I [ f ] O( h8 ) 53 2
2h 6 I Cn [ f (5) (b) f (5) (a)] 945 46
Romberg加速
梯形公式的加速
I T2 n 1 I Tn 4
辛普森公式
辛普森公式的加速
I S2n 1 I S n 16 I C2 n 1 I Cn 64
I
4 1 4 1 而 S n T2 n Tn I T2 n Tn 3 3 3 3 T2 n Tn 事后估计法 I T2 n 3
梯形法递推公式
将 [a, b] 分成 n 等分 [xi , xi+1] ，xi a i h,
h ba n n1 n1 h f ( x ) f ( x ) h f (a) 2 f ( x ) f (b) Tn i i i 1 2 2 i 0 i 1
① T1 =T0
② T2 =T0
③ S1 =T1
(0) (1) (2)
Romberg 算法是收敛的
⑥ C1 =T2
⑨ C2 =T2
(0) (1) (0)
④ T4 =T0
⑦ T8 =T0
⑤ S2 =T1
⑧ S4 =T1
⑩ R1 =T3





举例
例：计算定积分
T1 0.920735492