实验目的:
本实验的目的是通过对矿区面积的计算,掌握定积分的近似计算方法,对有关数值积分的有关理论和数值计算方法有所了解。
问题一:
问题重述:
已知函数y?f(x)的一些数据点如下:
分别用矩形,梯形和辛普生公式计算积分?f(x)dx的近似值。
24问题分析:
这个问题就是基本的计算,我们可以直接套用公式进行编程计算即可。
nb??(1) ?af(x)dx??f(xi?1)(xi?xi?1)i?1?nb?复合矩形求积公式,分为三种情况: ?(2) ?f(x)dx??f(xi)(xi?xi?1)
ai?1?n?bxi?1?xi)(xi?xi?1)?(3) ?af(x)dx??f(2i?1?梯形求积公式:
?baf(x)dx?ba辛普生求积公式: 实验程序:
?a?b[f(a)?f(b)] 2b?aa?bf(x)dx?[f(a)?f()?f(b)]
62function shiyan131
x=[2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0];
y=[1.65,1.56,1.38,1.12,0.77,0.34,-0.15,-0.7,-1.3,-1.91,-2.01]; n=length(x) for i=2:n
s1(i-1)=y(i-1)*(x(i)-x(i-1)); s2(i-1)=y(i)*(x(i)-x(i-1)); end
s11=sum(s1) s12=sum(s2) for i=2:(n-1)
s3(i-1)=y(i)*(x(i+1)-x(i-1)); end
s13=sum(s3)
s4=(x(n)-x(1))*(y(n)+y(1))/2
s5=(x(n)-x(1))*(y(1)+4*y((n+1)/2)+y(n))/6
运行结果:
复合矩形求积法: 方法一: s11 = 0.5520
方法二: s12 = -0.1800 方法三: s13 = 0.4440
梯形求积法: s4 = -0.3600 辛普生求积法: s5 = 0.3333
问题二: 问题重述:
煤矿的储量估计,下表给出了某露天煤矿在平面矩形区域(800m?600m)上,在纵横均匀的网格交点处测得的煤层厚度(单位:m)(由于客观原因,有些点无法测量煤层厚度,这里用/标出),其中的每个网格都为(10m?8m)的小矩形,试根据这些数据,来估算出该矩形区域煤矿的储藏量(体积) A B C D E F G H I J K 1 / / 12.5 13.5 17.2 / 8.8 14.7 8.0 13.0 / 2 / / / 15.6 18.2 13 6.4 8.9 9.2 11.7 / 3 / 12 13.5 13.5 17.8 16.9 13.2 / / / / 4 7.5 12.6 14.9 18.7 17.7 17.5 14.7 13 / / 6.5 5 8.9 7.8 12.4 13.5 15.7 17.6 11.7 9.6 9.2 9.5 8.6 6 / / / 13.7 13.6 16.5 12.5 8.7 9.7 / / 7 / / 8.6 11.8 12.5 11.3 13.4 / / / /
模型的建立:
(1)首先,由题目可知,矿区的总面积为: 800m?600m?48000m2,而每个小网格的面积为: 10m?8m?80m2,所以整个矿区被分成了格。我们可以画出整个矿区的分割图,如图1: 实验程序:
function shiyan141 for i=1:75 for j=1:80
x((i-1)*80+j)=8*i; y((i-1)*80+j)=10*j; end end
plot(x,y,'+')
480000?6000个小的网80运行结果:
图1
(2)我们借用书上求面积的思想,对体积同样适用,将体积V??进行分割后,求积分和,即: V??80008000dx?6000z(x,y)dydx?6000z(x,y)dy??48000080z(s)ds,将三维问题
转化为二维问题,然后借用书上的二维求积分的计算方法进行求解。
(3)要估计6000个区域当中各个区域的高度值,我们考虑将51个已知数据进行编号,并记高度值为f(s),设未知变量s为区域的面积,则s?[80,480000],将已知的51个数据列表如下: s 80 80*120 f (s) 12.5 13.5 s 80*840 80*960 f (s) 15.6 18.2 s 80*1680 80*1800 f (s) 12 13.5 s 80*2520 80*2640 f (s) 12.6 14.9 s 80*3360 80*3480 f (s) 6.5 8.9 s 80*4200 80*4320 f (s) 11.7 9.6 s 80*5040 80*5160 f (s) 16.5 12.5 s 80*5580 80*6000 f (s) 11.3 13.4 80*240 17.2 80*1080 13 80*1920 13.5 80*2760 18.7 80*3600 7.8 80*4440 9.2 80*5280 8.7 80*360 8.8 80*1200 6.4 80*2040 17.8 80*2880 17.7 80*3720 12.4 80*4560 9.5 80*5400 9.7 80*480 14.7 80*1320 8.9 80*2160 16.9 80*3000 17.5 80*3840 13.5 80*4680 8.6 80*5520 8.6 80*600 8.0 80*1440 9.2 80*2280 13.2 80*3120 14.7 80*3960 15.7 80*4800 13.7 80*5640 11.8 4800008080*720 13.0 80*1560 11.7 80*2400 7.5 80*3240 13 80*4080 17.6 80*4920 13.6 80*5760 12.5 (4) 对上表中的数据按照书上的求积分的思想进行处理。V??? 方法一: 采用样条函数积分法,实验程序如下:
function shiyan1421
f(s)ds
x(1)=80; for i=2:51
x(i)=80*120*(i-1); end
y=[12.5,13.5,17.2,8.8,14.7,8.0,13.0,15.6,18.2,13,6.4,8.9,9.2,11.7,12,13.5,13.5,17.8,16.9,13.2,
7.5,12.6,14.9,18.7,17.7,17.5,14.7,13,6.5,8.9,7.8,12.4,13.5,15.7,17.6,11.7,9.6,9.2,9.5,8.6, 13.7,13.6,16.5,12.5,8.7,9.7,8.6,11.8,12.5,11.3,13.4]; ss=spline(x,y); js=fnint(ss);
sd=ppval(js,[80,480000])
运行结果为:
sd = 5.9820?106 即为整个矿区的体积。
? 方法二: 采用分段插值法计算近似解,总共是51组数据,恰好分成三组,每组
17个,进行插值求解。s的区间分为: [80,80*120*16], [80*120*17, 80*120*33], [80*120*34, 80*120*50],我们也采用书上的方法,每组分别用三种不同的插值方法估计f (s)的形式,再进行积分求解,矿区的总的体积就是这三段积分的和。
实验程序:
function shiyan1422 x(1)=80; for i=2:51
x(i)=80*120*(i-1); end
y=[12.5,13.5,17.2,8.8,14.7,8.0,13.0,15.6,18.2,13,6.4,8.9,9.2,11.7,12,13.5,13.5,
17.8,16.9,13.2,7.5,12.6,14.9,18.7,17.7,17.5,14.7,13,
6.5,8.9,7.8,12.4,13.5,15.7,17.6,11.7,9.6,9.2,9.5,8.6,13.7,13.6, 16.5,12.5,8.7,9.7,8.6,11.8,12.5,11.3,13.4]; %第一区间: s1(1)=80; for i=2:17
s1(i)=80*120*(i-1); end
y1=[12.5,13.5,17.2,8.8,14.7,8,13,15.6,18.2,13,6.4,8.9,9.2,11.7,12,13.5,13.5]; t1=80:80:80*120*16; lglr1=lglrcz(s1,y1,t1); lglrjf1=80*trapz(lglr1) fdxx1=interp1(s1,y1,t1); fdxxjf1=80*trapz(fdxx1)
scyt1=interp1(s1,y1,t1,'spline'); scytjf1=80*trapz(scyt1) %第二区间: for i=1:17
s2(i)=x(i+17); y2(i)=y(i+17); end
t2=80*120*17:80:80*120*33; lglr2=lglrcz(s2,y2,t2); lglrjf2=80*trapz(lglr2) fdxx2=interp1(s2,y2,t2); fdxxjf2=80*trapz(fdxx2)
scyt2=interp1(s2,y2,t2,'spline'); scytjf2=80*trapz(scyt2) %第三区间: for i=1:17 s3(i)=x(i+34); y3(i)=y(i+34); end
t3=80*120*34:80:80*120*50; lglr3=lglrcz(s3,y3,t3); lglrjf3=80*trapz(lglr3) fdxx3=interp1(s3,y3,t3); fdxxjf3=80*trapz(fdxx3)
scyt3=interp1(s3,y3,t3,'spline'); scytjf3=80*trapz(scyt3) %整个区间 t(1)=80; for i=2:51
t(i)=80*120*(i-1); end
lglr=lglrcz(x,y,t); lglrjf=80*trapz(lglr) fdxx=interp1(x,y,t); fdxxjf=80*trapz(fdxx)
scyt=interp1(x,y,t,'spline'); scytjf=80*trapz(scyt) %三个区间的总和
lglrjfz=lglrjf1+lglrjf2+lglrjf3 fdxxjfz=fdxxjf1+fdxxjf2+fdxxjf3 scytjfz=scytjf1+scytjf2+scytjf3 %进行画图
subplot(2,2,1),plot(s1,y1,'*',t1,lglr1,'r',t1,fdxx1,'g',t1,scyt1,'b'), title('区间1')
subplot(2,2,2),plot(s2,y2,'*',t2,lglr2,'r',t2,fdxx2,'g',t2,scyt2,'b'), title('区间2')
subplot(2,2,3),plot(s3,y3,'*',t3,lglr3,'r',t3,fdxx3,'g',t3,scyt3,'b'), title('区间3')
subplot(2,2,4),plot(x,y,'*',t,lglr,'r',t,fdxx,'g',t,scyt,'b'),title('总区间')
运行结果:
其中总矿区体积1表示的是三个区间的体积和,而总矿区体积2表示的是整个区间上的体积和. 方式 区间1体积 区间2体积 区间3体积 总矿区体积1 总矿区体积2 拉格朗日-1489800 -4854900 -2888000 -9232800 4996400 插值积分 分段线性 1887300 2040500 1756800 5684560 4996400 插值积分 三次样条1881100 2034700 1745300 5661100 4996400 插值积分 从上表中的计算结果我们可以发现用拉格朗日积分的结果误差是比较大的,所以我们采用后两种积分,但是总体的积分还是不如分段积分值准确,所以我们采用分段积分,最后将后两种分段积分的值进行平均后得到,整个矿区的体积为:5672580m3,和方法一的计算结果相差不多. ? 积分结果画图后得到: