《数学实验11水塔水流量估计模型与数据插值》由会员分享，可在线阅读，更多相关《数学实验11水塔水流量估计模型与数据插值（70页珍藏版）》请在金锄头文库上搜索。
1、,Mathematics Laboratory,阮小娥博士,数学实验,办公地址：理科楼214,实验14 水塔水流量估计模型与数据插值,2、学会用matlab软件进行数据插值计算,1、掌握四种经典的插值方法：拉格朗日插值法、牛顿插值法、分段插值法、三次样条插值法,3、学会用数据插值、数据拟合方法建立数学模型并求解,实验问题,一条100米宽的河道见图1所示,为了测量其流量需要知道河道的截面积。为此从一端开始每隔5米测量出河床的深度如表1,试根据以上数据，估计出河道的截面积，进而在 已知河流速度（设为1m/s）的情况下计算出流量。 若沿河床铺设一条光缆，试估计光缆的长度。,一 数据插值,给定n个数据
2、点,称P(x)为插值函数，也即求解一条严格通过各数据点的曲线，用它来进行分析研究和预测，这种方法常称为数据插值法。,对于这类问题，选取一条何种类型的曲线作为插值函数是求解的关键。P(x)不同，就得到不同的插值函数。,1 拉格朗日插值,令,可以证明，对于n+1个不同结点，必存在唯一的次数不超过n的满足条件的多项式，这个多项式称为插值多项式，这种方法称为n次多项式插值（或代数插值）。,为了以后使用方便，先编制一个Lagrange插值函数程序： function p=lagrange(x,y) L=length(x); A=ones(L); %产生元素全是1的L阶方阵 for j=2:L A(:,j
3、)=A(:,j-1).*x; end X=inv(A)*y; for i=1:L p(i)=X(L-i+1); end,例4.1 已知观测数据,求其插值多项式曲线。,利用Matlab提供的计算以向量p为系数的多项式值 的命令polyval , 可以求得多项式函数在任意一点的值。 其使用格式为 y0=polyval(p,x0),其中p为多项式系数（从高次到底次排列）向量，x0为所 求值的点或向量，y0为返回x0处的函数值，从而可以绘制 多项式曲线。,程序如下： x=1 2 3 4 5; y=-1 1.5 2.1 3.6 4.9; plot(x,y,b.,markersize,30) axis(1
4、 5 -1 5) grid hold on p=lagrange(x,y); t=1:0.1:5; u=polyval(p,t); plot(t,u,r-,linewidth,3),从结果可以看到，所插值的4次多项式曲线较好地连接了5个 数据点，从而可以用此多项式曲线作为这5个数据的一个近 似变化。,2 牛顿插值法,拉格朗日插值法的最大缺点在于当增加差值点时，需要重新计算多项式的系数，没有承接性。下面的牛顿插值法就避免了这个问题。,差商具有如下性质：,（1）m阶差商是零阶差商的线性组合； （2）差商与插值结点的次序无关； （3）若f(x)是m次多项式，则 fx0,x1,xm=0,由差商公式，,
5、逐次代入得,称为牛顿插值公式，最后一项称为牛顿插值余项，记为Rn(x)，余项前的多项式称为插值多项式，记为Pn(x)。,牛顿插值多项式具有以下特点：,（1）在插值结点处与拉格朗日插值一样，误差为零； （2）多项式k次项的系数是f(x)的k阶差商； （3）增加插值节点时，只增加最后一项，不必像拉格朗日插值公式那样需要重新计算系数。,在做牛顿插值时，一般先做出差商表，然后套用公式。,3 样条插值法,特别地，如果m=1，则称为分段线性插值，即在每一个子区间上S(x)是线性函数，在整个区间上是分段线性函数。这种情况一般很难满足实际要求，通常使用最多的是3次样条插值函数，即S(x)在每一个子区间上是三次
6、多项式函数，而且在每个结点处满足二阶导数连续和相关的边界条件。,二 MATLAB软件实现数据插值,1 一维插值命令 yb=interp1(x,y,xb,method) x,y是同维数据向量，表示插值结点的横坐标和纵坐标。若x是向量，y是矩阵，则对y的每一列与x配对进行插值； xb表示待求函数值的插值结点向量，可以缺省； method是可选项，说明插值使用的方法。缺省时为线性插值，也可选择：nearest（最近插值），linear（线性），spline（三次样条），cubic（三次插值）。 命令返回值yb是插值曲线在xb处的纵坐标值。,2 二维插值命令 zb=interp2(x,y,z,xb,y
7、b,method) 根据同维数据向量x，y，z，按照指定的方法做插值运算，然后返回插值函数的竖坐标。,3 三维插值命令 vb=interp3(x,y,z,v,xb,yb,zb,method) 4 样条插值命令 yb=spline(x,y,xb) 该命令等同于 yb=interp1(x,y,xb,cubic),zb = griddata(x,y,z,xb,yb,method) x,y,z是同维数据向量，表示插值结点的横坐标、纵坐标和竖坐标。xb，yb表示待求函数值的插值结点向量，通常由meshgrid生成。 method是可选项，说明插值使用的方法。缺省时为线性插值，也可选择：nearest（最
8、近插值），linear（线性），cubic（三次插值），v4。 v4是一种插值算法，没有具体的名字，原文称为“MATLAB 4 griddata method”，是一种很圆滑的差值算法，效果不错。 命令返回值zb是插值曲面在(xb,yb)处的竖坐标值。,5 用griddata命令插值生成三维曲面,interp2和griddata都是二维插值，两者的区别是： interp2的插值数据必须是矩形域， 即已知数据点(x,y)组成规则的矩阵，或称之为栅格，可使用meshgid生成。 而griddata函数的已知数据点（X，Y）不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果。
9、,例4.2 已知观测数据,分别用拉格朗日、分段线性、3次样条进行插值，并绘出插值多项式曲线图。,(1)拉格朗日插值法,x=0:0.1:1； y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2; plot(x,y,b.,markersize,30) axis(0 1 -2 16) grid hold on p=lagrange(x,y); t=0:0.01:1; u=polyval(p,t); plot(t,u,r-,linewidth,3),结果显示，所插值出的10次多项式曲线，在数据点之间产生较大的起伏波动，与数据点的变化趋势
10、有明显的偏离，这时曲线并不能很好地反映数据点的变化规律。而且进一步实验发现，随着分点的增加，Lagrange插值出现大的起伏波动越明显，这就是插值问题中典型的“龙格现象”。 针对高次多项式插值时容易发生“龙格现象”，在实际插值时，常常采用分段低次插值方法，即在相邻两个数据点构成的子区间上分别进行低次插值, 整个区间上的插值函数将是一个分段的多项式函数。,(2)分段线性插值,x=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2; plot(x,y,
11、b.,markersize,30)； axis(0 1 -1 15)； grid hold on t=0:0.01:1; u=interp1(x,y,t); plot(t,u,r-,linewidth,3),结果分析,分段线性插值有效地回避了插值问题中的“龙格现象”，结果连线也大致描述了已知数据点的变化规律。但很明显，由分段直线连接的插值曲线在节点处不光滑，不可导。,(3)3次样条插值,x=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2; pl
12、ot(x,y,b.,markersize,30) axis(0 1 -1 16) grid hold on t=0:0.01:1; u=spline(x,y,t); plot(t,u,r-,linewidth,3),结果分析,从图可以看出，样条插值所得的曲线比较好地连接了已知的数据点，既有效地回避了插值问题中的“龙格现象”，又是连续光滑的，用此曲线来近似描述已知数据点的变化规律，能比较好地进行数据点之间的预测分析和求值。,例4.3 某海域上频繁地有各种吨位的船只经过，为保证船只的航行安全，有关机构在低潮时对水深进行了测量，得到了如下的测量数据： x=129.0 140.0 103.5 88.0
13、 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5; y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 60.5 121.0 3.0 56.5 116.5 84.0 43.5; Z=4 8 6 8 6 8 8 9 9 8 8 9 4 9; x,y为平面坐标，z为高度。请根据测量的数据描述该海域的地貌。,程序如下： clc; x=129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5; y=7.5
14、141.5 23.0 147.0 22.5 137.5 85.5 60.5 121.0 3.0 56.5 116.5 84.0 43.5; z=4 8 6 8 6 8 8 9 9 8 8 9 4 9; h=-z; xi=linspace(min(x),max(x),100); yi=linspace(min(y),max(y),100); X,Y=meshgrid(xi,yi); H=griddata(x,y,h,X,Y,v4); mesh(X,Y,H);,三 实验问题求解,1 画出河床观测点的散点图 x=0:5:100; y=0 2.41 2.96 2.15 2.65 3.12 4.23 5
15、.12 6.21 5.68 4.22 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0; axis(0 100 0 10) grid y1=10-y; plot(x,y1,b.,markersize,30),2 利用分段线性插值绘制河床曲线 根据已给的数据可以进行分段线性插值, 在此基础上，利用梯形法求积分命令trapz来计算河床截面积，同时计算每一段连接线长度之和来近似河床曲线长度。,x=0:5:100; y=0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 3.91 3.26 2.85 2.3
16、5 3.02 3.63 4.12 3.46 2.08 0; y1=10-y; plot(x,y1,b.,markersize,30); axis(0 100 0 10) grid hold on t=0:100; u=interp1(x,y1,t); plot(t,u,r,linewidth, 3); S=100*10-trapz(x,y1); p=sqrt(diff(x).2+diff(y1).2); L=sum(p); fprintf(S=%.2f,L=%.2fn,S,L),从运行结果得到： 河床截面面积：S=337.15， 河床曲线长度： L=102.09,3 利用样条插值计算截面积与曲线长度 另一方面, 为了提高河床曲线的模拟精度，可以根据已给的数据进行3次样条插值, 在此基础上，利用梯形法求积分命令
