Romberg求积分公式
龙贝格求 积分
龙贝格(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
Romberg算法
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
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
梯形法递推公式
Cha4.2.m
程序演示
Cha41.m: 复化N-C公式 Cha42.m: 梯形法的递推
龙贝格 (Romberg) 加速
梯形法递推公式算法简单，编程方便
但收敛速度较 慢 解决方法：龙贝格 (Romberg) 加速
复化积分公式的渐近状态
思想：利用余项公式与积分的定义
龙贝格求积公式
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)]}
Romberg数值积分
Romberg数值积分实验报告一、实验目的1、加深外推法的原理理解, 掌握Romberg外推法的计算方法。
2、用matlab软件实现Romberg数值积分来计算题目的运算。
二、基本理论及背景1、理论推导：对区间[a， b]，令h=b-a构造梯形值序列{T2K}。
T1=h[f(a)+f(b)]/2把区间二等分，每个小区间长度为h/2=(b-a)/2，于是T2 =T1/2+[h/2]f(a+h/2）把区间四(2)等分，每个小区间长度为h/2 =（b-a）/4，于是T4 =T2/2+[h/2][f（a+h/4）+f（a+3h/4).....................把[a，b] 2等分，分点xi=a+（b-a）/ 2 ·i （i =0，1，2 · · · 2k）每个小区间长度为（b-a）/ 2 .2、参考Romberg数值积分，实现积分的数值求解，完成下列题目：三、算法设计及实现1、算法设计(a)function y=fun1(x)y=sqrt(4-(sin(x)).^2)；(b)function y=fun2 (x)y=log(1+2*x)/3;(c)function y=fun3(x)y=exp(-x).*sin(x.^2)四、实验步骤1、○1打开matlab软件，新建Romberg.m文件,在窗口中编辑Romber数值积分函数程序代码,并保存在指定的文件夹下,在Current Directory窗口右边点击《Browse For Folder》按钮指向Romberg.m文件；○2在Command Window中编辑相应要计算的题目的数值函数及相应的题目的表达式。
2、输出结果和初步分析说明（见附件一）。
五、使用说明实验结果分析1、在Command Window窗口中编辑要调用的函数名与指定的函数名字不同导致出现错误，通过改正与函数名相同即可。
2、Romberg方法也称为逐次分半加速法。
龙贝格(Romberg)求积法
龙贝格（Romberg）求积法
复化求积方法对于提高计算精度是行之有效的 方法，但复化公式的一个主要缺点在于要先估计出 步长。若步长太大，则难以保证计算精度，若步长 太小，则计算量太大，并且积累误差也会增大。在 实际计算中通常采用变步长的方法，即把步长逐次 分半，直至达到某种精度为止。
1.1变步长的梯形公式 变步长复化求积法的基本思想是在求积过程中，
)
0.9397933
进一步二分求积区间,并计算新分点上的函数值
f
(
1 4
)
0.9896158
,
f
(
3 4
)
0.9088516
有
T4
1 2 T2
1 4
f
(
1 4
)
f
(
3 4
)
0.9445135
这样不断二分下去，计算结果如P139列表所 示。积分的准确值为0.9460831，从表中可
看出用变步长二分10次可得此结果。
的考积察分T值与nS等n 。份辛卜生公式 S n之间的关系。将
复化梯形公式
Tn
h 2
f
(a)
2
n1 k 1
f
(xk ) f (b)
梯形变步长公式
T2 n
Tn 2
h n1
2 k 0
f (xk1 ) 2
代入(6.9) T 表达式得
h
n1
n1
T
6 f (a) 4k0
f
(
x k
1
)
2
（
2
输入 a,b,ε
）
变
b-ah,
h 2
f(a)+f(b)T1
7.3Romberg积分(共64张PPT)
第五页，共六十四页。
表7-1
O(h)
O(h2 )
O(h3 )
O(h4 )
例1 设
带余项的差分(chà fēn)公式为
第六页，共六十四页。
导出具有误差为 解令
用 h/2代替(dàitì)h,得
(11) 的外推公式.
h2
(12)
为消去含 h2的项,用4乘(12)式减去 (11)式,得
第七页，共六十四页。
第二页，共六十四页。
把(1)式改写(gǎixiě)为 (3)
用h/2代替(3)式中的h,得
(4) 用2乘(4)式再减去(3)式,消去含h的项,得
(5)
令
,且记
第三页，共六十四页。
那么(5)式可写为
(6) 这里, 逼近 的误差为
再用 h/2 代替(dàitì) h , 使(6)式变为
(7) 用4乘(7)式减去(6)式,消去含 的项,得
则所谓 gn (x)与 的根相同,即是指这两个正
交多项式的根有如下的关系.
第三十二页，共六十四页。
性质5 (1) 区间[a,b]上带权函数 (的x) 正交
多项式序列
与
对应元素之间
只相差一个比例常数.
(2)区间[a,b]上带权函数
首项 系数 (x)
(shǒu xiànɡ)
为1的
正交多项式序列
唯一.
常见的正交多项式有Legendre(勒让德)多
7.5.1 引言 求积公式
当求积系数
、求积节点
(1) 都可以
自由选取时,其代数精确度最高可以达到多少 次?
下面的引理可以回答上述(shàngshù)问题.
第二十四页，共六十四页。
引理1 当求积系数
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 ················
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 计算的近似值。
数学实验题目2 Romberg积分法
数学实验题目2 Romberg 积分法摘要考虑积分()()b aI f f x dx =⎰欲求其近似值，可以采用如下公式：（复化）梯形公式 110[()()]2n i i 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 →∞=这是一个收敛速度更快的一个数值求积公式，我们称为龙贝格积分法。
龙贝格求 积分
龙贝格（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++k k h 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求积分公式
《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，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。
Romberg求积法
安徽中医药大学题目：FOnberg求积法c语言编程姓名：杨撞撞学号：13713042班级：13医软（1）班目录1简介2计算公式3算法描述4程序流程图5算法程序表示6算法结果截图1. 简介龙贝格求积公式也称为逐次分半加速法。
它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。
作为一种外推算法，它在不增加计算量的前提下提高了误差的精度.在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。
这样，前一次分割得到的函数值在分半以后仍可被利用，且易于编程。
2. 计算公式梯形公式----- 二复化辛普森公式------ A复化科特斯公式V龙贝格求积公式其对应的公式为：T2n=1/2（Tn+Hn）（梯形公式）Sn=4/（4-1）T2n-1（4-1）Tn （辛普森公式）Cn=4八2/（4八2-1）S2n-1/（4八2-1）Sn （柯特斯公式）Rn=4八3/（4八3-1）C2n-1/（4八3-1）Cn （龙贝格求积法公式）3. 算法描述3.1龙贝格算法基本描述先算出T0 （ 0）,从而计算出T0（1）,以此类推，直到计算出|T0（0）-Tn-1 （0）|<e即可利用加速推算公式推算出结果。
3.2龙贝格算法程序包步骤1.输入积分上限2 输入程序下限3输入区间等分数4输入要求的函数5计算出所求函数的积分，分别是：复化梯形求积结果辛普森求积结果柯特斯求积结果龙贝格求积结果4. 程序流程图例题：用Romberg 方法计算积分I= 01sin （x）/xdx 的相关算法流程图表示如下图开始 欢迎界面择函 择/——输入1——输入信息：积分:卜线，等分段积.4.x2*sinxs Ynx 2* cosx Yx *sin 2xYsi nx/xX合法？3 输入0Y*辛普林 求积科特斯 求积龙贝格 求积输出结 果5.算法程序表示#include<stdio.h> #include<math.h> #defi ne A(x)(si n(x)/x)//宏定义若干常用函数A,B,C,D,E,G#define B(x)(cos(x*x+2*x+1))#define C(x)(atan(sqrt(x*x+1)))#define D(x)(sqrt(exp(x)+sin(2*x)))#define E(x)(x*x*x+3*x*x+5)#define G(x)(log10(x)/pow(2,x))double t[20],s[20],c[20],r[20];// 定义全局数组double dh,fan,a,b,m; // 定义全局变量int jj=0;char hs;double F(double x) // 用switch 调用若干被积函数{switch(hs){case 'A':fan=A(x);break;case 'B':fan=B(x);break;case 'C':fan=C(x);break;case 'D':fan=D(x);break;case 'E':fan=E(x);break;case 'G':fan=G(x);break;default :printf(" 输入错误！");}return(fan);// 返回被积函数值} double H(int i) // 求和函数并返回和SUM{int j;double zh,SUM=0.0; II定义求和变量SUM并赋初值for(j=1;j<=pow(2,i-1);j++){ zh=(a+((2*j-1)*(b-a))Ipow(2,i));SUM=SUM+F(zh); II 调用F(x) 函数}SUM=(b-a)*SUMIpow(2,i);return(SUM);}double Txing(int k) II 梯形公式if(k==O) dh=t[jj]=((b-a)/2)*(F(a)+F(b));//分半次数为零时T 形公式求积 else {dh=0.5*Txi ng(k-1)+H(k); //Txing函数递归调用循环输出并返回dht[++jj]=dh; }m=pow(2,jj);prin tf("T[%0.0lf]=%0.7lf\t",m,t[jj]);//输出并返回dhreturn(dh); }double Simpson(int k) // 辛普森公式 错误！未找到引用源int i,j;b —俱 V 1.「 "‘+弓厂『[°十(雄-耳{廉G +甩)],Txin g(k); // 调用梯形公式for(i=0;i<=k-1;i++)s[i]=(4.0*t[i+1]-t[i])/3.0; II 森公式prin tf("\n");for(j=0;j<=k-1;j++) //{m=pow(2,j);prin tf("S[%0.0lf]=%0.7lf",m,s[j]);prin tf("\t");}return(s[k-1]); // 返回最后一个值}double Cotes(i nt k) //{int i,j;Simps on (k);for(i=0;i<=k-2;i++)c[i]=(16.0*s[i+1]-s[i])/15.0; //递推辛普循环输出s[k-1]科特斯公式递推科特斯公式prin tf("\n");for(j=0;jv=k-2;j++) II环输出{m=pow(2,j);prin tf("C[%0.0lf]=%0.7lf\t",m,c[j]);}return(c[k-2]); // 返回最后一个值}double Romberg(i nt k) //公式」-亠厂{int i,j; //斯公式Cotes(k);for(i=0;i<=k-3;i++) r[i]=(64.0*c[i+1]-c[i])/63.0; //隆贝格公式prin tf("\n");c[k-2]隆贝格调用科特递推for(j=0;j<=k-3;j++) // 环输出{m=pow(2,j); printf("R[%0.0lf]=%0.7lf\t",m,r[j]);printf("\n");}return(r[k-3]); // 返回最后一个值r[k-3] }main(){int k;char y;printf(" 请从以下公式中选择积分函数：\n"); printf("A: sin (x)/(x)\tB:cos(xA2+2x+1)\tC:ata n( sqrt(x A2+1))\n\n");tG:log10(x)/pow(2,x)\n\n");printf(" 请选择函数F(x)(大写)：\n");scanf("%c",&hs);。
微积分的数值计算方法Romberg算法
则由(1)(2)(3)式,有 梯形公式 1次二等分后 的梯形值
ba T1 [ f (a) f (b)] T0 (0) 2 1 ba 1 T2 T1 f (a h) T0 (1) 2 2 2
k 1
ba 若n 2 , h k 1 2 记Tn T0 (k 1)表示k 1次二等分后的所求的梯形值 ba x j a jh a j k 1 2 1 x 1 x j h a ( j 1 ) b a a (2 j 1) b a k 1 k j 2 2 2 2 2
k 1
梯形公式
第k次二等分后所求的梯形值
第k-1次二等分后所求的梯形值
上式称为递推的梯形公式
二、加速公式(误差补偿手段)
由复化梯形公式的余项公式:
I Tn ch ,
2
I T2 n 1 I Tn 4
h I T2n c 2
2
3I 4T2 n Tn 0
因此(1)(2)(3)式可化为如下递推公式
ba [ f (a ) f (b)] T0 (0) 2
2 1 1 b a T0 (k ) T0 (k 1) k f (a (2 j 1) b k a ) 2 2 j 0 2 -------(4) k 1, 2,
T0 (2)
T2 (0)
T2 (1)
T0 (3)
T1 (2)
T3 (0)
n 1 4 1 b a 1 * I T ( Tn f ( x 1 )) Tn j 3 2 2n j 0 3 2
1 4(b a) n1 Tn f (x 1 ) j 3 6n j 0 2
Romberg积分及初值问题
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
--------(3)
3
n 1时，h b a
则由(1)(2)(3)式,有
ba T1 [ f ( a ) f (b)] 2 1 ba 1 T2 T1 f ( a h) 2 2 2
T0 (0) T0 (1)
k 1, 2 ,
若n 2 k 1
记Tn T0 (k 1)
1 b a 2 1 ba T0 (k ) T0 (k 1) k f (a (2 j 1) k ) 2 2 2 j 0 k 1, 2 ,
k 1
ba T0 (0 ) [ f ( a ) f (b )] 2
上式称为递推的梯形公式
5
三种公式之间的关系 （3，-3）（-1,2）（0,1）（1，-1） 4T2 n Tn （3，-4） Sn 4 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)
5.3龙贝格(Romberg)算法
一、复合梯形公式的递推化
将定积分I = ∫ f ( x )dx的积分区间[ a , b ]分割为n等份 a b−a xk = a + jh , j = 0 ,1,L , n 各节点为 h= n 复合梯形(Trapz)公式为
n−1 b−a Tn = [ f ( a ) + 2 ∑ f ( x j ) + f (b )] 2n j =1
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; }
13
其中外推加速公式可简化为
1 Tm ( k − 1) = m [ 4 m Tm − 1 ( k ) − Tm − 1 ( k − 1)] 4 −1
--------(11)
并且m可以推广到 m = 1,2 ,L
Romberg算法求解步骤
k = 1 , 2 ,L
Romberg算法的代 数精度为m的两倍 Romberg算法的收敛 阶高达m+1的两倍
第6章 数值积分与数值微分
龙贝格(Romberg)算法 算法 龙贝格
1
龙贝格(Romberg)算法 算法 龙贝格
综合前几节的内容,我们知道 梯形公式,Simpson公式,Cotes公式的代数精度分别为 1次,3次和5次 复合梯形、复合Simpson、复合Cotes公式的收敛阶分别为 2阶、4阶和6阶 无论从代数精度还是收敛速度,复合梯形公式都是较差的 有没有办法改善梯形公式呢?
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《MATLAB程序设计实践》课程考核
1、编程实现以下科学计算算法，并举一例应用之。
“Romberg求积分公式”
2、编程解决以下科学计算和工程实际问题。
1）、给定半径的为r，重量为Q的均质圆柱，轴心的初始速度为v0，初始角速度为w0且v0>r*w0，地面的摩擦系数为f，问经过多少时间后，圆柱将无滑动地滚动，求此时的圆柱轴心的速度。
2）、在一丘陵地带测量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。
100200300400 100636697624478
200698712630478
300680674598412
400662626552334
一、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;end
s=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);
end
t(k,1)=h*(f(a)/2+w+f(b)/2);
for l=2: k
for i=1:(k-l+1)
t(i,l)=(4^(l-1)*t(i+1,l-1)-t(i,l-1))/(4^(l-1)-1);
end
end
s=t(1,k);
s0=(t(1,k-1));
k=k+1;
n=k;
else s=s0;
n=-k;
end
end
2）、改进的Romberg求积函数
function [s,eer]=rbg2(a,b,eps)
if nargin<3,eps=1e-6;end
m=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+*(b-a)/2^(m-1));
end
r(m,1)=(b-a)*c/2^(m-1);
for j=2:m
for 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);
end
t(1,j)=r(1,j-1)+2*(4^(j-2)-1)*(t(1,j-1)-r(1,j-1))/(4^(j-1)-1);
end
end
err=abs(t(1,m)-r(1,m))/2;
s=t(1,m);
3）定义函数如下：
function f=f(x);
f=x.^3;
4）运行命令及结果
>> rbg1(0,2)
>> rbg2(0,2)
3、流程图
二、圆柱体问题
1、问题分析
圆柱体水平方向受到地面的摩擦阻力f*Q，该摩擦力对轴心速度起减速作用，同时又产生一个力矩，对角速度起加速作用。
综上，等到轴心速度v=w*r时，圆柱体将无摩擦运动。
2、源程序：M文件
r=input('r=');
Q=input('Q=');
g=input('g=');
f=input('f=');
v0=input('v0=');
w0=input('w0=');
if v0<r*w0;end;
j=Q*r^2/2/g;
F=f*Q;
beta=F*r/j;
a=-F/(Q/g);
t=(v0-w0*r)/(beta*r-a)
v=v0+a*t
3、执行命令
>> move
r=1
Q=100
g=
f=
v0=3
w0=2
t =
v =
三、高程
1、题中已给出4*4=16个数据，分别对应16个坐标位置上的高程，现只需采用插值的方法，向其中填补数值，便可拟合对应的曲面，考虑到找到适合曲面的二元函数比较复杂，并且插值之后的数据量够大（10000个），具有一定的代表性，因此在求解丘陵最高点及其高程的时候，可以将所有数据进行比较，取其最大值所对应的x，y值作为最高点。
具体程序如下
2、源程序：M文件
function [s,x0,y0]=high(N)
x=[100 200 300 400];
y=[100 200 300 400]';
z=[636 697 624 478;698 712 630 478;680 674 598 412;662 626 552 334];
xx=linspace(100,400,N);
yy=linspace(100,400,N)';
zh=interp2(x,y,z,xx,yy,'cubic') mesh(xx,yy,zh)
s=0;
for i=1:N^2
if zh(i)>s
s=zh(i);
n=mod(i,N);
m=(i-n)/N;
end
end
x0=100+300/N*m
y0=100+300/N*n
3、执行命令
>> high(100)
运行结果如下：
作出拟合曲面为
求得结果：
zh =
Columns 1 through 7
………………………………
Zh为100*100的矩阵，此处只列出部分值，其余省略。
并解得最高点和最高程为：
x0 =
166
y0 =
196
ans =
即最高点在坐标（166,196）处，且高程为。