龙贝格求积法
(由上节梯形公式)
h 6
f
n1
a 2
k 1
f
xk
n1
4
k 0
f
x
k
1 2
f
b
Sn
（由上节Simpson公式）
故（4.4.6）式成立.
同理，由上节近似式
I
S2n
1 15
S
2n
Sn
类似推导可得：
C2n
h p
)
Fk
q pk 1
(h)
,k
0,1, 2,L
定义的序列{Fk(h)}有
Fk (h) F (0) an(n)1h pn1
a h (n) pn2 n2
a h (n) pn3 n3
L
,
其中an(n)k (k 1, 2, 3,L )与h 无关，q>1.
Richardson外推法应用非常广泛且有效，下面介绍应用 于数值积分的情形。
)
p
O(hs
),
削去截断误差的主项，得新的算法
q pF( h) F(h)
F1(h)
q qp 1
F (0) O(hs ),
我们称这个计算过程为Richardson外推法。 F1(h)逼近F(0)的截 断误差是 O(hs ).
只要知道F(h)的更加完整的关于h幂的展开式，而无需知 道展开式中各个系数的具体数值，就能重复使用Richardson外 推法，直到截断误差达到容许误差。用归纳法可以证明下面更 一般的定理。
可将上述外推思想推广到一般情况。设F(h)是计算F(0)的一 种近似算式，带截断误差的表达式为
F(h) F(0) aphp O(hs ), s p,
其中，a p 与h无关。如果我们用h和h/q (q>1)两种步长分别计算
F(h)和F(h /q)，则有
F
(
h q
)
F
(0)
ap
(
h q
① T1 ② T2 ④ T4
③ S1 ⑤ S2
⑥ C1
⑨ C2
⑩ R1
⑦ T8
⑾ T16
⑧S4 ⑿ S8
⒀ C4
⒁ R2
Sn
43T2n
1 3 Tn , Cn
16 15
S2n
1 15
Sn , Rn
64 63 C2n
1 63 Cn ,
T1,T2、T4、T6 ……由梯 形公式、复化梯形公式
dx
1
t 2 1
二、Romberg算法
变步长梯形求积法算法简单，但精度较差，收敛速度较慢， 但可以利用梯形法算法简单的优点，形成一个新算法，这就是龙 贝格求积公式。龙贝格公式又称逐次分半加速法。
Romberg算法是在积分区间逐次分半的过程中，对用复化梯 形法产生的近似值进行加权平均，以获得准确度较高的近似值的 一种方法，具有公式简练，使用方便，结果较可靠的优点。
§4.4 外推原理与Romberg求积方法
4.4.1 外推原理
在科学与工程计算中，很多算法与步长h有关，特别是数值 积分、数值微分和微分方程数值解的问题。对于这些算法，我 们可以通过外推技巧提高计算精度。
例1 计算的近似值。
由函数 sin x 的Taylor展开式有
n sin
n
3
3!n2
Tn Sn 即：
Sn
4 3
T2
n
1 3
Tn
4T2n Tn 4 1
(4.4.6)
这表明在收敛缓慢的梯形数列 T2k 的基础上，若将T2n 与 Tn
按(4.4.6)作线性组合就可产生收敛速度较快的Simpson序列：
S2k ：S1 、S2 、S4...
关于(4.4.3.6)的证明：由(4.4.1)可知：
.
可能是更好的结果。
（4.4.4）
1 Tn T2n 3 (T2n Tn )
这就是说，用梯形法二分前后两个积分值 Tn 和T2n作线性组合
即
41 Tn 3 T2n 3 Tn
有可能比
T2n
更好地接近于积分
b
a
f
xdx
的真值
I
(4.4.5)
这种新近似值 Tn实质上又是什么呢？可以验证：
T2
4
1
1 3
T1
3.133333
C1
16 15
S2
1 15
S1
3.142118
S2 3 T4 3 T2 3.141569
T8
1 2 T4
1 8
f
(
1 8
)
f
(
3 8
)
f
(
5 8
)
f
(
7 8
)
3.138989
S4
4 3
T8
1 3
T4
3.141593
C2
C k2 2
R k3 2
O O.9207355 1 O.9397933 2 O.9445135 3 0.9456909
O.9461459 O.9460869 O.9460833
O.9460830 O.9460831
O.9460831
我们看到，这里利用二分3次的数据(它们的精度都很差，只有二三位
是有效数字)，通过三次加速求得 R1 =0.9460831，这个结果的每一位
Romberg序列为止。
2、Tmk 可用二维数组来存放并参加运算，也可用一维数组。
3、 对于积分限为无穷的积分 ，可利用变量代换化成有限区
间的积分然后再进行计算。例如： dx
1 1 x2 x3
令
x1 t
则
dx
1 t2
dt
代入得：
dx
1 1 x2 x3
0
4T2n Tn
3
Tn
2ba n
n k 1
f
a 2k
1
1 2
b
n
a
1 3
1 h
3
2
f
a
n1
2
k 1
f
xk
f
b
n
2h
k 1
f
a
k
1 2
h
值小于事先给定的精度要求。取最后一次的数值积分值作为积分的
近似值。
四、几点说明：
4m
1
1、在上面“加工”过程中的系数4m 1 和4m 1 ，当
m
4
时， 1
4m 1
1 255
而另一个 系数则接近于1，也就是
新公式与原公式差别不大，Tmk Tmk1 但工作量却大增。
因此，在实际计算中常规定 m 3，即计算到出现
16 15
S2n
1 15
Sn
42
S2n Sn 42 1
（4.4.7）
即将Simpson序列 S2k 按(4.4.7) 作线性组合就可产生收敛
速度更快的新序列 C2k : C1,C2,C4,K
------------Cotes序列
由近似式
I
C2n
1 63
(C2
n
Cn )
16 15
S4
1 15
S2
3.1415946
64 1 R1 63 C2 63 C1 3.141586292
1 sin x
算例结2 果用见R表om4-b5(ekr代g算表法二计分算次I数)。0 计x算值dx.的得误到差的不梯超形过值，计
0.510-6.
表4-5
k
T2k
S k 1 2
m=2在Simpson序列中，m=3在Cotes序列中，Tmk 引入上面记号
后，Romberg算法所用到的各个计算公式可统一化为：
T (0) 1
ba
2
f
(a)
f
(b)
T1i
1 2
T1i
1
b 2i
a
2i1 j1
f
(a
2
j 1 2i (b
a) , i
证 T2n 计算结果的误差很小。这样直接用计算结果来估计误差的方法
通常称作误差的事后估计法。
按式(4.4.3),积分近似值与 T2n
的误差大致等于
1 3
T2
n
Tn
因此如果
用这个误差值作为 T2n 的一种补偿，可以期望所得到的
T
T2 n
1 3
(T2n
Tn )
4 3
T2
n
1 3
Tn
f
(
1 2
)
1 3 2
1 16 25
3.1
S1
4 3
T2
1 3
T1
3.133333
T4
1 2
T2
1 4
f
龙贝格求 积分
龙贝格(Romberg ）求积法1.算法理论Romberg 求积方法是以复化梯形公式为基础，应用Richardson 外推法导出的数值求积方法。
由复化梯形公式 )]()(2)([2222b f h a f a f h T +++=可以化为)]()]()([2[212112h a f h b f a f hT +++=＝)]([21211h a f h T ++一般地，把区间［a,b ］逐次分半k －1次，（k ＝1，2,……,n)区间长度（步长)为kk m a b h -=，其中mk ＝2k －1。
记k T ＝)1(k T由)1(k T ＝]))12(([21211)1(1∑=---++km j k k k h j a f h T 从而⎰badxx f )(＝)1(kT－)(''122k f h a b ξ- (1)按Richardson 外推思想，可将（1）看成关于k h ，误差为)(2k h O 的一个近似公式，因而，复化梯形公式的误差公式为⎰badxx f )(－)1(k T ＝......4221++kkh K h K ＝∑∞=12i i k i h K （2）取1+k h ＝k h 21有 ⎰ba dx x f )(－)1(1+k T ＝∑∞=+121221i ik ii hK （3）误差为)(2jh O 的误差公式 )(j kT＝)1(-j kT＋141)1(1)1(------j j k j k T T2。
误差及收敛性分析（1）误差，对复化梯形公式误差估计时，是估计出每个子区间上的误差,然后将n 个子区间上的误差相加作为整个积分区间上的误差。
(2）收敛性,记h x i =∆，由于∑=++=ni i i n x f x f h f T 01))]()([2)(＝))()((21101∑∑-==∆+∆n i ni i i i i x x f x x f上面两个累加式都是积分和，由于)(x f 在区间],[b a 上可积可知，只要],[b a 的分划的最大子区间的长度0→λ时,也即∞→n 时,它们的极限都等于积分值)(f I 。
龙贝格求积公式
2
f( ) 2
T3 , 1 T3 , 2
T3 , 3
T4 , 1 T5 , 1
M
T4 , 2 T5 , 2
M
T4, 3 T5 , 3
M
T4 , 4 T5 , 4
M
用龙贝格方法计算积分的步骤为：
（1）：准备初值，先用梯形公式计算积分近
似 （值 2）：：T1按 变b 2步a[长f (a梯) 形f (公b)]式计算积分近似值：
4 0.918741799 0.916327874 0.916297224 0.916294351
5 0.916905342 0.916293190 0.916290077 0.916290776
T5 , 4 0.916290776
例3：取=0.00001,用龙贝格方法计算积分
1
I
4
dx
01 x2
解：由题意
f(x)=4/(1+x2) a=0 b=1 f(0)=4 f(1)=2 由梯形公式得
T1=1/2[f(0)+f(1)]=3 计算f(1/2)=16/5 用变步长梯形公式得
T2=1/2[T1+f(1/2)]=3.1 由加速公式得
S1=1/3(4T2-T1)=3.133333333
求出f(1/4) f(3/4) 进而求得 T4=1/2{T2+1/2[f(1/4)+f(3/4)]}
Simpson加速公式： Cn
42 S2n Sn 42 1
I
C2n
1 43 1(C2n
Cn )
43C2n Cn 43 1
Cotes加速公式：Rn
43C2n Cn 43 1
Romberg 值序列
计算方法-4.6-4.7龙贝格、高斯求积公式
4 1 T1 ( k 1) T0 ( k ) T0 ( k 1) 3 3 16 1 T2 ( k 1) T1 ( k ) T1 ( k 1) 15 15 64 1 T3 ( k 1) T2 ( k ) T2 ( k 1) 63 63
k 1, 2 ,
引入 T1 (k 1),令
S n S 2 k 1
即
Sn S2 k 1 T1 ( k 1)
S 2 n T1 ( k )
--------(6)
当然
2016/8/14
10
因此由复合Simpson公式的余项
I S2 n 1 ( S2 n Sn ) 15
Cn
可得 令 即
1次,3次和5次
复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为
2阶、4阶和6阶 无论从代数精度还是收敛速度,复合梯形公式都是较差的 有没有办法改善梯形公式呢?
2016/8/14
2
一、复合梯形公式的递推化
将定积分I f ( x)dx 的积分区间 [a , b]分割为n等份 a ba x a jh , j 0 , 1 , , n 各节点为 h j n 复合梯形(Trapz)公式为
n 1 ba Tn [ f ( a) 2 f ( x j ) f (b)] 2n j 1
b
--------(1)
如果将 [a, b]分割为2n等份，而h (b a) /n不变, 则
n 1 n 1 ba T2 n [ f ( a ) 2 f ( x j ) 2 f ( x 1 ) f (b)] j 4n j 1 j 0 2 --------(2)
外推 加速 公式
龙贝格算法
龙贝格积分1. 算法原理采用复化求积公式计算时，为使截断误差不超过ε，需要估计被积函数高阶导数的最大值，从而确定把积分区间[]b a ,分成等长子区间的个数n 。
首先在整个区间[]b a ,上应用梯形公式，算出积分近似值T1；然后将[]b a ,分半，对 应用复化梯形公式算出T2；再将每个小区间分半，一般地，每次总是在前一次的基础上再将小区间分半，然后利用递推公式进行计算，直至相邻两个值之差小于允许误差为止。
实际计算中，常用ε≤-n n T T 2作为判别计算终止的条件。
若满足，则取n T f I 2][≈；否则将区间再分半进行计算，知道满足精度要求为止。
又经过推导可知，∑=-++=ni i i n n x x f h T T 112)2(221，在实际计算中，取kn 2=，则k a b h 2-=，112)1*2(2++--+=+k i i ab i a x x 。
所以，上式可以写为∑=++--+-+=+kk i k k ab i a f a b T T 211122)2)12((2211k开始计算时，取())()(21b f a f ab T +-=龙贝格算法是由递推算法得来的。
由梯形公式得出辛普森公式得出柯特斯公式最后得到龙贝格公式。
根据梯形法的误差公式，积分值n T 的截断误差大致与2h 成正比，因此步长减半后误差将减至四分之一，即有21114n n T T -≈-将上式移项整理，知2211()3n n n T T T -≈-由此可见，只要二分前后两个积分值n T 和2n T 相当接近，就可以保证计算保证结果计算结果2n T 的误差很小，这种直接用计算结果来估计误差的方法称作误差的事后估计法。
按上式，积分值2n T 的误差大致等于21()3n n T T -，如果用这个误差值作为2n T 的一种补偿，可以期望，所得的()222141333n n n n n T T T T T T =+-=-应当是更好的结果。
4.4龙贝格求积公式
4 1 T1 ( k 1) T0 ( k ) T0 ( k 1) 3 3 16 1 T2 ( k 1) T1 ( k ) T1 ( k 1) 15 15 64 1 T3 ( k 1) T2 ( k ) T2 ( k 1) 63 63
k 1, 2 ,
1 h n1 Tn f ( x 1 ) j 2 2 j 0 2 1 b a n 1 1 Tn f (a ( j )h) 2 2n j 0 2 1 b a n 1 ba Tn f ( a ( 2 j 1) ) 2 2 n j 0 2n
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
T0 (0 ) T0 (1) T0 ( 2 ) T0 ( 3) T1 (0 ) T1 (1) T1 ( 2 ) T2 (0 ) T2 (1)
T3 (0 )
龙贝格积分
例：计算
I T2n 1 I Tn 4
令
--------(7)
即 当然
Cn C2 k 1 T2 (k 1) C 2 n T2 (k )
--------(8)
同样由复合Cotes公式的余项
I C2 n 1 (C2 n Cn ) 63
64 1 64 1 得 I C2 n Cn T2 ( k ) T2 ( k 1) 63 63 63 63
43 C 2n C n Rn 3 4 1
Romberg
T1 = T0( 0 )
T2 = T0( 1 )
<?
算法： T4 = T0( 2 )
T8 = T0
(3)
S1 = T1( 0 ) S2 = T1
数值分析63 复化求积公式龙贝格求积公式讲解
只增加了一个分点
1 xk?1/ 2 ? 2 ( xk ? xk?1)
设hn=(b? a)/n, xk=a+kh n (k=0,1,? ,n),在[xk, xk+1] 上用梯形公式得
T1 ?
hn 2
?f
(
xk
)
?
f ? ( xk ? 1 )
复化求积的基本想法 :
将积分区间 [a, b]n等分, 步长
h?
b
? n
a
,
分点为
xk=a+kh (k=0,1,…,n) , 则由定积分性质知
? ? ? I ?
b
n?1
f ( x )dx ?
xk?1 f ( x )d x
a
k ? 0 xk
每个子区间 上的积分
?xk?1 f ( x )dx xk
用低阶求积公式 , 然后把所有区间的 计算结果求和 ,
注2: 同样也可用 | S4m-S2m |<ε 来控制计算的精度 . 这就是下面要介绍的 龙贝格求 积公式 .
6.4 龙贝格求积公式
6.4.1 梯形公式的递推化
复化求积方法可提高求积精度，实际计算时若
精度不够可将步长逐次分半 . 设将区间 [a, b]分为n等
分，共有 n+1个分点，如果将求积区间再分一次，则 分点增至 2n+1个，我们将二分 前后两个积分值 联系
果T8=0.9456909 只有2位有效数字，而应用复化辛普 森公式计算的结果 S4= 0.9460832 却有6位有效数字 .
注：为了利用余项公式估计误差，要求 f(x)=sin x/x 的高阶导数，由于
龙贝格求积算法
龙贝格求积算法
龙贝格求积算法（Romberg Integration Algorithm）是用于数
值积分的一种高效的迭代方法。
它通过连续的二分、四分、八分等等
区间的方式，逐渐逼近最终的积分值，从而提高计算的精度。
该算法的基本思想是利用Richardson外推技术，结合复合梯形
法则，逐渐缩小区间并增加采样点数，以得到更精确的积分值。
下面
我们来介绍龙贝格求积算法的步骤：
1. 将积分区间[a, b]进行二分，得到初始的两个子区间；
2. 对每个子区间应用复合梯形公式进行数值积分，可以得到初始的近似积分值；
3. 利用Richardson外推技术，对不同精度的积分值进行线性组合，得到更高精度的积分值；
4. 重复步骤2和3，将积分区间不断地二分，并逐步增加采样点数，直到达到所需的精度要求。
龙贝格求积算法的主要优点是在保持高精度的能够有效减少计算量。
该算法还可以通过预先计算一些常见函数在一些固定的点上的值，以进一步提高计算速度。
总结起来，龙贝格求积算法通过利用复合梯形法则和Richardson 外推技术，逐渐逼近积分值的精确结果。
它是一种高效且精确的数值
积分方法，广泛应用于科学计算和工程领域。
6b复合求积公式龙贝格算法
步长折半：[xi , xi+1/2] ， [xi +1/2 , xi+1]
n1
xi 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 13
1 I T2 n (T2 n Tn ) 3
I Tn 4( I T2n )
3I 4T2n Tn
1 3
验后误差估计式 I T2 n (T2 n Tn )
当
T2n Tn 时，T2n即为所求的近似值。
1 (T2 n Tn ) 3
是T2n 的修正项，它与T2n 之和比T2n 、 Tn更接近与真值，即它是一种补偿。
|| T T2-T|< 2-T1|<
输出T2
16
举例
计算精度满足 | T2n Tn | 107
I [ f ]=0.946083070367
例：用梯形法的递推公式计算定积分 解：
1
0
sin( x ) dx , 要求 x
k
0 1 2 3 4 5 6 7 8 9 10
T (k)
梯形法递推公式
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
计算方法--龙贝格积分法
结论:从运行结果可见， 当设定误差 eps=1e-5 时， 为了满足误差要求， 运行 n=14 时才能满足， 这样精度比较高，但是代价就是运行的时候速度有点慢。
按照题目要求， 要求尽可能的使结果精确一些， 因此我这里采用的是根据误差要求进行的计 算，计算程序如下： function [R,T,n]=Romberg3(myfunc,a,b,eps) %This function uses Romberg to calculate integration %myfunc is the integrand,a is integral lower bound %b is integral upperbound %eps is tolerance format long; n=5; while 1 T(n,:)=0; T(1,1)=(b-a)*(myfunc(a)+myfunc(b))/2; for k=1:n-1 sum=0; for i=1:2^k sum=sum+myfunc(a+(2*i-1)*(b-a)/(2^(k+1))); end T(k+1,1)=0.5*T(k,1)+(b-a)*sum/(2^(k+1)); end for i=2:n T(i,2)=T(i,1)+1/3*(T(i,1)-T(i-1,1)); end for i=3:n T(i,3)=T(i,2)+1/15*(T(i,2)-T(i-1,2)); end for i=4:n T(i,4)=T(i,3)+1/(4^3-1)*(T(i,3)-T(i-1,3)); end if abs(T(n,4)-T(n-1,4))<eps break; else n=n+1; end end R=T(n,4); 以上便是这次编写的龙贝格积分的程序。 在matlab的command window运行如下： >> myfuc=inline('log(1+x)/(1+x^2)','x'); >> eps=1e-5;
龙贝格求积法
h p
)
Fk
q pk 1
(h)
,k
0,1, 2,L
定义的序列{Fk(h)}有
Fk (h) F (0) an(n)1h pn1
a h (n) pn2 n2
a h (n) pn3 n3
L
,
其中an(n)k (k 1, 2, 3,L )与h 无关，q>1.
Richardson外推法应用非常广泛且有效，下面介绍应用 于数值积分的情形。
16 15
S4
1 15
S2
3.1415946
64 1 R1 63 C2 63 C1 3.141586292
1 sin x
算例结2 果用见R表om4-b5(ekr代g算表法二计分算次I数)。0 计x算值dx.的得误到差的不梯超形过值，计
0.510-6.
表4-5
k
T2k
S k 1 2
§4.4 外推原理与Romberg求积方法
4.4.1 外推原理
在科学与工程计算中，很多算法与步长h有关，特别是数值 积分、数值微分和微分方程数值解的问题。对于这些算法，我 们可以通过外推技巧提高计算精度。
例1 计算的近似值。