数学模型与数学建模实验五

1、实验报告五学院名称：理学院 专业年级： 姓 名： 学 号：课 程：数学模型与数学建模 报告日期：2015年12月8日一、实验题目例2.2.1 水库库容量与高程 设一水库将河道分为上、下游两个河段，降雨的开始时刻为8时，这是水位的高程为168m，水库容量为，预测上游的流量，d取值如表2.2.1所示。 表2.2.1 上有流量的预测812162430444856360054007800920010100350025001600已知水库中水的容量与水位高程H（m）的数值关系为表2.2.2 表2.2.2 水库库容量与水位高程的关系23.9324.0624.1224.3324.4724.624.75168

2、.75168.8168.85168.9168.95169169.05 如果当日从8时开始，水一直保持的泄流量，根据所给数据，预报从降雨时刻到56h以内每小时整点时刻水库中水的库容量与水位高程。例2.2.2 地下含沙量 某地区有优质细沙埋在地下，某公司拟在此处采沙，已得到该地区钻探资料图的一角如下表，在每个格点上有三个数字列，都是相对于选定基点的高度（m），最上面的数字是覆盖表面的标高，中间的数字是沙层顶部的标高最下面的数字是沙层底部的标高，每个格子都是正方形，边长50m。画星号处，即沼泽表层地带，没有钻探数据。试估计整个矩形区域内的含沙量。ABCDEFGH022.420.05.8*22.518

3、.40.523.017.80.413.218.00.4123.019.96.023.120.03.223.220.01.623.419.81.023.519.91.124.020.01.024.019.80.824.019.40.9223.119.82.223.219.71.423.419.40.623.420.00.523.520.10.324.220.3-0.224.120.3-0.124.120.50.0二、实验目的 插值模型是数据挖掘的另一类模型，插值（Interpolation）的目的是根据能够获得的观测数据推测缺损的数据，此时观测数据被视为精确的基准数据，寻找一个至少满足条件的函数

4、，使得，在本节我们强调的是插值模型的应用，而不是插值方法的构造。三、问题陈述2.2.1 一维插值 例2.2.1 水库库容量与高程2.2.2 二维插值 例2.2.2 地下含沙量2.2.3 泛克里金插值四、模型及求解结果 2.2.1 一维插值 一元函数差值公式为其中是满足条件的函数，依据插值的公式，如最近邻差值，线性插值、分段三次Hermite插值等，分别取阶梯函数、线性函数、三次多项式函数等，相应的数学表达式可以查阅本科生数值计算教材。下面先通过简单的Matlab一维插值令interpret1了解相应的计算结果。例：2.2.1 水库库容量与高程 为了给出每小时的报告，需要补充每小时整点时刻上有流

5、量的数据，以及相应不同库容量的水位高程。 假设（a）已知数据准确（b）相邻两个时刻之间的流量变化是现行的（c）相邻两个水位高程之间的高程对水的库容量的变化也是线性的 首先，利用Matlab线性插值令，确定每小时的上游流量q（t）.由程序在Matlab中运行的结果为：q = 1.0e+004 * Columns 1 through 9 0.3600 0.4050 0.4500 0.4950 0.5400 0.6000 0.6600 0.7200 0.7800 Columns 10 through 18 0.7975 0.8150 0.8325 0.8500 0.8675 0.8850 0.902

6、5 0.9200 0.9350 Columns 19 through 27 0.9500 0.9650 0.9800 0.9950 1.0100 0.9629 0.9157 0.8686 0.8214 Columns 28 through 36 0.7743 0.7271 0.6800 0.6329 0.5857 0.5386 0.4914 0.4443 0.3971 Columns 37 through 45 0.3500 0.3000 0.2500 0.2410 0.2320 0.2230 0.2140 0.2050 0.1960 Columns 46 through 49 0.1870

7、0.1780 0.1690 0.1600 然后确定每个时刻t的水库容量，因为，水库容量=原库存量+流入量-泄流量，即：这里我们遇到数值积分，被积函数没有解析表达式，只有一个数列表示，表示在i整点时刻的流量，利用Matlab逐点积分指令，可以得到水库容量在每一刻的值。 最后确定每时刻t水库的水位高程h（t），因为最大水库容量已经远远超出了已知数据范围，需要利用外插方法补充数据，确定水库高程对水库容量的依赖关系h=H（v）。最后利用函数复合得到水位高程确定每小时的上游流量q（t）利用函数复合得到水位高程2.2.2 二维插值二维插值大致可以分为两种，规则点插值和散乱点插值，前者即利用在方形网格点给定

8、的数据，在加密网格点上补充相应函数值，插值公式为其中是满足的二元函数，如双线性函数、双三次多项式函数、样条函数等，一句不同的插值方式选取，相应的表达式可以从本科生数值计算教材中查阅。这里先通过Matlab二维插值指令interp2了解相应的计算结果散乱点插值是要利用在不规则排列的观测点上的数据确定其他点上的函数值，插值公式形式为常用到的插值方法，如距离加权反比插值，Kriging插值，需要查阅专业书籍，下面先通过Matlab二维插值指令griddata了解相应的计算结果例2.2.2 地下含沙量 假设地下沙层连续变化，于是可以利用积分运算矩形区域内的含沙量。首先通过插值补充缺失的数据，采用一维样

9、条插值。为了提高梯形积分精度，利用二维样条插值加密数据，得到边长为0.5m的长方形网格上的数据。最后，将二重积分化为累次积分，利用梯形公式得到积分的近似值。v = 6.4290e+005如果在一维插值补齐缺失数据中用线性插值代替上面的样条插值，将得到沙层体积634525，和以上结果有些差别。 对于加密二维网格，应该加到多密？注意，从理论上讲，网格越密计算精度越高，但是从计算角度看，网格越密计算误差累计也越大，所以需要通过逐步加密网格，确定精度最优的积分值。 如果直接利用二维插值，例如采用Matlab的散乱点插值指令，仍用二次梯形公式，会得到相近的结果。v = 6.3572e+0052.2.3

10、泛克里金插值考虑到沙子是一种特殊的地址，下面试用地质学常用的泛克里金插值方法计算。先介绍克里金插值方法，这本身就是运用统计学方法建模的一个有趣的例子。1951年南非地质学家Krige讲处矿藏的贮藏量看成是随机函数的一个实现，提出了依据观测值，寻找基函数，以获得的具有形式 的最小方差无偏估计，取的条件期望 给出未知点估值的插值方法，即考虑形如 的随机函数，其中在地质学中称为飘移，在一些随机分析的场合也成为趋势，是均值为的未知的随机变量，是多项式基函数，例如这里取，是满足，的随机函数是给定的且满足当当的核函数。这里按通常取法，取高斯核函数根据无偏估计要求， 即： 从而得到无偏条件，对任意的有： 在

11、这个约束条件下极小化方差 =：其中，。引入拉格朗日乘子，由拉格朗日函数的极小值点必为其驻点的结论，得到关系式 其中，于是，在未知点上的随机函数值可以用它的最小方差无偏估计的条件期望近似，从而得到泛克里金插值函数： = =其中练习：证明上面构造的泛克里金函数是连续的，且通过已知点，即插值函数满足： 用泛克里金插值得到的边长为0.5m的加密网格上的沙层厚度，通过二次梯形公式，得到沙层体积近似值为637320五、程序代码2.2.1 一维插值>> x=0:4:20;>> y=37 51 45 74 83 88;>> xx=0:1:20;>> y1=int

12、erp1(x,y,xx,'nearest'); %最近邻插值，间断函数>> y2=interp1(x,y,xx); %线性插值，连续函数>> y3=interp1(x,y,xx,'cubic'); %分段三次Hermite插值，一阶连续函数>> y4=interp1(x,y,xx,'splinet'); %样条插值，二阶倒数连续的函数>> subplot(2,2,1),plot(x,y,'kd',xx,y1),title('nearest')>>

13、 subplot(2,2,2),plot(x,y,'kd',xx,y2),title('linear')>> subplot(2,2,3),plot(x,y,'kd',xx,y3),title('cubic')>> subplot(2,2,4),plot(x,y,'kd',xx,y4),title('spline')例2.2.1 水库库容量与高程利用Matlab线性插值令，确定每小时的上游流量q（t）>> T=8,12,16,24,3

14、0,44,46,56;>> Q=3600,5400,7800,9200,10100,3500,2500,1600;>> t=8:56;>> q=interp1(T,Q,t,'linear')%得到每小时上游流量>> plot(T,Q,'kd',t,q),title('time-flow')水库容量在每一刻的值>> v=21.9+36*10(-6)*cumtrapz(t,q-103*ones(size(t); %每小时水库容量>> vmax=max(v); %得到最大的水库容量

15、为30.7524（108立方米）利用函数复合得到水位高程>> V=21.9 23.93 24.06 24.12 24.33 24.47 24.3 24.75;>> H=168 168.75 168.8 168.85 168.9 168.95 169 169.05;>> h=interp1(V,H,v,'linear','extrap'); %得到每小时水位高程>> hmax=max(h);%最大水位高程171.0508米>> plot(V,H,'kd',v,h),title('v

16、olume-atitude')2.2.2 二维插值二维插值指令>> x=0:4:16;y=0:4:16;>> z=620 730 800 850 870;760 880 970 720 1050880 1080 630 1250 1280;980 1180 1320 1450 14201060 1230 1390 1500 1500;>> X,Y=meshgrid(0:16,0:16);>> Z1=interp2(x,y,z,X,Y,'nearest');Z2=interp2(x,y,z,X,Y);>> Z3=

17、interp2(x,y,z,X,Y,'cubic');Z4=interp2(x,y,z,X,Y,'spline');>> subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title

18、('spline')散乱点插值>> x=1 2 3 4 5;>> y=4 1 3 2 5;>> z=4 1 6 2 5;>> X,Y=meshgrid(0:0.2:5);>> Z1=griddata(x,y,z,X,Y,'nearest');>> Z2=griddata(x,y,z,X,Y);>> Z3=griddata(x,y,z,X,Y,'cubic');>> Z4=griddata(x,y,z,X,Y,'v4');>>

19、 subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title('v4')例2.2.2 地下含沙量通过一维插值补充缺失的数据>> D=23.0,23.1,23.2,23.4,23.5,24.0,24

20、.0,24.0;23.1,23.3,23.4,23.4,23.5,24.2,24.1,24.1;>> E=19.9,20.0,20.0,19.8,19.9,20.0,19.8,19.6;19.8,19.7,19.4,20.0,20.1,20.3,20.3,20.5;>> F=6.0,3.2,1.6,1.0,1.1,1.0,0.8,0.9;2.2,1.4,0.6,0.5,0.3,-0.2,-0.1,0.0;>> P=1,6,7,8;K=1:8;>> A=22.4,22.5,23.0,23.2;>> A1=interp1(P,A,K,&#

21、39;spline');a=A1;D>> B=20.0,18.4,17.8,18.0;>> B1=interp1(P,B,K,'spline');b=B1;E>> C=5.8,0.5,0.4,0.4;>> C1=interp1(P,C,K,'spline');c=C1;F二维样条插值加密数据，得到边长为0.5m的长方形网格上的数据>> M,N=meshgrid(1:8,0:2);>> delta=0.01;X,Y=meshgrid(1:delta:8,0:delta:2);>&

22、gt; a1=interp2(M,N,a,X,Y,'spline');>> b1=interp2(M,N,b,X,Y,'spline');>> c1=interp2(M,N,c,X,Y,'spline');>> mesh(X,Y,a1),hold on,>> mesh(X,Y,b1),hold on,mesh(X,Y,c1)最后，将二重积分化为累次积分，利用梯形公式得到积分的近似值。>> f=b1-c1;>> v=trapz(trapz(f,1),2)*delta2*2500

23、v = 6.4290e+005（v为沙层体积） 如果直接利用二维插值，例如采用Matlab的散乱点插值指令，仍用二次梯形公式，会得到相近的结果。>> x0=1 6 7 8 1:8 1:8;>> y0=zeros(1,4) ones(1,8) 2*ones(1,8);>> f0=B E'-C F'>> f=griddata(x0,y0,f0,X,Y,'cubic');>>>> v=trapz(trapz(f,1),2)*delta2*2500v = 6.3572e+005（v为沙层体积）2.2

24、.3 泛克里金插值泛克里金插值的MATLAB指令为：function f=krige2(x0,y0,f0,X,Y)%x0,y0,f0·Ö±ðÊÇÒÑÖªÊý¾ÝµÄºá¡¢×Ý×ø±êºÍ¹Û²âÖµ£¬±í´

25、ï³ÉÐÐÏòÁ¿%X,Y²åÖµµãµÄºá¡¢×Ý×ø±ê£¬¿ÉÒÔÊÇÐÐÏòÁ¿£¬Ò²¿ÉÒÔ&#

26、202;ǾØÕó±í´ïÍø¸ñµã×ø±ê¡£%Êä³öfÊǶÔÓ¦X,YµÄº¯ÊýÖµn=length(x0);syms x y;%¶¨Òå·ûºÅ±äÁ¿h0=(x-x0').2+(y-y0').2;%¾àÀ뺯