Romberg求积分公式

《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

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)';