数学建模实验 五DOC.docx
《数学建模实验 五DOC.docx》由会员分享，可在线阅读，更多相关《数学建模实验 五DOC.docx（9页珍藏版）》请在冰豆网上搜索。
数学建模实验五DOC
西北农林科技大学实验报告
学院名称：
理学院专业年级：
2013级信计1班
姓名：
学号：
课程：
数学模型与数学建模报告日期：
2015年12月8日
实验五插值模型与误差估计
实验目的
插值模型是数据挖掘的另一类模型，插值的目的是根据能够获得的观测数据推测缺损的信息，此时观测数据
被视为精确地基准数据，寻找一个至少满足连续条件的函数
，使得
，i=1,2,...,n。
在本节我们强调的是插值模型的应用，而不是插值方法的构造。
一、一维插值
一元函数插值公式为
其中
是满足
的函數，依据插值的方式，如最近邻插值、线性插值、分段三次Hermite插值等，分别取阶梯函数、线性函数、三次多项式函数等，相应的数学表达可以查阅教材。
下面先通过简单的MATLAB一维插值指令了解。
>> x=0:
4:
20;
>> y=[37 51 45 74 83 88];
>> xx=0:
1:
20;
>> y1=interp1（x,y,xx,'nearest'）;
>> y2=interp1（x,y,xx）;
>> y3=interp1（x,y,xx,'cubic'）;
>> y4=interp1（x,y,xx,'spline'）;
>> plot（y1,'DisplayName','y1','YDataSource','y1'）;figure（gcf）
>> plot（y2,'DisplayName','y2','YDataSource','y2'）;figure（gcf）
>> plot（y3,'DisplayName','y3','YDataSource','y3'）;figure（gcf）
>> plot（y4,'DisplayName','y4','YDataSource','y4'）;figure（gcf）
例一水库库容量与高程
1实验题目
设一水库将河道分为上、下游两个河段，降雨的开始时刻为8时，这时水位的高程为168m，水库容量为21.9×10^8m³.预测上游流量Q（t）（m³/s）d取值如表2.2.1所示
表2.2.1上游流量Q（t）的预测
t/h
8
12
16
24
30
44
48
56
Q/m³s-1
3600
5400
7800
9200
10100
3500
2500
1600
已知水库中水的库容量V（108m3）与水位高程H（m）的数值关系为表2.2.2
表2.2.2水库库容量与水位高程关系
V/108m3
23.93
24.06
24.12
24.33
24.47
24.6
24.75
H/m
168.75
168.8
168.85
168.9
168.95
169
169.05
如果从当日8时起,水库一直保持每秒1000立方米的泄流量请按所给数据1000立方米的泄流量,请按所给数据,预报当日20时水库的库容量与水位高程.
2实验内容
解为了给出每小时的预报，需要补充每小时整点时刻上游流量的数据，以及相应库容量的水位高程数据
假设:
（a）已知数据准确
（b）相邻两个时刻之间的流量变化是线性的.
（c）相邻两个水位高程之间的高程对水的库容量的变化也是线性的.
首先，利用MATLAB线性插值指令，确定每小时的上游流量q（t）
>> T=[8,12,16,24,30,44,48,56];
>> Q=[3600,5400,7800,9200,10100,3500,2500,1600];
>> t=8:
56;
>> q=interp1（T,Q,t,'linear'）
>> plot（q,'DisplayName','q','YDataSource','q'）;figure（gcf）
然后确定每个时刻t的水库容量v（t）.因为，水库容量=原库存量+流入量-泄流量（108m3/s）,即
这里我们遇到数值积分，被积函数q（s）没有解析表达式，只是由一个数列表示，qi表示在i时刻的流量。
利用MATLAB得到水库容量v在每一时刻t∈[8,56]的值
>> v=21.9+36*10^（-6）*cumtrapz（t,q-10^3*ones（size（t）））;
>> plot（v,'DisplayName','v','YDataSource','v'）;figure（gcf）
>> vmax=max（v）;
最后确定每个时刻t水库的水位高程h（t）.因为最大的水库容量已经远远超出已知数据范围。
需要利用外插方法补充数据，确定水库高程对水库容量的依赖关系h=H（v），最后利用函数复合得到水位高程h（t）=H（v（t））
>> V=[21.9 23.93 24.06 24.12 24.33 24.47 24.6 24.75];
>> H=[168 168.75 168.8 168.85 168.9 168.95 169 169.05];
>> h=interp1（V,H,v,'linear','extrap'）;
>> plot（h,'DisplayName','h','YDataSource','h'）;figure（gcf）
>> hmax=max（h）;
2、二维插值
二维函数插值大致可以分为两种，规则点插值和散乱点插值，前者即利用在方形网格点
给定的数据
，i=1,2,...,n,j=1,2,...,m,在加密网格点上补充相应函数值，插值公式为
其中
是满足
的二元函数。
先通过MATLAB二维插值指令了解，指令如下：
>> x=0:
4:
16;
y=0:
4:
16;
>> z=[620 730 800 850 870;760 880 970 720 1050;880 1080 630 1250 1280;980 1180 1320 1450 1420;1060 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=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（'spline'）
散乱点插值是要利用在不规则的观测点
上的数据zj,j=1,2,...,n,确定其他点上的函数值，插值公式形式为
，
常用到的插值方法。
>> 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'）;
>> 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'）
例二、地下含沙量
1.实验题目
某地优质细沙埋入地下，某公司拟在此处采沙。
已得到该地区钻探资料图的一角，在每个格点上有三个数字列，都是相对于选定基点的高度，最上面的数字是覆盖表面的标高，中间的数字是沙层顶部的标高，最下面的数字是沙层底部的标高，每个格子都是正方形，边长50m。
画*处，没有数据。
A
B
C
D
E
F
G
H
0
22.4
22.5
23.0
23.2
20.0
18.4
17.8
18.0
5.8
0.5
0.4
0.4
1
23.0
23.1
23.1
23.4
23.5
24.0
24.0
24.0
19.9
20.0
20.0
19.8
19.9
20.0
19.8
19.6
6.0
3.2
1.6
1.0
1.1
1.0
0.8
0.9
2
23.1
23.2
23.4
23.4
23.5
24.2
24.1
24.1
19.8
19.7
19.4
20.0
20.1
20.3
20.3
20.5
2.2
1.4
0.6
0.5
0.3
-0.2
-0.1
0.0
2.实验内容
解假设地下沙层连续变化，于是可以利用积分运算计算矩形区域内的含沙量，首先通过插值补充缺失的数据，采用一维样条插值
>> D=[23.0,23.1,23.2,23.4,23.5,24.0,24.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,'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.5,0.4];
C1=interp1（P,C,K,'spline'）;
c=[C1;F];
>> [M,N]=meshgrid（1:
8,0:
2）;
>> delta=0.01;
>> [X,Y]=meshgrid（1:
delta:
8,0:
delta:
2）;
>> 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）*delta^2*2500;
如果在一维插值补充缺失数据中用线性插值代替上面的样条插值，将得到沙层体积634525m3,和以上结果有些差别。
如果利用二维插值，就会得到相近的结果
>>x0=[16781:
81:
8];
>>y0=[zeros（1,4）ones（1,8）2*ones（1,8）];
>>f0=[BE]'-[CF]';
>>f=griddata（x0,y0,f0,X,Y,'cubic'）;
>>v=trapz（trapz（f,1）,2）*delta^2*2500；
例2.2.3 范克里金插值
考虑沙子是一种特殊的地质，下面试用地质学常用的范克里金插值方法计算。
先介绍克里金插值方法，这本身就是运用统计学方法建模的一个有趣的例子。
1951年南非地质学家Krige将x=（x1,x2） ∈ R2处矿藏的贮藏量f（x）看成是随机函数F（x）的一个实现，提出了依据观测值（xj , fj）,j=1,2,…,n,寻找基函数｛λj（x）,j=1,2,…,n｝,以获得F（x）的具体形式
的最小方差无偏估计，取F*（x）的条件期望f*（x）=E（F*（x）|F（xj）=fj,j=1,2,…,n）给出未知点x估值的插值方法，既考虑形如
的随机函数，其中M（x）在地质学中称为漂移（drift），在一些随机分析的场合也称为趋势（trend），Pα是均值为E（Pα）=Pα的未知的随机变量，｛Φα｝是多项式基函数，例如这里取Φ1（x）=1,Φ2（x）=x1,Φ3（x）=x2,Φ4（x）=x12,Φ5（x）=x22,Φ6（x）=x1x2,R（x）是满足E（R（x））=0,E（R（x）R（x+h））=σ（h）的随机函数，σ（h）是给定的且满足当h->∞时，有σ（h）->0;当h->0时，有σ（h）->1的核函数。
这里按通常取法，取高斯核函数σ（h）=exp（-h2）。
根据无偏估计要求，
即
从而得到无偏条件，对任意的α有
在这个约束条件下极小化方差
其中A=（σ（xj-xk））n*n,D（x）（σ（xj-xk））n*1）,
σ（0）=1.
引入拉格朗日乘子
由拉格朗日函数
升级会员 