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