龙贝格求积公式
龙贝格求 积分
龙贝格(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
龙贝格求积公式
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)]}
龙贝格算法
2）理查森外推加速 从梯形公式出发， 将区间[a,b]逐次二分可以提高求积公式精度， 当[a,b] 分为 n 等份时，若记T������ = T ℎ ，当区间[a,b]分为 2n 等份时，则有 T2������ = T
ℎ 2
。再有泰勒公式展开为： T ℎ = ������ + ������1 ℎ2 + ������2 ℎ+ ⋯ +
建立一个命名为Romberg.m的function文件：
function[T]=Romberg(f,a,b,e) T=zeros(10,10); T(1,1)=(b-a)/2*(f(a)+f(b)); for k=2:10, sum=0; for i=1:2^(k-2), x=a+(2*i-1)*(b-a)/2^(k-1); sum=sum+f(x); end T(k,1)=T(k-1,1)/2+[(b-a)/2^(k-1)]*sum; for j=2:k, end %第一列用递推梯形公式 %定义龙贝格函数 %定义10阶的零元矩阵
1 3/2 ������ 0
������������。
算，并取������ ������, ������ ≈ ������ ；否则令 k=k+1，转（3）继续计算。 6）下图为我按照自己的算法所设计的示意表： 算法设计表：
k 1 2 3 4 …
h b-a (b-a)/2 (b-a)/4 (b-a)/8 …
������ ������, 1 T(1,1) T(2,1) T(3,1) T(4,1)
������ ������, 2
������ ������, 3
������ ������, 4
������ ������, 5
数值分析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 的高阶导数，由于
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
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 ················
龙贝格积分算法实验
实验题目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 →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
该方法的计算可按下表进行0,0T 0,1T 0,2T … 0,m T 1,0T 1,1T … 1,1m T - 2,0T … 2,2m T - … … ,0m T很明显，龙贝格计算过程在计算机上实现时，只需开辟一个一维数组，即每次计算的结果,m k T ，可存放在0,k T 位置上，其最终结果,0m T 是存放在0,0T 位置上。
4-4Romberg积分法
T ( h) I a1h2 a2 h4 a3 h6 ... ak h2 k ...
式中系数 ak ( k 1,2,...) 与h无关.
李查逊外推加速法的处理过程：
由 T ( h) I a1h2 a2 h4 a3 h6 ... ak h2 k ... 可知T(h)≈I 是二阶收敛的 减
对辛甫生公式加工： 对柯特斯公式加工：
实质：将粗糙的梯形公式值逐步加工成精度较
高的公式。
现在的问题是：
能否对龙贝格公式再加工 取得较高精度的公式？
李查逊（Richardson）
外推加速法
三、李查逊(Richardson)外推加速法
李查逊外推加速法基于如下原理: 定理 设
T ( h) Tn 则成立 f ( x ) C [a , b ] ，
n1
只需计算新增分点的 函数值
dx 至具有7位有 举例 计算积分值 I 0 x 效数字。
详见书上p110例5,计算时要注意公式中步长的含义.
1 sin
x
注：本例提供了精确值供作比较； 事后估计法；事先估计法
请观察：本例最后结果二分多少次？ 共有多少个分点？ 答：二分10次,共有210 +1=1025个分点,将区间等分1024份。
若构造
则又可以进一步消去展开式中的 h4 项，而有
T2 ( h) I 1h6 2 h8 ...
这样构造出的 T2(h) ，其实就是柯特斯公式， 它与积分值I的逼近阶为六阶。 如此继续下去，每加速一次，误差的量级 便提高2阶，这就是李查逊外推加速法。
定义：利用李查逊（Richardson）外推加速法将低 阶求积公式通过加密节点，再加工成高阶的求积公 式，统称为龙贝格（Romberg）公式。
龙贝格求积法
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 计算的近似值。
(龙贝格求积)
针对这类问题的算法技巧是在不同区间上预测被积函数 变化的剧烈程度确定相应的步长.
这种方法称为自适应积分方法.
设给定精度要求 0，计算积分
b
I a f (x)dx
的近似值.先取步长 h b a ，应用辛普森公式有
I b f (x)dx S(a, b) b a ( h )4 f (4) (),
注意到每个子区间等分分为设将区间一次如果将求积区间再二分1个分点共有子区间上的积分值为用复化梯形公式求得该分点经过二分只增加了一个注意到不难导出下列递推公式逐次分半算法据梯形公式计算得9588510求出将区间二等分进一步二分求积区间9445135计算积分值8414709例题1精确至三位有效数字例题25695注意到57例如
例如：对区间[a, b], 令h b a构造梯形值序列{T2k }
h
T1
[ 2
f (a)
f
(b)]
逐次分半算法
把区间二等分，每个小区间长度为 h/2=(b-a）/2，
于是
T2 =T1/2+[h/2]f(a+h/2） 把区间四(22)等分，每个小区间长度为h/2 2 =（b-a）
/4，于是
T4 =T2/2+[h/22][f(a+h/4)+f(a+3h/4)]
················
把[a,b] 2k 等分，分点xi=a+(b-a)/ 2k ·i （i =0， 1，2 ···2k）每个小区间长度为(b-a)/ 2k ，由归
纳法可得
T T 1
2k 2
2 k1
ba 2k
注意到每个子区间
[ xk , xk1 ]经过二分只增加了一个分点
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，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
// 龙贝格积分公式.cpp : 定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include<iostream>
#include<math.h>
#define F(x) (exp(x)*cos (x)) //函数举例。
using namespace std;
//------------步长及4，16，64.......的实现----------
double pf (int i)
{
int s=1;
for (int j=0;j<i;j++)
s*=2;
return s/2;
}
//---------定义一个求t1,t2......的函数-------------
double gcc (double a, double b,double aa[][20],int i)
{
double s,h,x;
h=(b-a)/pf(i);
s=0;
x=a+h/2;
do
{
s+=F(x);
x+=h;
}
while (x<b);
aa[0][i]=aa[0][i-1]/2+h/2*s;
return 0;
}
//----------------------主函数---------------------------
int main()
{
double aa[20][20]={0},e,a,b,h;
int j,i,n;
cout <<"请输入积分区间：\na= ";
cin >>a;
cout <<"b= ";
cin >>b;
cout <<"请输入精度：e=";
cin >>e;
aa[0][0]=(b-a)*(F(a)+F(b))/2.0;
gcc(a,b,aa,1);
aa[1][0]=(4*aa[0][1]-aa[0][0])/3;
for (i=1;i<20;i++)
{
gcc(a,b,aa,i+1);//求下一个要用的t。
for ( n=i,j=1; n>0;n--,j++)//加速公式的实现。
aa[j][n]=(pf(2*j+1)*aa[j-1][n+1]-aa[j-1][n])/(pf(j*2+1)-1);
if (fabs(aa[i][1]-aa[i][0])<e)//判断是否达到精度。
break;
else
aa[i+1][0]=(pf(2*i+3)*aa[i][1]-aa[i][0])/(pf(2*i+3)-1);
}
cout <<"积分结果为："<<aa[i][1];
return 0;
}。