数学实验：实验14 水塔水流量估计模型与数据插值

MathematicsLaboratory阮小娥博士ExperimentsinMathematics数学实验实验14水塔水流量估计模型与数据插值2、学会用matlab软件进行数据插值计算1、掌握四种经典的插值方法：拉格朗日插值法、牛顿插值法、分段插值法、三次样条插值法3、学会用数据插值、数据拟合方法建立数学模型并求解实验问题一条100米宽的河道见图1所示,为了测量其流量需要知道河道的截面积。为此从一端开始每隔5米测量出河床的深度如表1

试根据以上数据，估计出河道的截面积，进而在已知河流速度（设为1m/s）的情况下计算出流量。若沿河床铺设一条光缆，试估计光缆的长度。一数据插值给定n个数据点称P(x)为插值函数，也即求解一条严格通过各数据点的曲线，用它来进行分析研究和预测，这种方法常称为数据插值法。对于这类问题，选取一条何种类型的曲线作为插值函数是求解的关键。P(x)不同，就得到不同的插值函数。1拉格朗日插值

令

可以证明，对于n+1个不同结点，必存在唯一的次数不超过n的满足条件的多项式，这个多项式称为插值多项式，这种方法称为n次多项式插值（或代数插值）。为了以后使用方便，先编制一个Lagrange插值函数程序：functionp=lagrange(x,y)L=length(x);A=ones(L);%产生元素全是1的L阶方阵forj=2:L

A(:,j)=A(:,j-1).*x';endX=inv(A)*y';fori=1:L

p(i)=X(L-i+1);end例4.1

已知观测数据x

12345y-11.52.13.64.9求其插值多项式曲线。利用Matlab提供的计算以向量p为系数的多项式值的命令polyval

,可以求得多项式函数在任意一点的值。其使用格式为

y0=polyval(p,x0)其中p为多项式系数（从高次到底次排列）向量，x0为所求值的点或向量，y0为返回x0处的函数值，从而可以绘制多项式曲线。程序如下：x=[12345];y=[-11.52.13.64.9];plot(x,y,'b.','markersize',30)axis([15-15])gridholdonp=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次多项式，则

f[x0,x1,…,xm]=0由差商公式，逐次代入得称为牛顿插值公式，最后一项称为牛顿插值余项，记为Rn(x)，余项前的多项式称为插值多项式，记为Pn(x)。牛顿插值多项式具有以下特点：（1）在插值结点处与拉格朗日插值一样，误差为零；（2）多项式k次项的系数是f(x)的k阶差商；（3）增加插值节点时，只增加最后一项，不必像拉格朗日插值公式那样需要重新计算系数。结点零阶差商一阶差商二阶差商三阶差商在做牛顿插值时，一般先做出差商表，然后套用公式。3样条插值法

特别地，如果m=1，则称为分段线性插值，即在每一个子区间上S(x)是线性函数，在整个区间上是分段线性函数。这种情况一般很难满足实际要求，通常使用最多的是3次样条插值函数，即S(x)在每一个子区间上是三次多项式函数，而且在每个结点处满足二阶导数连续和相关的边界条件。二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,yb,'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（最近插值），linear（线性），cubic（三次插值），v4。

v4是一种插值算法，没有具体的名字，原文称为“MATLAB4griddatamethod”，是一种很圆滑的差值算法，效果不错。

命令返回值zb是插值曲面在(xb,yb)处的竖坐标值。5用griddata命令插值生成三维曲面

interp2和griddata都是二维插值，两者的区别是：interp2的插值数据必须是矩形域，即已知数据点(x,y)组成规则的矩阵，或称之为栅格，可使用meshgid生成。而griddata函数的已知数据点（X，Y）不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好的效果。例4.2已知观测数据

分别用拉格朗日、分段线性、3次样条进行插值，并绘出插值多项式曲线图。x00.10.20.30.40.50.60.70.80.91y-0.4471.9783.286.167.087.347.669.569.489.311.2(1)拉格朗日插值法x=0:0.1:1；y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,'b.','markersize',30)axis([01-216])gridholdonp=lagrange(x,y);t=0:0.01:1;u=polyval(p,t);plot(t,u,'r-','linewidth',3)

结果显示，所插值出的10次多项式曲线，在数据点之间产生较大的起伏波动，与数据点的变化趋势有明显的偏离，这时曲线并不能很好地反映数据点的变化规律。而且进一步实验发现，随着分点的增加，Lagrange插值出现大的起伏波动越明显，这就是插值问题中典型的“龙格现象”。针对高次多项式插值时容易发生“龙格现象”，在实际插值时，常常采用分段低次插值方法，即在相邻两个数据点构成的子区间上分别进行低次插值,整个区间上的插值函数将是一个分段的多项式函数。(2)分段线性插值x=[00.10.20.30.40.50.60.70.80.91];y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,‘b.’,‘markersize’,30)；axis([01-115])；gridholdont=0:0.01:1;u=interp1(x,y,t);plot(t,u,'r-','linewidth',3)结果分析分段线性插值有效地回避了插值问题中的“龙格现象”，结果连线也大致描述了已知数据点的变化规律。但很明显，由分段直线连接的插值曲线在节点处不光滑，不可导。(3)3次样条插值x=[00.10.20.30.40.50.60.70.80.91];y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,'b.','markersize',30)axis([01-116])gridholdont=0:0.01:1;u=spline(x,y,t);plot(t,u,'r-','linewidth',3)结果分析从图可以看出，样条插值所得的曲线比较好地连接了已知的数据点，既有效地回避了插值问题中的“龙格现象”，又是连续光滑的，用此曲线来近似描述已知数据点的变化规律，能比较好地进行数据点之间的预测分析和求值。例4.3某海域上频繁地有各种吨位的船只经过，为保证船只的航行安全，有关机构在低潮时对水深进行了测量，得到了如下的测量数据：x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];

y=[7.5141.523.0147.022.5137.585.560.5121.03.056.5116.584.043.5];

Z=[48686889988949];

x,y为平面坐标，z为高度。请根据测量的数据描述该海域的地貌。程序如下：clc;x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.560.5121.03.056.5116.584.043.5];z=[48686889988949];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=[02.412.962.152.653.124.235.126.215.684.22…3.913.262.852.353.023.634.123.462.080];axis([0100010])gridy1=10-y;plot(x,y1,'b.','markersize',30)2利用分段线性插值绘制河床曲线根据已给的数据可以进行分段线性插值,在此基础上，利用梯形法求积分命令trapz来计算河床截面积，同时计算每一段连接线长度之和来近似河床曲线长度。x=0:5:100;y=[02.412.962.152.653.124.235.126.215.684.22…3.913.262.852.353.023.634.123.462.080];y1=10-y;plot(x,y1,'b.','markersize',30);axis([0100010])gridholdont=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=%.2f\n',S,L)从运行结果得到：河床截面面积：S=337.15，河床曲线长度：

L=102.093利用样条插值计算截面积与曲线长度另一方面,为了提高河床曲线的模拟精度，可以根据已给的数据进行3次样条插值,在此基础上，利用梯形法求积分命令trapz来计算河床截面积，同时对样条曲线加密，计算每一段连接线长度之和来近似河床曲线长度。x=0:5:100;y=[02.412.962.152.653.124.235.126.215.684.22...3.913.262.852.353.023.634.123.462.080];y1=10-y;plot(x,y1,'b.','markersize',30)axis([0100210])gridholdont=0:0.5:100;u=spline(x,y1,t);plot(t,u,'r','linewidth',3)；S=100*10-trapz(t,u);p=sqrt(diff(t).^2+diff(u).^2);L=sum(p);fprintf('S=%.2f,L=%.2f\n',S,L)河床截面面积：S=339.50，河床曲线长度：L=102.23.数据插值模型试验：水塔水流量估计一．实验问题某一地区的用水管理机构要求各社区提供每小时以加仑计的用水量以及每天所用水的总量。由于许多社区没有测量流入或流出水塔的水量装置，所以他们只能通过测量水塔每小时的水位高度来代替每小时的用水量(误差不超过0.05)。当水塔中的水位下降到最低水位L(27.00ft)时水泵就自动启动向水塔输水,当水塔水位达到限定最高水位H(35.00ft)时，水泵停止供水。水泵供水期间无法测量水泵的供水量。水泵每天输水一次或两次,每次约二小时。

已知某一小镇的水塔是一高为40ft，直径为57ft的正圆柱，下表记录的是某一天水塔水位的真实数据(1ft=0.3024米(m))。试估计任何时刻(包括水泵正在输水时间)从水塔流出的水流量和一天的用水总量。二、问题分析

流量是单位时间内流出水的体积，由于水塔是正圆柱形，截面面积是常数，所以，流量很容易根据水位相对时间的变化率算出。

水泵不工作时段的用水量可以由测量记录直接得到:由表1中记录的水位下降高度乘以水塔的截面面积就是这一时段的用水量。

问题的难点在于:如何估计水泵供水时段的流量。

水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。而作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量数据越准确越好。这个数值可以用来检验数据插值或拟合结果的准确性。二、问题分析流量大体上可由两种方法计算:

直接对表1中的用水量用数值插值方法算出各时段的流量，用它们拟合其它时刻或连续时间的流量;

先用表中数据拟合水位--时间函数,求导数即可得到连续时间的流量。有了任何时刻的流量，就不难计算一天的总用水量。

三、问题求解

为了表示方便，将所给表1中的数据全部化为国际标准单位，时间用小时(h),高度用米(m)。3.1模型假设(1)流量只取决于水位差，与水位本身无关。(2)Torricelli定律：从小孔流出液体的流速正比于液面高度的平方根。在假设出口的水位为零的前提下，题目给出水塔的最低和最高水位分别是8.1648m(27×0.3024)和10.7352m(35.50×0.3024)。因为：sqrt(10.7352/8.1648)=1.1467≈1所以，可忽略水位对流速的影响。(3)流量可以看作是时间的连续函数。因此，为了计算简单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值。(4)水塔截面面积为3.2流量估计方法由表2给出的数据，用MTALAB软件做出时间—水位散点图1：为了计算水箱水流量与时间的关系，最简单的方法就是：

根据给出的散点数据图，将数据分为三段，然后对每一段数据用数据插值或者数据拟合的方法得到流量与时间的近似函数关系。具体方法如下：中点流速=(左端点的水位-右端点的水位)/时间区间长度每段数据首尾点的流速用下面的公式计算：经过计算，得到时间与流速之间的关系数据表3数据表3对应的时间--流速散点图如下：下面分别用数据插值和数据拟合两种方法来估计水塔水流量。(1)数据插值法采用拉格朗日、分段线性和三次样条三种数据插值方法。根据表3的数据，对水泵不工作时段1,2通过采取数据插值方法可以得到任意时刻的流速，从而知道任意时刻的流量。对于水泵工作时段1，应用前后时期的流速进行插值。由于最后水泵工作时段2的数据太少，我们将它与水泵不工作时段3合并，一同进行数据插值处理（简称混合时段）。以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。t=[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97];v=[29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.27];t0=0:0.1:8.97;lglr=lglrcz(t,v,t0);%拉格朗日插值，需要定义m文件lglrjf=0.1*trapz(lglr)%梯形法算数值积分fdxx=interp1(t,v,t0);%分段线性插值fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,‘spline’);%三次样条插值sancytjf=0.1*trapz(scyt)plot(t,v,'*',t0,lglr,'r',t0,fdxx,'g',t0,scyt,'b')gtext('lglr')gtext('fdxx')gtext('scyt')计算结果：lglrjf=145.6231，fdxxjf=147.1430，sancytjf=145.6870functiony=lglrcz(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;end

y(i)=s;end由表2知，第1未供水时段的总用水高度为146(=968-822)。上述三种插值方法计算的结果与实际值146相比都比较接近。计算结果：lglrjf=145.6231，fdxxjf=147.1430，sancytjf=145.6870(2)数据拟合法1)拟合水位--时间函数2)确定流量—时间函数3)一天总用水量的估计4)流量及总用水量的检验从表2中测量记录看，一天有两次供水时段和三次未供水时段，分别对第1,2未供水时段的测量数据直接作多项式拟合，可得到水位函数。第3未供水时段数据过少,拟合效果不是很理想。应当注意:根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以3-6次为宜）。1)拟合水位--时间函数t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95…12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96…20.84,23.88,24.99,25.66];h=[9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82…10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59…10.35,10.18];figure(1)c1=polyfit(t(1:10),h(1:10),3);tp1=0:0.1:8.9;x1=polyval(c1,tp1);%变量x1中存放了以0.1为步长算出的各个时刻的水位高度。plot(tp1,x1)figure(2)c2=polyfit(t(11:21),h(11:21),3);tp2=10.9:0.1:20.9;x2=-polyval(c2,tp2);plot(tp2,x2)第1、2未供水时段时间与水位图程序：2)确定流量—时间函数第1,2未供水时段第1供水时段第2供水时段＋第3未供水时段＝混合时段第1,2未供水时段的流量可直接对水位函数求导得到，用三次多项式拟合曲线，程序如下：c1=polyfit(t(1:10),h(1:10),3);c2=polyfit(t(11:21),h(11:21),3);a1=polyder(c1);a2=polyder(c2);tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x113=-polyval(a1,[0:0.01:8.97]);wgsysl1=100*trapz(tp1,x113);x14=-polyval(a1,[7.93,8.97]);%(为后面的计算准备数据)x23=-polyval(a2,tp2);x114=-polyval(a2,[10.95:0.01:20.84]);wgsysl2=100*trapz(tp2,x114);x24=-polyval(a2,[10.95,12.03]);%(为后面的计算准备数据)x25=-polyval(a2,[19.96,20.84]);%(为后面的计算准备数据)subplot(1,2,1)plot(tp1,x13*100)subplot(1,2,2)plot(tp2,x23*100)上下曲线比较发现，用三次多项式拟合效果不是很好，故改用五次多项式重新拟合。对于第1供水时段的流量,则用前后时期的流量进行拟合得到。为使流量函数在t=8.97和t=10.93连续，只取4个点，用3次多项式拟合得到第1供水时段的时间-流量图。dygsdsj=[7.93,8.97,10.95,12.03];dygsdls=[x14,x24];nhjg=polyfit(dygsdsj,dygsdls,3);nhsj=7.93:0.1:12.03;nhlsjg=polyval(nhjg,nhsj);gssj1=8.97:0.01:10.95;gs1=polyval(nhjg,[8.97:0.01:10.95]);gsysl1=100*t