数学建模——水塔流量问题

212实验十四水塔流量问题【实验目的】1．了解有关数据处理的基本概念和原理。2．初步了解处理数据插值与拟合的基本方法，如样条插值、分段插值等。3．学习掌握用MATLAB命令处理数据插值与拟合问题。【实验内容】某居民区有一供居民用水的圆形水塔，一般可以通过测量其水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间是无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约两小时。水塔是一个高12.2米、直径17.4米的正圆柱。按照设计，水塔水位降到约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作。某一天的水位测量记录如表1所示，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量。表1水位测量启示录（//表示水泵启动）时刻（h）水位（cm）09680.929481.849312.959133.878984.988815.908697.018527.938398.97822时刻（h）水位（cm）9.98//10.92//10.95108212.03105012.95102113.8899414.9896515.9094116.8391817.93892时刻（h）水位（cm）19.0486619.9684320.8482222.01//22.96//23.88105924.99103525.911018【实验准备】在生产实践和科学研究中，常常遇到这样的问题：由实验或测量得到的一批离散样点，需要确定满足特定要求的曲线或曲面（即变量之间的函数关系或预测样点之外的数据）。如果要求曲线（面）通过所给的所有数据点（即确定一个初等函数通过已知各数据，一般用多项式或分段多项式），这就是数据插值。在数据较少的情况下，这样做能够取得好的效果。但是，如果数据较多，那么插值函数是一个次数很高的函数，比较复杂。如果不要求曲线（面）通过所有的数据点，而是要求它反映对象整体的变化趋势，可得到更简单实用的近似函数，这就是数据拟合。函数插值和曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。1．数据插值的基本方法拉格朗日插值若知道函数y＝)(xf在互异的两个点0x和1x处的函数值0y和1y，而想估计该函数在另一点处的函数值，最自然的想法是作过点（0x，0y）和点（1x，1y）的直线y＝)(1xL，用)(1L作为准确值的近似值，如果得到的结果误差太大，还可增加一点)(xf的函数值，即已知y＝)(xf在互异的三个点0x，1x和2x处的函数值0y，1y和2y，可以构造过这三点的二次曲线y＝)(2xL，用)(2L作为准确值)(f的近似值。213一般的，若已知y＝)(xf在互异的n＋1个点0x，1x，…，nx处的函数值0y，1y，…，ny，则可以考虑构造一个过这n＋1个点的次数不超过n的多项式)(xLn)(xLn＝mxa0＋11mxa＋…＋xam1＋ma（1）通过所有n＋1个点，即满足)(knxL＝ky，k＝0，1，…，n（2）然后用)(nL作为准确值)(f的近似值。这样构造出来的多项式)(xLn称为)(xf的n次拉格朗日插值多项式或插值函数。分段插值多项式历来都被认为是最好的逼近工具之一，它插值光滑，但不具有收敛性，会随着节点数目增多而次数升高，一般不宜采用高次多项式（如m＞7）插值，否则逼近的效果往往是不理想的，甚至发生龙格振荡（当节点数目n不断增大时，)(xLn在区间中部趋于)(xf，但对于区间两端的x，)(xLn并不趋于)(xf，也称龙格现象）。在插值范围较小，用低次插值往往就能奏效。最直观的办法就是将各数据点用折线连接起来，这种增加节点，用分段低次多项式插值的化整为零的处理方法称作分段插值法，即不去寻求整个插值区间上的一个高次多项式，而是把区间划分为若干个小区间。如果a＝0x＜1x＜…＜nx＝b（3）那么分段线性插值公式为)(xP＝11iiiiyxxxx＋iiiiyxxxx11，1ix＜x≤ix，i＝0，1，…，n（4）分段线性插值通常有较好的收敛性和稳定性，算法简单，克服了龙格现象，其缺点是不如拉格朗日插值多项式光滑。样条插值分段线性插值函数在节点的一阶导数一般不存在，且不光滑，这就导致了样条插值函数的提出。在机械制造、航海、航空工业中，经常需要解决下列问题：已知一些数据点（0x，0y），（1x，1y），（nx，ny），如何全部通过这些数据点作一条比较光滑的曲线呢？绘图员解决了这一问题，首先把数据描绘在平面上，再把一根富有弹性的细直条（称为样条）弯曲，使其一边通过这些数据点，用压铁固定其形状，沿样条边绘出一条光滑的曲线，往往要用几根样条，分段完成上述工作，同时也应让连接点处保持光滑。对绘图员用样条画出的曲线，进行数学模拟，就导出了样条函数的概念。如今已经成为了一个应用极为广泛的数学分支。现在数学上所说的样条，实质上指分段多项式的光滑连接。设有区间［a，b］的一个划分如（3）式，称分段函数)(xS为k次样条函数，若它有：（1）)(xS在每个小区间上的次数不超过k多项式；（2）)(xSi＝iy（3）)(xS在区间［a，b］上有k－1阶连续的导数；用样条函数作出的插值称为样条插值，工程上广泛采用三次样条插值。2．曲线拟合的基本方法曲线拟合问题是指：已知平面上n个点（ix，iy），i＝0，1，…，n，ix互不相同，寻求函数y＝)(xf，使)(xf在某种准则下与所有数据点最为接近，即曲线拟合得最好。线性最小二乘法是解决曲线拟合最常用的方法，其基本思路是，令)(xf＝)(11xra＋)(22xra＋…＋)(xramm（5）其中)(xrk是事先选定的一组函数，系数ka（k＝0，1，…，m,mn）待定。寻求ka，使得残差平方和Q＝niixf12i)y)((（6）214达到最小。这里的建模原理实质上与实验七中的回归分析是一致的。3．数据插值与拟合的MATLAB命令对于多项式插值和拟合，有一个方便的方法p=polyfit(x,y,m)用m次多项式拟合向量数据（x，y），返回多项式（1）的降冥系数。当m≥n－1时，polyfit实现多项式插值，p返回多项式的系数向量；y=polyval(p,x)求polyfit所得的多项式在x处的插值y，它是准确值f（x）的近似值；对于一维和二维插值，其命令格式如下yi=interp1(x,y,xi,method)求解一维插值问题，x，y分别表示数据点的横、纵坐标向量，且x要单调。xi为插值点，它不能走出x的范围。method为可选参数，对应四种插值方法：最近邻点插值：'nearest'；线性插值：'linear'；（是缺省值）三次样条插值：'spline'；分段三次插值：'cubic'ZI=interp2(X,Y,Z,XI,YI,method)求解二维插值，X，Y分别是m、n维族自变量，均单调递增；Z是m×n维矩阵，标明相应于所给数据的函数值。向量XI,YI给定网格点的横、纵坐标，ZI返回函数在网格坐标（XI，YI）处的函数值。XI，YI应是方向不同的向量，即一个是行向量，一个是列向量。method的定义与interp1命令中的定义是一致的；ZI=griddata(x,y,z,XI,YI,method)插值基点为散乱的节点，各参数定义与二维插值中一致，只不过向量x，y散乱无序，同时method方法中多了一种MATLAB提供的网格插值方法V4；有关上述命令的详细内容可在MATLAB帮助文件中查阅。对于线性最小二乘拟合，用得较多的是多项式拟合，使用前面所讲的polyfit命令；若要寻求f（x）任意的非线性函数，则称为非线性最小二乘拟合，MATLAB提供了两个求解命令：curfit和leastsq。二者都要事先定义M－函数文件，但定义方式稍有不同：c=curvefit('fun',c0,xdata,ydata，options)求含参量非线性函数fun中的参变量c，使残差平方和（6）最小，xdata，ydata为数据点的横、纵坐标向量，c0为参变量的迭代初始值，options见实验一表1（它可以缺省）；非线性函数fun的M－函数文件定义方式为：fun(c,xdata)；c=leastsq('fun',c0,options)求含参量非线性函数fun中的参变量c，使得各数据点函数值fun的平方和最小，命令中各参数定义与curvefit命令中一致，非线性函数fun的M－函数文件定义方式为：fun(c,xdata,ydata)，这里的fun已经与数据点的函数值向量ydata作差；【实验方法与步骤】1．引例问题的分析流量是单位时间流出的水的体积，由于水塔是圆柱形，横截面积是常数，在水泵不工作时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量拟合得到，作为拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。我们可以考虑先用表中数据拟合水位～时间函数，215然后对之求导即可得到各时段的流量。有了任意时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如由某一时段水位下降量乘以水塔的截面积（水塔截面积是常数S＝(17.4/2)2＝237.8（m2））就得到这一时段的用水量。这个数值可以还可以用来检验拟合效果。流量是时间的连续函数，只取决于水位差，与水位本身无关，与水泵是否工作无关。按照Torricelli定律从小孔流出的液体的速度正比于水面高度的平方根，题目给出水塔的最高和最低水位分别为10.8米和8.2米（设出水口的水位为0），因为2.8/8.10＝1.15，可以忽略水位对流速的影响。简单起见，计算中将流量定义为单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的）。水泵第1次供水时段为t＝9.0到t＝11.0（小时），第2次供水时段为t＝20.8到t＝23.0（小时）。这是根据最高和最低水位分别为10.8米和8.2米，及表1的水位测量记录作出的假设，其中前3个时刻直接取自实测数据（精确到0.1小时），最后1个时刻来自每次供水约两小时的已知条件（从记录看，第2次供水时段应在记录的22.96小时之后不久结束）。水泵工作时单位时间的供水量大致为常数，这个常数应该大于单位时间的平均流量。首先考虑拟合水位～时间函数，从表1测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和三个水泵不工作时段（简称第1时段t＝0到t＝8.97，第2时段t＝10.95到t＝20.84，第3时段t＝23以后）。对第1、2时段的测量数据可直接分别作多项式拟合，得到水位函数。为使拟合曲线比较光滑，多项式次数不要太高，一般在3～6次。由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合。接着确定流量～时间函数，对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内。最后一天总用水量等于两个水泵不工作时段和两个供水时段（将第3时段包含在第2供水时段内）用水量之和，它们都可以由流量对时间的积分再乘以水塔截面积得到。2．MATLAB命令求解拟合第1、2时段的水位，并导出流量，t，h为时刻和水位测量记录（水泵启动的4个时刻不输入），程度代码如下：t=[00.921.842.953.874.985.907.017.938.9710.9512.0312.9513.8814.9815.9016.8317.9319.0419.9620.8423.8824.9925.91];h=[968948931913898881869852839822108210501021994965941918892866843822105910351018];c1=polyfit(t(1:10),h(1:10),3);%用3次多项式