Romberg求积分公式.docx
《Romberg求积分公式.docx》由会员分享，可在线阅读，更多相关《Romberg求积分公式.docx（16页珍藏版）》请在冰豆网上搜索。
Romberg求积分公式
《MATLAB程序设计实践》课程考核
1、编程实现以下科学计算算法，并举一例应用之。
“Romberg求积分公式”
2、编程解决以下科学计算和工程实际问题。
1）、给定半径的为r，重量为Q的均质圆柱，轴心的初始速度为v0，初始角速度为w0且v0>r*w0，地面的摩擦系数为f，问经过多少时间后，圆柱将无滑动地滚动，求此时的圆柱轴心的速度。
2）、在一丘陵地带测量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟合一曲面确定合适的模型，并由此找出最高点和该点的高程。
100
200
300
400
100
636
697
624
478
200
698
712
630
478
300
680
674
598
412
400
662
626
552
334
一、Romberg求积分公式
1、算法说明：
此算法可自动改变积分步长,使其相临两个值的绝对误差或相对误差小于预先设定的允许误差.Romberg加速法公式
在等距节点的情况下，通过对求积区间（a，b）的逐次分半，由梯形公式出可逐次提高求积公式精度，这就是Romberg求积的基本思路，由于梯形公式余项只有
精度，即
，但当节点加密时可组合成
其精度达到
，如果再由
与
组合成
则可使误差精度达到
，于是
依赖于x，若
在
上各阶导数存在，将
展开，可将
展成
的幂级数形式，即
，记
的计算精度，可利用外推原理逐次消去式
右端
只要将步长h逐次分半，利用
及
组合消去
，重复同一过程最后可得到递推公式
，此时
.说明用
其误差阶为
，这里
表示m次加速。
计算时用序列
表示区间分半次数，即
具体计算公式为
，就是Romberg求积方法。
2、程序代码：
M文件
1）、Romberg加速法
function[s,n]=rbg1（a,b,eps）
ifnargin<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）
fori=1:
（2^（k-1）-1）
w=w+f（a+i*h）;
end
t（k,1）=h*（f（a）/2+w+f（b）/2）;
forl=2:
k
fori=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;
elses=s0;
n=-k;
end
end
2）、改进的Romberg求积函数
function[s,eer]=rbg2（a,b,eps）
ifnargin<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;
forj=1:
2^（m-1）
c=c+f（a+*（b-a）/2^（m-1））;
end
r（m,1）=（b-a）*c/2^（m-1）;
forj=2:
m
fork=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）定义函数如下：
functionf=f（x）;
f=x.^3;
4）运行命令及结果
>>rbg1（0,2）
>>rbg2（0,2）
3、流程图
否
是
否
是
是
否
否
是
否
是
输入左右端点a，b
输入精度值eps
nargin>3
eps<1e-6
S=10；S0=0；K=2
t（1,1）=（b-a）*（f（a）+f（b））/2
abs（s-s0）>eps
h=（b-a）/2^（k-1）w=0
h~=0
1<=i<=（2^（k-1）-1）
w=w+f（a+i*h）
t（k,1）=h*（f（a）/2+w+f（b）/2）
2<=l<=k
1<=i<=（k-l+1）
否
s=s0;n=-k;
流程图1
i=i+1
l=l+1
i=1
=
l=2
开始
i=1
t（i,l）=（4^（l-1）*t（i+1,l-1）-t（i,l-1））/（4^（l-1）-1）;
s=t（k,1）;s0=（t（k-1,1））;k=k+1;n=k;
输出结果
结束
i=i+1
是
是
否
是
否
是
是
否
否
是
输入左右端点a,b
输入精度值eps
nargin>3
eps<1e-6
m=1
t（1,1）=（b-a）*（f（a）+f（b））/2
（（abs（t（1,m）-r（1,m））/2）>eps
c=0;m=m+1;j=1
1<=j<=2^（m-1）
j=2
否
j<=m
流程图2
j=j+1
r（1,1）=0
c=c+f（a+*（b-a）/2^（m-1））
j=j+1
r（m,1）=（b-a）*c/2^（m-1）
k=1
k<=m-j+1
r（k,j）=r（k+1,j-1）+（r（k+1,j-1）-r（k,j-1））/（4^（j-1）-1）
k=k+1
t（1,j）=r（1,j-1）+2*（4^（j-2）-1）*（t（1,j-1）-r（1,j-1））/（4^（j-1）-1）
输出结果
结束
err=abs（t（1,m）-r（1,m））/2
s=t（1,m）
二、圆柱体问题
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='）;
ifv0j=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=[100200300400];
y=[100200300400]';
z=[636697624478;698712630478;680674598412;662626552334];
xx=linspace（100,400,N）;
yy=linspace（100,400,N）';
zh=interp2（x,y,z,xx,yy,'cubic'）
mesh（xx,yy,zh）
s=0;
fori=1:
N^2
ifzh（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=
Columns1through7
………………………………
Zh为100*100的矩阵，此处只列出部分值，其余省略。
并解得最高点和最高程为：
x0=
166
y0=
196
ans=
即最高点在坐标（166,196）处，且高程为。
升级会员 