公式和关系Romberg求积分公式

《MATLAB程序设计实践》课程考核
1、编程实现以下科学计算算法，并举一例应用之。“Romberg求积分公式”
2、编程解决以下科学计算和工程实际问题。
1）、
给定半径的为r，重量为Q的均 质圆柱，轴心的初始速度为v0，初始角速度为w0
且v0>r*w0，地面的摩擦系数为f，问经过多 少时间后，圆柱将无滑动地滚动，求此时
的圆柱轴心的速度。
2
）、在一丘陵地带测 量高程，x和y方向每隔100m测一个点，得高程数据如下，试拟
合一曲面确定合适的模型，并由此找 出最高点和该点的高程。

100
200
300
400
100
636
698
680
662
200
697
712
674
626
300
624
630
598
552
400
478
478
412
334

一、Romberg求积分公式
1、算法说明 ：此算法可自动改变积分步长,使其相临两个值的绝对误差或相对误差小于预
先设定的允许误差.Rom berg加速法公式

在等距节点的情况下，通过对求积区间（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+(j-0.5)*(b-a)2^(m-1));

end

r(m,1)=(b-a)*c2^(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）定义f.m函数如下：
function f=f(x);

f=x.^3;

4）运行命令及结果
>> rbg1(0,2)
>> rbg2(0,2)

3、流程图






流程图1
开始
输入左右端点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？
是
i=1
否
1<=i<=(2^(k-1)-1)
i=i+1
是
w=w+f(a+i*h)
t(k,1)=h*(f(a)2+w+f(b)2)
l=2
2<=l<=k
是
l=l+1
否
i=1
1<=i<=(k-l+1
)
否
否
s=s0;
n=-k;

是
i=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
输出结果
结束
输入左右端点a,b
流程图2
输入精度值eps
nargin>3
否
eps<1e-6
是
m=1
t(1,1)=(b-a)*(f(a)+f(b))2
r(1,1)=0
((abs(t(1,m)-r
否
(1,m))2)>eps
是
c=0;m=m+1;j=1
1<=j<=2^(m-1)
否
是
j=j+1
c=c+f(a+(j-0.5)*(b-a)2^(m-1))
r(m,1)=(b-a)*c2^(m-1)
j=2
否
j<=m
是
k=1
j=j+1
k<=m-j+1
否
是
k=k+1 r(k,j)=r(k+1,j-1)+(r(k+1,j-
1)-r(k,j- 1))(4^(j-1)-1)
t(1,j)=r(1,j-1)+2*(4^(j-2)-1)*< br>(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=');
if v0
j=Q*r^22g;

F=f*Q;

beta=F*rj;

a=-F(Qg);

t=(v0-w0*r)(beta*r-a)

v=v0+a*t

3、 执行命令
>> move
r=1
Q=100
g=9.81
f=0.1
v0=3
w0=2
t =
0.3398
v =
2.6667
三、高程
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+300N*m

y0=100+300N*n
3、执行命令
>> high(100)
运行结果如下：
作出拟合曲面为











求得结果：
zh =
Columns 1 through 7
636.0000 639.8173 643.5115 647.0826 650.5308 653.8558 657.0579
639.0542 642.7968 646.4175 649.9164 653.2935 656.5488 659.6822
642.0349 645.7037 649.2519 652.6794 655.9864 659.1727 662.2384
644.9421 648.5381 652.0146 655.3717 658.6094 661.7276 664.7265
647.7759 651.2999 654.7057 657.9932 661.1625 664.2136 667.1464
650.5363 653.9892 657.3251 660.5440 663.6458 666.6305 669.4982
653.2231 656.6060 659.8729 663.0240 666.0591 668.9784 671.7818
655.8365 659.1502 662.3491 665.4332 668.4027 671.2574 673.9974
658.3765 661.6218 664.7536 667.7717 670.6763 673.4673 676.1448
660.8430 664.0209 667.0864 670.0395 672.8801 675.6083 678.2240
663.2360 666.3475 669.3476 672.2365 675.0140 677.6802 680.2352
665.5556 668.6015 671.5372 674.3627 677.0781 679.6832 682.1781
667.8017 670.7829 673.6551 676.4182 679.0722 681.6172 684.0530
669.9743 672.8918 675.7014 678.4029 680.9965 683.4821 685.8597
672.0735 674.9282 677.6760 680.3169 682.8510 685.2781 687.5983
674.0992 676.8920 679.5790 682.1602 684.6355 687.0051 689.2688
676.0514 678.7833 681.4103 683.9326 686.3502 688.6630 690.8711
677.9302 680.6020 683.1700 685.6344 687.9950 690.2520 692.4053
679.7355 682.3481 684.8581 687.2654 689.5700 691.7720 693.8713
681.4674 684.0217 686.4745 688.8256 691.0751 693.2230 695.2692
683.1258 685.6228 688.0192 690.3150 692.5103 694.6049 696.5990
684.7107 687.1513 689.4923 691.7338 693.8756 695.9179 697.8607
686.2222 688.6073 690.8938 693.0817 695.1711 697.1619 699.0542
687.6602 689.9907 692.2236 694.3589 696.3967 698.3369 700.1796
689.0248 691.3016 693.4818 695.5654 697.5524 699.4429 701.2368
690.3159 692.5399 694.6683 696.7011 698.6383 700.4799 702.2259
691.5335 693.7057 695.7832 697.7661 699.6543 701.4479 703.1469
692.6777 694.7989 696.8264 698.7603 700.6004 702.3469 703.9997
693.7484 695.8196 697.7980 699.6837 701.4767 703.1769 704.7844
694.7456 696.7677 698.6979 700.5364 702.2831 703.9379 705.5010
695.6694 697.6433 699.5262 701.3184 703.0196 704.6300 706.1494
696.5197 698.4463 700.2829 702.0295 703.6862 705.2530 706.7297
697.2966 699.1768 700.9679 702.6700 704.2830 705.8070 707.2419
698.0000 699.8347 701.5813 703.2397 704.8099 706.2920 707.6860
698.5943 700.3851 702.0886 703.7048 705.2338 706.6754 708.0298
699.0484 700.7973 702.4598 704.0359 705.5255 706.9287 708.2455
699.3689 701.0779 702.7013 704.2391 705.6913 707.0580 708.3390
699.5625 701.2334 702.8196 704.3209 705.7375 707.0693 708.3164
699.6359 701.2705 702.8211 704.2877 705.6703 706.9689 708.1835
699.5958 701.1956 702.7122 704.1456 705.4958 706.7628 707.9465
699.4488 701.0154 702.4995 703.9012 705.2204 706.4571 707.6114
699.2017 700.7364 702.1893 703.5606 704.8501 706.0580 707.1841
698.8610 700.3651 701.7882 703.1303 704.3914 705.5715 706.6706
698.4335 699.9082 701.3025 702.6165 703.8503 705.0038 706.0770
697.9259 699.3721 700.7387 702.0257 703.2331 704.3610 705.4092
697.3449 698.7635 700.1033 701.3641 702.5461 703.6492 704.6734
696.6970 698.0889 699.4026 700.6381 701.7954 702.8745 703.8754
695.9890 697.3549 698.6432 699.8540 700.9873 702.0430 703.0212
695.2276 696.5680 697.8316 699.0182 700.1280 701.1609 702.1170
694.4195 695.7349 696.9740 698.1370 699.2237 700.2343 701.1686
693.5713 694.8620 696.0771 697.2167 698.2807 699.2692 700.1821
692.6897 693.9559 695.1472 696.2636 697.3051 698.2718 699.1635
691.7814 693.0232 694.1908 695.2841 696.3033 697.2482 698.1189
690.8530 692.0705 693.2143 694.2846 695.2813 696.2045 697.0541
689.9113 691.1042 692.2243 693.2713 694.2455 695.1468 695.9752
688.9630 690.1311 691.2270 692.2507 693.2021 694.0813 694.8883
688.0146 689.1576 690.2291 691.2290 692.1573 693.0141 693.7992
687.0729 688.1903 689.2369 690.2125 691.1173 691.9512 692.7141
686.1445 687.2358 688.2568 689.2077 690.0883 690.8987 691.6390
685.2362 686.3006 687.2954 688.2208 689.0766 689.8629 690.5797
684.3546 685.3913 686.3591 687.2582 688.0884 688.8498 689.5424
683.5064 684.5144 685.4543 686.3262 687.1299 687.8655 688.5331
682.6983 683.6766 684.5875 685.4311 686.2073 686.9162 687.5577
681.9369 682.8843 683.7651 684.5793 685.3269 686.0079 686.6223
681.2289 682.1442 682.9936 683.7772 684.4949 685.1468 685.7328
680.5811 681.4628 682.2794 683.0310 683.7175 684.3389 684.8952
680.0000 680.8466 681.6290 682.3471 683.0009 683.5904 684.1157
679.4545 680.2651 681.0123 681.6959 682.3161 682.8729 683.3661
678.9091 679.6833 680.3948 681.0437 681.6299 682.1535 682.6144
678.3636 679.1011 679.7767 680.3904 680.9423 681.4324 681.8607
677.8182 678.5185 679.1578 679.7361 680.2533 680.7096 681.1048
677.2727 677.9355 678.5382 679.0806 679.5629 679.9850 680.3469
676.7273 677.3522 677.9179 678.4242 678.8711 679.2587 679.5869
676.1818 676.7686 677.2969 677.7666 678.1779 678.5306 678.8248
675.6364 676.1846 676.6751 677.1080 677.4832 677.8008 678.0607
675.0909 675.6002 676.0527 676.4483 676.7872 677.0692 677.2944
674.5455 675.0155 675.4295 675.7876 676.0897 676.3359 676.5261
674.0000 674.4304 674.8056 675.1258 675.3909 675.6009 675.7557
673.4545 673.8449 674.1810 674.4629 674.6906 674.8641 674.9833
672.9091 673.2591 673.5557 673.7990 673.9890 674.1255 674.2087
672.3636 672.6729 672.9297 673.1340 673.2859 673.3852 673.4321
671.8182 672.0864 672.3030 672.4680 672.5814 672.6432 672.6534
671.2727 671.4995 671.6756 671.8009 671.8755 671.8994 671.8726
670.7273 670.9123 671.0474 671.1327 671.1682 671.1539 671.0897
670.1818 670.3246 670.4185 670.4635 670.4595 670.4066 670.3048
669.6364 669.7367 669.7890 669.7932 669.7494 669.6576 669.5177
669.0909 669.1483 669.1587 669.1218 669.0379 668.9068 668.7286
668.5455 668.5597 668.5277 668.4494 668.3250 668.1543 667.9374
668.0000 667.9706 667.8959 667.7759 667.6107 667.4001 667.1442
667.4545 667.3812 667.2635 667.1014 666.8949 666.6441 666.3488
666.9091 666.7914 666.6303 666.4258 666.1778 665.8863 665.5514
666.3636 666.2013 665.9965 665.7491 665.4592 665.1268 664.7519
665.8182 665.6108 665.3619 665.0714 664.7393 664.3656 663.9503
665.2727 665.0200 664.7266 664.3926 664.0179 663.6026 663.1466
664.7273 664.4288 664.0906 663.7127 663.2952 662.8379 662.3409
664.1818 663.8372 663.4539 663.0318 662.5710 662.0714 661.5331
663.6364 663.2453 662.8165 662.3498 661.8454 661.3032 660.7232
663.0909 662.6530 662.1783 661.6668 661.1184 660.5332 659.9112
662.5455 662.0604 661.5395 660.9827 660.3900 659.7615 659.0971
662.0000 661.4674 660.8999 660.2975 659.6602 658.9881 658.2810
????????????
Zh为100*100的矩阵，此处只列出部分值，其余省略。
并解得最高点和最高程为
：
x0 =
166
y0 =
196
ans =
718.1243
即最高点在坐标（166,196）处，且高程为718.1243m。