数学建模与数学实验：矿区储藏量和面积的计算问题研究

数学建模与数学实验：矿区储藏量和面积的计算
问题研究
研究目标
本实验的目的是通过对矿区面积的计算,掌握定积分的近似计算方法,对有关数值积分的有关理论和数值计算方法有所了解。

解决问题
1.计算积分4
2
() f x dx
的近似值。

2.矿区储量
问题1：计算积分4
2()f x dx ⎰的近似值。

已知函数()y f x =的一些数据点如下:
分别用矩形,梯形和辛普生公式计算积分4
2()f x dx ⎰的近似值。

[问题分析]
这个问题就是基本的计算,我们可以直接套用公式进行编程计算即可。

复合矩形求积公式,分为三种情况:
111
11111(1) ()()()
(2) ()()()
(3) ()()()2n b i i i a i n b i i i a i n b i i
i i a i f x dx f x x x f x dx f x x x x x f x dx f x x --=-=--=⎧
=-⎪⎪
⎪
=-⎨⎪
⎪+=-⎪⎩
∑⎰∑⎰∑⎰ 梯形求积公式: ()[()()]2
b
a a b
f x dx f a f b +=
+⎰ 辛普生求积公式: ()[()()()]62
b
a b a a b
f x dx f a f f b -+=++⎰
[实验程序]
⏹ function 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
问题2：矩形矿区储藏量
煤矿的储量估计,下表给出了某露天煤矿在平面矩形区域(800m ⨯600m)上,在纵横均匀的网格交点处测得的煤层厚度(单位:m)(由于客观原因,有些点无法测量煤层厚度,这里用/标出),其中的每个网格都为(10m ⨯8m)的小矩形,试根据这些数据,来估算出该矩形区域煤矿的储藏量(体积)。

[问题分析]
这个问题时在计算矿区面积的基础上,扩展到计算矿区的体积。

按照求体积的方法,我们可以知道,整个矿区的体积应为: 800
600
00(,)V dx z x y dy =⎰⎰,其中(,)z x y 表示的是在(,)x y 点处煤层的高度。

但是题目中只给出了51个网格点处的高度,所以我们可以采用书上介绍的求积分的近似值的方法来进行求解。

[模型建立]
(1)首先,由题目可知,矿区的总面积为: 280060048000m m m ⨯=,而每个小网格的面积为:
210880m m m ⨯=,所以整个矿区被分成了
480000
600080
=个小的网格。

我们可以画出整个矿区的分割图,如图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,'+')
[运行结果]
(2)我们借用书上求面积的思想,对体积同样适用,将体积800600
00(,)
V dx z x y dy
=⎰⎰进行分
割后,求积分和,即: 800
600
480000
0080
(,)()V dx z x y dy z s ds ==⎰⎰⎰,将三维问题转化为二维问题,然
后借用书上的二维求积分的计算方法进行求解。

(3)要估计6000个区域当中各个区域的高度值,我们考虑将51个已知数据进行编号,并记高度值为()f s ,设未知变量s 为区域的面积,则[80,480000]s ∈,将已知的51个数据列表如下:
(4)对上表中的数据按照书上的求积分的思想进行处理。

480000
80()V f s ds =⎰
方法一：采用样条函数积分法
⏹function shiyan1421
⏹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 = 6
即为整个矿区的体积。

5.982010
方法二：采用分段插值法计算近似解
总共是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表示的是整个区间上的体积和.
从上表中的计算结果我们可以发现用拉格朗日积分的结果误差是比较大的,所以我们采用后两种积分,但是总体的积分还是不如分段积分值准确,所以我们采用分段积分,最后将后两种分段积分的值进行平均后得到,整个矿区的体积为:56725803m,和方法一的计算结果相差不多.
积分结果画图后得到:。