Romberg求积分公式.docx
- 文档编号:25740863
- 上传时间:2023-06-12
- 格式:DOCX
- 页数:16
- 大小:148.36KB
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=');
ifv0 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=[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)处,且高程为。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Romberg 积分 公式
![提示](https://static.bdocx.com/images/bang_tan.gif)