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

实验报告五学院名称：理学院专业年级：姓名：学号：课程：数学模型与数学建模报告日期：2015年12月8日一、实验题目例2.2.1水库库容量与高程设一水库将河道分为上、下游两个河段，降雨的开始时刻为8时，这是水位的高程为168m，水库容量为38109.21m，预测上游的流量smtQ/3，d取值如表2.2.1所示。表2.2.1上有流量tQ的预测ht/81216243044485613/smQ360054007800920010100350025001600已知水库中水的容量3810mV与水位高程H（m）的数值关系为表2.2.2表2.2.2水库库容量与水位高程的关系3810/mV23.9324.0624.1224.3324.4724.624.75mH/168.75168.8168.85168.9168.95169169.05如果当日从8时开始，水一直保持sm/10003的泄流量，根据所给数据，预报从降雨时刻到56h以内每小时整点时刻水库中水的库容量与水位高程。例2.2.2地下含沙量某地区有优质细沙埋在地下，某公司拟在此处采沙，已得到该地区钻探资料图的一角如下表，在每个格点上有三个数字列，都是相对于选定基点的高度（m），最上面的数字是覆盖表面的标高，中间的数字是沙层顶部的标高最下面的数字是沙层底部的标高，每个格子都是正方形，边长50m。画星号处，即沼泽表层地带，没有钻探数据。试估计整个矩形区域内的含沙量。ABCDEFGH022.420.05.8************22.518.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）的目的是根据能够获得的观测数据推测缺损的数据，此时观测数据niiiyx1,被视为精确的基准数据，寻找一个至少满足条件的函数xyy，使得nixyyii,,2,1,，在本节我们强调的是插值模型的应用，而不是插值方法的构造。三、问题陈述2.2.1一维插值例2.2.1水库库容量与高程2.2.2二维插值例2.2.2地下含沙量2.2.3泛克里金插值四、模型及求解结果2.2.1一维插值一元函数差值公式为niiixyxy1其中xi是满足条件ijix的函数，依据插值的公式，如最近邻差值，线性插值、分段三次Hermite插值等，分别取阶梯函数、线性函数、三次多项式函数等，相应的数学表达式可以查阅本科生数值计算教材。下面先通过简单的Matlab一维插值令interpret1了解相应的计算结果。例：2.2.1水库库容量与高程为了给出每小时的报告，需要补充每小时整点时刻上有流量的数据，以及相应不同库容量的水位高程。假设（a）已知数据准确（b）相邻两个时刻之间的流量变化是现行的（c）相邻两个水位高程之间的高程对水的库容量的变化也是线性的首先，利用Matlab线性插值令，确定每小时的上游流量q（t）.由程序在Matlab中运行的结果为：q=1.0e+004*Columns1through90.36000.40500.45000.49500.54000.60000.66000.72000.7800Columns10through180.79750.81500.83250.85000.86750.88500.90250.92000.9350Columns19through270.95000.96500.98000.99501.01000.96290.91570.86860.8214Columns28through360.77430.72710.68000.63290.58570.53860.49140.44430.3971Columns37through450.35000.30000.25000.24100.23200.22300.21400.20500.1960Columns46through490.18700.17800.16900.1600然后确定每个时刻t的水库容量tv，因为，水库容量=原库存量+流入量-泄流量sm/1038，即：81036001036009.21588tdssqtvt这里我们遇到数值积分，被积函数sq没有解析表达式，只有一个数列表示，iq表示在i整点时刻的流量，利用Matlab逐点积分指令yxcumptrapz,，可以得到水库容量tvv在每一刻56,8t的值。最后确定每时刻t水库的水位高程h（t），因为最大水库容量已经远远超出了已知数据范围，需要利用外插方法补充数据，确定水库高程对水库容量的依赖关系h=H（v）。最后利用函数复合得到水位高程)(tvHth确定每小时的上游流量q（t）利用函数复合得到水位高程)(tvHth2.2.2二维插值二维插值大致可以分为两种，规则点插值和散乱点插值，前者即利用在方形网格点),(iiyx给定的数据mjnizij,,2,1,,,2,1,，在加密网格点上补充相应函数值，插值公式为nimjijijyxlzyxz11),(),(其中ijl是满足jsikskijyxl),(的二元函数，如双线性函数、双三次多项式函数、样条函数等，一句不同的插值方式选取，相应的表达式可以从本科生数值计算教材中查阅。这里先通过Matlab二维插值指令interp2了解相应的计算结果散乱点插值是要利用在不规则排列的观测点),(jjyx上的数据njzj,,2,1,确定其他点上的函数值，插值公式形式为ijjjinjjjyxlyxlzyxz),(,),(),(1常用到的插值方法，如距离加权反比插值，Kriging插值，需要查阅专业书籍，下面先通过Matlab二维插值指令griddata了解相应的计算结果例2.2.2地下含沙量假设地下沙层连续变化，于是可以利用积分运算矩形区域内的含沙量。首先通过插值补充缺失的数据，采用一维样条插值。为了提高梯形积分精度，利用二维样条插值加密数据，得到边长为0.5m的长方形网格上的数据。最后，将二重积分化为累次积分，利用梯形公式得到积分的近似值。v=6.4290e+005如果在一维插值补齐缺失数据中用线性插值代替上面的样条插值，将得到沙层体积6345253m，和以上结果有些差别。对于加密二维网格，应该加到多密？注意，从理论上讲，网格越密计算精度越高，但是从计算角度看，网格越密计算误差累计也越大，所以需要通过逐步加密网格，确定精度最优的积分值。如果直接利用二维插值，例如采用Matlab的散乱点插值指令，仍用二次梯形公式，会得到相近的结果。v=6.3572e+0052.2.3泛克里金插值考虑到沙子是一种特殊的地址，下面试用地质学常用的泛克里金插值方法计算。先介绍克里金插值方法，这本身就是运用统计学方法建模的一个有趣的例子。1951年南非地质学家Krige讲221,Rxxx处矿藏的贮藏量xf看成是随机函数xF的一个实现，提出了依据观测值njyxjj...,3,2,1,,，寻找基函数njxj...,2,1,，以获得xF的具有形式jjjxFxxF*的最小方差无偏估计，取xF*的条件期望njfxFxFExfjj....,2.1,|**给出未知点x估值的插值方法，即考虑形如xRxPxRxMxFK1的随机函数，其中xM在地质学中称为飘移，在一些随机分析的场合也成为趋势，P是均值为pPE)(的未知的随机变量，}{是多项式基函数，例如这里取1)(1x，12)(xx，23)(xx，214)(xx，225)(xx，216)(xxx，xR是满足0))((xRE，)())()((hhxRxRE的随机函数)(h是给定的且满足当;0)(hh时，有当1)(0hh时，有的核函数。这里按通常取法，取高斯核函数2exp)(hh根据无偏估计要求，1*jjjxEFxxEFxEF即：KnjJjKXPxxP111从而得到无偏条件，对任意的有：jjjxxx在这个约束条件下极小化方差jjjxFxxFE2)(=jjjjjxPxFxxxPxFE2))()((=jjjxRxxRE2)(=jjjjkkjkjxRExRxRExRxREx22=：)0()()(2)()(xDxxAxTT其中nnkjxxA))((，1))(()(njxxxD，1))(()(njxx，1)0(。引入拉格朗日乘子K...,21，由拉格朗日函数KnjjjTTxxxxDxxAxL11)()(212：，的极小值点必为其驻点的结论，得到关系式xxDxAT0其中nnkjxxA))((，,))(()(Knjxx1))(()(njxx，,)(1K1))(()(njxxxD，,))(()(1Kxx于是，在未知点x上的随机函数值xF可以用它的最小方差无偏估计的条件期望njfxFxFExfjj....,2.1,|**近似，从而得到泛克里金插值函数：jTTTjjfxfxfxxf0)()(*=00~)()(1fAxxDTTT=gxxDTT)()(其中1,1)0()0(~,KKKnjOOff练习：证明上面构造的泛克里金函数是连续的，且通过已知点，即插值函数满足：njfxfjj,...,2,1,)(用泛克里金插值得到的边长为0.5m的加密网格上的沙层厚度，通过二次梯形公式，得到沙层体积近似值为6373203m五、程序代码2.2.1一维插值x=0:4:20;y=[375145748388];xx=0:1:20;y1=interp1(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')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,30,44,46,56];Q=[3600,5400