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 ················
- 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）处，且高程为。