龙贝格求积公式
龙贝格求 积分
龙贝格(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 。
Romberg求积公式
)
T3
(
h 2k
)
T4
(
h 2k
)
0 (1)T1 h
(3)T2 h
(6)T3 h (10)T4 h
1
(2)
T1
(
h 2
)
(5)
T2
(
h 2
)
(9)
T3
(
h 2
)
(14)
T4
(
h 2
)
2
(4)
T1
(
h 4
)
(8)
T2
(
h 4
)
(13)
T3
(
h 4
)
3
(7)
T1
(
h 8
)
(12)
T2
(
h 8
)
4
(11)
数值计算
Romberg求积公式
1.1外推法基本思想
以较小的计算量为代价,达到提高数值结 果的精度是外推法的中心思想.
设f (x) C 2k2[x h , x h ],用中心差商 22
f (x h) f (x h)
G(h)
2
2
h
逼近f ' (x)。
由Taylor展式起截断误差为
f '(x) G(h) 2h2 4h4 ... 2kh2k E(h)
4E1(h) E1(h) 4 1
,
r2 j
1 (4( j1) 3
1)2 j
( j 2,3,...k)。
此时f ' (x) G1(h) O(h4 )
同
理可
以有G1
(h)和G1
(
h 2
)构造G2
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 ,L
因此有如下递推公式 b−a [ f ( a ) + f (b )] T0 (0) = 2
1 T0 (k ) = T0 (k − 1) + hk 2
2 k −1 −1 j =0
∑ f (a + (2 j + 1)h )
k
k = 1, 2 ,L
上式称为递推的梯形公式
思考
递推梯形公式加上一个控制精度,即 可成为自动选取步长的复合梯形公式
带权）正交。 不大于n 不大于 的多项式 P(x) （带权）正交。
k =0
n
x 证明： 证明： “⇒”0 … xn 为 Gauss 点, 则公式 ∫a ρ( x) f ( x)dx ≈ ∑Ak f ( xk ) k=0 次代数精度。 至少有 2n+1 次代数精度。 对任意次数不大于 的多项式 不大于n 对任意次数不大于 ⇔ 求w(x) Pm(x)， Pm(x) w(x)的次数 ， 的次数 求 Gauss 点 不大于2n+1，则代入公式应精确成立： 精确成立： 不大于 ，则代入公式应精确成立 n 0 b ∫ ρ ( x ) Pm ( x ) w ( x )dx = ∑ Ak Pm ( x k ) w ( x k ) = 0
外推加速公式
由复合梯形公式的余项公式
I − T2 n 1 ≈ I − Tn 4 1 I − T2 n ≈ (T2 n − Tn ) 3 4 1 I ≈ T2 n − Tn 3 3 1 f ( x 1 )) − Tn j+ 3 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 =+-=-应当是更好的结果。
龙贝格求积公式
S2k
4T2k1 T2k 41
C2k
42 S2k1 S2k 42 1
R2k
4 C 3 2k 1
C2k
43 1
Romberg积分法的一般公式
Tm, j
4
T j 1 m,
j 1
Tm 1,
j 1
4 j1 1
其中
j 2,3, 4;m j
➢ Romberg积分思想
由上节分析知，用复化梯形公式计算积分值 I
T2n 的误差大约为：
1 3 (T2n Tn )
令
I
T2n
1 3 (T2n
Tn )
4T2n Tn 3
由复化梯形公式知
1
b a n1
T2n 2 Tn
2n
k0
f
(
x
k
1
)
2
4T2n Tn
3
1 3 Tn
2(b a) n1 3n k0
Tm,1 T2m1 (m 1) Tm,3 C2m3 (m 3)
Tm,2 S2m2 (m 2) Tm,4 R2m4 (m 4)
Romberg积分表
1
b a n1
T1 , 1
T2n 2 Tn
2n
k0
f
(
x
k
1
)
2
T2 , 1 T2 , 2
1 ba ba
T2,1 2 T1,1
=3.131176471 S2=1/3(4T4-T2)=3.141568628 C1=1/15(16S2-S1)=3.142117648 计算f(1/8) f(3/8) f(5/8) f(7/8)进而求得 T8=1/2{T4+1/4[f(1/8)+f(3/8)+f(5/8)+f(7/8)]}
数值分析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 的高阶导数，由于
龙贝格算法
龙贝格算法一、问题分析1、1龙贝格积分题目要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y（x)=l sin （tx），则给定雨篷得长度后，求所需要平板材料得长度）.二、方法原理2、１龙贝格积分原理龙贝格算法就是由递推算法得来得。
由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式.在变步长得过程中探讨梯形法得计算规律.设将求积区间［ａ，b]分为ｎ个等分，则一共有n+1个等分点，n.这里用表示复化梯形法求得得积分值，其下标n 表示等分数。
先考察下一个字段[]，其中点,在该子段上二分前后两个积分值显然有下列关系将这一关系式关于k从０到n－１累加求与，即可导出下列递推公式需要强调指出得就是,上式中得代表二分前得步长，而梯形法得算法简单，但精度低,收敛速度缓慢,如何提高收敛速度以节省计算量，自然式人们极为关心得.根据梯形法得误差公式，积分值得截断误差大致与成正比，因此步长减半后误差将减至四分之一,既有将上式移项整理,知由此可见，只要二分前后两个积分值与相当接近，就可以保证计算保证结果计算结果得误差很小，这种直接用计算结果来估计误差得方法称作误差得事后估计法。
ﻩ按上式,积分值得误差大致等于，如果用这个误差值作为得一种补偿，可以期望,所得得应当就是更好得结果。
ﻩ按上式，组合得到得近似值直接验证，用梯形二分前后得两个积分值与按式组合，结果得到辛普生法得积分值。
再考察辛普生法。
其截断误差与成正比.因此,若将步长折半,则误差相应得减至十六分之一。
既有由此得不难验证,上式右端得值其实就等于,就就是说,用辛普生法二分前后得两个积分值与，在按上式再做线性组合,结果得到柯特斯法得积分值，既有重复同样得手续，依据斯科特法得误差公式可进一步导出龙贝格公式应当注意龙贝格公式已经不属于牛顿—柯特斯公式得范畴.在步长二分得过程中运用公式加工三次，就能将粗糙得积分值逐步加工成精度较高得龙贝格,或者说,将收敛缓慢得梯形值序列加工成熟练迅速得龙贝格值序列，这种加速方法称龙贝格算法。
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
龙贝格积分 python
龙贝格积分 python龙贝格积分 python一、什么是龙贝格积分？龙贝格积分（Romberg integration）是一种数值积分方法，它是对梯形法的递推加速处理。
梯形法是一种比较简单的数值积分方法，但它的精度不高，而且需要很多次计算。
龙贝格积分通过递推计算，可以大大提高计算精度，并且减少计算次数。
二、如何实现龙贝格积分？在 Python 中实现龙贝格积分可以使用以下代码：```pythondef romberg(f, a, b, n):"""Calculate the Romberg Integration of f(x) from a to b with n iterations"""R = [[0] * (n+1) for _ in range(n+1)]h = b - aR[0][0] = 0.5 * h * (f(a) + f(b))for i in range(1, n+1):h = 0.5 * hsum = 0for k in range(1, 2**i, 2):sum += f(a + k*h)R[i][0] = 0.5 * R[i-1][0] + sum*hfor j in range(1, i+1):R[i][j] = (4**j*R[i][j-1] - R[i-1][j-1]) / (4**j - 1)return R[n][n]```三、代码解析1. 定义函数 romberg(f, a, b, n)，其中 f 为被积函数，a 和 b 分别为积分上下限，n 为迭代次数。
2. 创建一个n+1 行n+1 列的二维数组R，用于存储递推计算的结果。
3. 计算初始值 R[0][0]，即使用梯形法计算第一次迭代的结果。
4. 进行 n 次迭代，每次将步长 h 减半，并且计算新的递推值。
具体过程如下：a. 计算当前步长 h。
9-3-4s外推法与龙贝格求积公式
h2 = 1
⎧T1, 0 = ( h / 2 )[ f ( a ) + f ( b )], h = b − a ⎪ ⎨ T m , k +1 − T m , k , ⎪T m + 1, k = T m , k + 1 + 4m −1 ⎩
例题2 利用Romberg序列，近似计算 若T1,0=4, T2,0=5,求f (0.5).
1
3
x
cos xdx
使精度达到10-6.
引言
n+1个节点的插值求积公式
（一）定理：
n k =0
∫
b
a
f ( x) dx ≈ ∑ Ak f ( xk )
n+1个节点的插值型求积公式 ∫ ρ ( x ) f ( x ) dx ≈ a 代数精度最高不超过2n+1次。
b
∑A
k =0
n
k
f ( xk )
的
P261
Romberg方法计算
∫ f ( x)dx 步骤如下：
b a
（1）构造复合梯形序列{T1, k},k=0,1,…, 定义
T1,0= T1 =h0[f(a)+f(b)] /2
T1,k = T 2k , hk =
b−a 2k
h0=b-a
T1,1= T2 = (h1/2)[f(a)+f(b)+2f(a+h1)]=1/2[T1,0 +h0 f(a+h1)] , h1=(b-a)/2 T1,2 = T4 =(1/2){T1,1 +h1 [f(a+h2)+ f(a+3h2)] } , h2=(b-a)/4 ················
5.3龙贝格(Romberg)算法
|R1-C2|<ε
否 k=4 计算 T
、S2k -1、C2k −2、R2k −3 2
k
输出积分值 I
结束
| R2k −3 − C2k −2 |< ε
否 k=k+1
是
17
龙贝格算法源程序
#include "stdio.h" #include "math.h" main() { double f(double x); double t0fun(int k,double t0[50]);
k −3
if(fabs(t3[k-3]-t2[k-2])<0.000000001) // 若 | R2 { printf("k=%3d,I=%12.9f\n",k,t3[k-3]); break; } } } double f(double x) { double y; y=sqrt(1.0+x*x); return y; }
龙贝格romberg算法等份分割为的积分区间将定积分不变等份而分割为如果将上式称为递推的梯形公式递推梯形公式加上一个控制精度即可成为自动选取步长的复合梯形公式具体的方法请同学们完成思考因此由复合simpson公式的余项可得151615161516自己证明当然同样由复合cotes公式的余项63646364636415166364加速公式以上整个过程称为龙贝格romberg算法1012龙贝格算法可用下表表示
1 b − a n −1 1 b − a n −1 1 = Tn + ∑ f ( x j + 1 ) = 2 Tn + 2 n ∑ f ( a + ( j + 2 )h) 2 2n j =0 j =0 2
9-3-4s外推法与龙贝格求积公式
T1 , k = T 2 k =
可以证明
2 1 b−a [T1 , k + h k −1 ∑ f (a + h k ( 2 i − 1) )], h k = 2 2k i =1 k = 1, 2 ,
k −1
∑
∫
b
a
2 4 6 f ( x)dx − T1,k = a1hk + a2 hk + a3 hk +
+ ai hk2i +
1
(2)对复合梯形序列{T1,k} 作Richardson外推 Richardson外推法 若 p1 < p 2 < < pm <
Romberg积分法
T1 , k = T 2 k =
,且
p2
F * − F ( h ) = a1 h
则有
p1
+ a2 h
+
+ a m h pm +
⎧T = ( h / 2 )[ f ( a ) + f ( b )], h = b − a ⎪ 1, 0 ⎪ T m , k + 1 − (1 / 2 ) 2 m T m , k T ⎪ 1 − Tm ,k , = T m , k + 1 + m , k +m ⎨T m + 1 , k = 1 − (1 / 2 ) 2 m 4 −1 ⎪ ⎪ b 2m ⎪ ∫a f ( x ) dx = T m , k + O ( h k ), m = 1, , k , k = 0 ,1, n . ⎩ k 是 2 分区间 [ a , b ]的次数； m 是外推的次数。
其中η∈(a,b),wn+1(x)=(x-x0)(x-x1)…..(x-xn)。 （2）高斯求积公式的系数Ak(k=0,1,…,n)大于零。
龙贝格算法
龙贝格算法 一、问题分析1.1龙贝格积分题目要求学生运用龙贝格算法解决实际问题（塑料雨篷曲线满足函数y(x)=l sin(tx),则给定雨篷的长度后，求所需要平板材料的长度）。
二、方法原理2.1龙贝格积分原理龙贝格算法是由递推算法得来的。
由梯形公式得出辛普生公式得出柯特斯公式最后得到龙贝格公式。
在变步长的过程中探讨梯形法的计算规律。
设将求积区间[a ，b]分为n 个等分，则一共有n+1个等分点,k x a kh =+,0,1,b ah k n-==，n 。
这里用n T 表示复化梯形法求得的积分值，其下标n 表示等分数。
先考察下一个字段[1,k k x x +]，其中点()11212k k k xx x ++=+，在该子段上二分前后两个积分值()()112k k hT f x f x +=+⎡⎤⎣⎦ ()()21124k k k h T f x f x f x ++⎡⎤⎛⎫=++⎢⎥ ⎪⎢⎥⎝⎭⎣⎦显然有下列关系2112122k h T T f x +⎛⎫=+⎪⎝⎭将这一关系式关于k 从0到n-1累加求和，即可导出下列递推公式12102122n n n k k h T T fx -+=⎛⎫=+ ⎪⎝⎭∑ 需要强调指出的是，上式中的b ah n-=代表二分前的步长，而 1212k xa k h +⎛⎫=++ ⎪⎝⎭梯形法的算法简单，但精度低，收敛速度缓慢，如何提高收敛速度以节省计算量，自然式人们极为关心的。
根据梯形法的误差公式，积分值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 =+-=- 应当是更好的结果。
数值计算方法 龙贝格求积公式 - 龙贝格求积公式
类似地可验证：
3
S2
4
4
1
T4
4
1
1
T2
龙
贝 格 算 法
Sn
4
4
1
T2
n
4
1
1
Tn
即
Sn
4 3 T2n
1 3
Tn
Sn __ 辛 普 森 积 分 值
注 意 复 化 辛 普 森 求 积 公式 的 余 项
3
Rn
I
Sn
ba 180
h 2
4
f
(4) ( )
O(h4
)
龙 贝 格
R2 n
I
S2n
T4
1 2
T2
1 4
[
f
(
1 4
)
f
(
3 4
)]
0.9445135
典型例题
例2
用复化梯形公式计算积分I 1 dx
0 (1 x ) x
精确至三位有效数字
I
1 dx 0 (1 x )
x t x
1 2dt 0 1 t2
1
g(t )dt
0
T1
1 [g(0) 2
g(1)] 1.5
T2
1 2
[T1
1
g(0.5)] 1.55
T4
1 2
[T2
1 ( g(0.25) 2
g(0.75))] 1.5656
典型例题
T8
1 2 [T4
1 ( g(0.125) 4
g(0.375)
g(0.625)
g(0.875))]
1.5695
注意到:
T8 T4
Romberg求积分公式
《MATLAB程序设计实践》课程考核1、编程实现以下科学计算算法，并举一例应用之。
“Romberg求积分公式”2、编程解决以下科学计算和工程实际问题。
1）、给定半径的为r，重量为Q的均质圆柱，轴心的初始速度为v0，初始角速度为w0且v0>r*w0，地面的摩擦系数为f，问经过多少时间后，圆柱将无滑动地滚动，求此时的圆柱轴心的速度。
2）、在一丘陵地带测量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。
一、Romberg求积分公式1、算法说明：此算法可自动改变积分步长,使其相临两个值的绝对误差或相对误差小于预先设定的允许误差.Romberg加速法公式在等距节点的情况下，通过对求积区间（a，b）的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有精度，即，但当节点加密时可组合成其精度达到，如果再由与组合成则可使误差精度达到，于是依赖于x，若在上各阶导数存在，将展开，可将展成的幂级数形式，即，记的计算精度，可利用外推原理逐次消去式右端只要将步长h逐次分半，利用及组合消去，重复同一过程最后可得到递推公式，此时.说明用其误差阶为，这里表示m次加速。
计算时用序列表示区间分半次数，即具体计算公式为，就是Romberg求积方法。
2、程序代码：M文件1）、Romberg加速法function [s,n]=rbg1(a,b,eps)if nargin<3,eps=1e-6;ends=10;s0=0;k=2;t(1,1)=(b-a)*(f(a)+f(b))/2;while (abs(s-s0)>eps)h=(b-a)/2^(k-1);w=0;if (h~=0)for i=1:(2^(k-1)-1)w=w+f(a+i*h);endt(k,1)=h*(f(a)/2+w+f(b)/2);for l=2: kfor i=1:(k-l+1)t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);endends=t(1,k);s0=(t(1,k-1));k=k+1;n=k;else s=s0;n=-k;endend2）、改进的Romberg求积函数function [s,eer]=rbg2(a,b,eps)if nargin<3,eps=1e-6;endm=1;t(1,1)=(b-a)*(f(a)+f(b))/2;r(1,1)=0;while ((abs(t(1,m)-r(1,m))/2)>eps)c=0;m=m+1;for j=1:2^(m-1)c=c+f(a+(j-0.5)*(b-a)/2^(m-1));endr(m,1)=(b-a)*c/2^(m-1);for j=2:mfor k=1:(m-j+1)r(k,j)=r(k+1,j-1)+(r(k+1,j-1)-r(k,j-1))/(4^(j-1)-1);endt(1,j)=r(1,j-1)+2*(4^(j-2)-1)*(t(1,j-1)-r(1,j-1))/(4^(j-1)-1);endenderr=abs(t(1,m)-r(1,m))/2;s=t(1,m);3）定义f.m函数如下：function f=f(x);f=x.^3;4）运行命令及结果>> rbg1(0,2)>> rbg2(0,2)3、流程图二、圆柱体问题1、问题分析圆柱体水平方向受到地面的摩擦阻力f*Q，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。
7.数值积分-Newton-Cotes公式+龙贝格算法
c(n) k
k 0,1,2,,n，若记其绝对值的和为
n
|
c(n) k
|,
k0
则可以证明
sup{ n } .
n
（2.10）
显然，当 f ( x) 1时，对所有 n 1，都有 I ( f ) In( f ),
n
即
c(n) k
1
k0
结论：当n充分大时，
c(n) k
(k 0,1,2,,n)
当n = 8时，出现了负系数，从而影响稳定性和 收敛性，因此实用的只是低阶公式。
Newton-Cotes公式
b
a
f
( x)dx
(b
n
a)C
(n) j
f
(
x
j)
j0
• 柯特斯系数
n
1 1/2 1/2
2 1/6 4/6 1/6
3 1/8 3/8 3/8 1/8
4 7/90 16/45 2/15 16/45 7/90
但用n阶牛顿-柯特斯公式计算时会出现如下的计算结果
I
41 4 1 x2 dx 2argtan4 2.6516
n
In
2
5.4902
4
2.2776
6
3.3288
8
1.9411
10
3.5956
由上表可以看出：此时数值求积过程是发散的。
注意： 对于n阶Newton-Cotes公式的Cotes系数
（2.6）
其中：
K2(t)
1
72 1
(t (b
a t
)3 )3
(3t a (b 2a
2b), 3t ),
72
a t ab
龙贝格积分算法实验
实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式： （复化）梯形公式 11[()()]2n ii i hT f x f x-+==+∑2()12b a E h f η-''=-[,]a b η∈ （复化）辛卜生公式 11102[()4()()]6n i i i i hS f x f x f x -++==++∑4(4)()1802b a h E f η-⎛⎫=- ⎪⎝⎭ [,]a b η∈ （复化）柯特斯公式 111042[7()32()12()90n i i i i hC f x f x f x -++==+++∑31432()7()]i i f xf x +++6(6)2()()9454b a h E f η-⎛⎫=- ⎪⎝⎭[,]a b η∈ 这里，梯形公式显得算法简单，具有如下递推关系121021()22n n n i i h T T f x -+==+∑因此，很容易实现从低阶的计算结果推算出高阶的近似值，而只需要花费较少的附加函数计算。
但是，由于梯形公式收敛阶较低，收敛速度缓慢。
所以，如何提高收敛速度，自然是人们极为关心的课题。
为此，记0,k T 为将区间[,]a b 进行2k等份的复化梯形积分结果，1,k T 为将区间[,]a b 进行2k等份的复化辛卜生积分结果，2,k T 为将区间[,]a b 进行2k等份的复化柯特斯积分结果。
根据李查逊（Richardson ）外推加速方法，可得到1,11,,0,1,2,40,1,2,41m m k m km k m k T T T m -+-=-⎛⎫=⎪=-⎝⎭可以证明，如果()f x 充分光滑，则有,lim ()m k k T I f →∞= （m 固定）,0lim ()m m T I f →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
T3,2
4T3,1 T2,1 3
0.916787624
T3,3
16T3,2 T2,2 3
0.916478228
具体结果见下表
k Tk ,1
1 1.05
Tk , 2
Tk , 3
Tk , 4
I 0.9162907318L
2 0.953571429 0.921428571
3 0.925983575 0.916787624 0.916478228
梯形加速公式：
Sn
4T2n Tn 3
4 4 1T2n
1 4 1Tn
上述公式说明：
利用复化梯形公式前后两次积分近似值 Tn和 T2n，按
照上式作出的线性组合得到了具有更高精度的积分值。
Romberg积分公式正是由此思想产生
类似于梯形加速公式的处理方法，得到：
I
S2n
1 42 1(S2n
Sn )
42 S2n Sn 42 1
0.9461459 0.9460869 0.9460834
C2k 2
K 2k 3
0.9460830 0.9460831
0.9460831
例3：利用Romberg 积分法式计算积分I 1.5 1 dx
要求精确到小数点后面7位。
0 x 1
解： f ( x) 1 x 1
根据Romberg 积分法计算得
解：由题意
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)]}
Cotes值序列 C2k Romberg值序列 R2k
S2k
4T2k1 T2k 41
C2k
42 S2k1 S2k 42 1
R2k
4 C 3 2k 1
C2k
43 1
Romberg积分法的一般公式
Tm, j
4
T j 1 m,
j 1
Tm 1,
j 1
4 j1 1
其中
j 2,3, 4;m j
注意事项：
（1）使用复化梯形公式、Simpson公式，首先要确定步长 h ；
（2）而步长要根据余项确定，这就涉及到高阶导数的估计； （3）高阶导数的估计一般比较困难，且估计值往往偏大； （4）计算机上实现起来不方便，通常采用“事后估计法”。
三、积分步长的自动选取： ➢基本思想： 将积分区间逐次分半
➢终止法则： 前后两次近似值的误差小于已知精度
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
把区间再二分，重复上述步骤算得
T16=3.140941613 S8=3.141592652 C4=3.141592662 R2=3.141592640 由于 |R1-R2|<=0.00001,计算可停止，取
011+4x2 dx R2=3.14159
➢ Romberg积分思想
由上节分析知，用复化梯形公式计算积分值 I
T2n 的误差大约为：
1 3 (T2n Tn )
令
I
T2n
1 3 (T2n
Tn )
4T2n Tn 3
由复化梯形公式知
1
b a n1
T2n 2 Tn
2n
k0
f
(
x
k
1
)
2
4T2n Tn
3
1 3 Tn
2(b a) n1 3n k0
主讲人：
I2n In
➢具体过程（以复化梯形公式为例）
1、首先将区间 [a, b]n等分： h b a
n
h
n1
Tn
2
f (a) 2
k 1
f ( xk )
f (b)
h
2、再将区间[a, b]2n等分，即步长减半： h1 2
T2n
h1 2
f
(a)
n1
2
k 1
n1
f ( xk ) 2
令
h
b
2i
a
(i 0,1,2, )
计 （算 3）：：T2n按 12加Tn速 h2公ni01式f (求xi12积) 分n（ 2为i 便于编程，写
为下列形式）
梯形加速：Sn T2n (T2n Tn)/3 辛普生加速：Cn S2n (S2n Sn) /15 柯特斯加速： Rn C2n (C2n Cn ) / 63 （4）：精度控制
1.5
T1,1
[ f (0) f (1.5)] 1.05 2
1 T2,1 2 [T1,1 1.5 f (0.75)] 0.953571429
T2,2
4T2,1 T1,1 3
0.921428571
T3,1
1 2 {T2,1
0.75[
f
(0.75)
f
(1.125)]}
0.925983575
f
(
x
k
1
)
2
b a
n
n1
6n
f (a) 2
k 1
f ( xk )
f (b) 4
k0
f
(
x
k
1 2
)
S b a
n1
n
Tn
2n
f (a) 2 f (xk ) f (b)
k 1
h
n1
n1
Sn (
f
)
6
f
(a)
4
k0
f
(
x
k
1 2
)
2
k
1
f
(xk )
f
(b)
由此得到近似关系式
I
T2n
1 4 1(T2n
Tn )
误差控制条件
1 4 1 (T2n
Tn )
上述条件满足，程序终止；否则，继续分半计算。
➢对于复化梯形公式
1 I T2n 4 1(T2n Tn )
➢对于复化Simpson公式、Cotes公式可以类似得到
I
S2n
1 42 1(S2n
Sn
)
1 I C2n 43 1(C2n Cn )
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 值序列
通过上述3个积分值序列求积分近似值的方法， 称之为Romberg积分法。
4个积分值序列：
梯形值序列 T2k Simpson值序列 S2k
k0
f
(
x
k
1
)
f
2
(b)
1
b a n1
2 Tn
2n
k0
f
(
x
k
1
)
2
f
(
x
f
( xk
) 2
3、终止条件： 由复化梯形公式的余项知
f ( x) 变化不大时
I
Tn
ba 12
(b
a )2 n
f
(1 )
I
T2n
ba 12
(b a 2n
)2
f
(2 )
I Tn 4 I T2n
应用步长逐次减半得到的复化梯形值、复化 Simpson值、复化Cotes值与精确值的比较
S I
T2n
1 4 1(T2n
Tn
)
n
I
S2n
1 42 1(S2n
Sn
)
Cn
加速 收敛
1
I C2n 43 1(C2n Cn ) Rn
Romberg积分法/*Romberg Integration Method */
=3.131176471 S2=1/3(4T4-T2)=3.141568628 C1=1/15(16S2-S1)=3.142117648 计算f(1/8) f(3/8) f(5/8) f(7/8)进而求得 T8=1/2{T4+1/4[f(1/8)+f(3/8)+f(5/8)+f(7/8)]}
=3.138988495 S4=1/3(4T3-T4)=3.141592503 C2=1/15(16S4-S2)=3.141594095 R1=1/63(64C2-C1)=3.141585784
当 R2n Rn , (为精度）时终止计算 并取 R2n为近似值，否则，将步长折半，转 （2）执行。实际计算时的加工流程图7.4.1 示：
例7.4.2 用龙贝格算法加工例7.4.1得到的近似 值
解：用算法图7.4.2计算结果见表7.4.1,表中
k表示二分次数 ：
k
T2k
S 2k 1
0 0.9207355 1 0.9397933 2 0.9445135 3 0.9456909
Tm,1 T2m1 (m 1) Tm,3 C2m3 (m 3)
Tm,2 S2m2 (m 2) Tm,4 R2m4 (m 4)