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

1、实验十四水塔流量问题【实验目的】1了解有关数据处理的基本概念和原理。2初步了解处理数据插值与拟合的基本方法，如样条插值、分段插值等。3学习掌握用MA TLAB命令处理数据插值与拟合问题。【实验内容】某居民区有一供居民用水的圆形水塔， 一般可以通过测量其水位来估计水的流量。 但面临的困难是， 当水塔水位下降到设定的最低水位时， 水泵自动启动向水塔供水， 到设定的最高水位时停止供水，这段时间是无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约两小时。水塔是一个高12.2 米、直径17.4 米的正圆柱。按照设计，水塔水位降到约 8.2 米时，水泵自动启动，水位升到约10.8 米时水泵停

2、止工作。某一天的水位测量记录如表 1 所示，试估计任何时刻 （包括水泵正供水时） 从水塔流出的水流量，及一天的总用水量。表 1水位测量启示录（ /表示水泵启动）时刻（ h）00.921.842.953.874.985.907.017.938.97水位（ cm）968948931913898881869852839822时刻（ h）9.9810.9210.9512.0312.95 13.88 14.9815.9016.8317.93水位（ cm）/108210501021994965941918892时刻（ h）19.0419.9620.8422.0122.9623.8824.9925.91水位

3、（ cm）866843822/105910351018【实验准备】在生产实践和科学研究中，常常遇到这样的问题：由实验或测量得到的一批离散样点，需要确定满足特定要求的曲线或曲面（即变量之间的函数关系或预测样点之外的数据）。如果要求曲线 （面）通过所给的所有数据点 （即确定一个初等函数通过已知各数据，一般用多项式或分段多项式） ，这就是 数据插值 。在数据较少的情况下，这样做能够取得好的效果。但是，如果数据较多， 那么插值函数是一个次数很高的函数，比较复杂。如果不要求曲线 （面）通过所有的数据点，而是要求它反映对象整体的变化趋势，可得到更简单实用的近似函数，这就是 数据拟合 。函数插值和曲线拟合都

4、是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。1数据插值的基本方法拉格朗日插值若知道函数 y f (x) 在互异的两个点x0 和 x1 处的函数值 y0 和 y1 ，而想估计该函数在另一点处的函数值， 最自然的想法是作过点（ x0 ， y0 ）和点（ x1 ， y1 ）的直线 y L1 ( x) ，用 L1() 作为准确值的近似值，如果得到的结果误差太大，还可增加一点f ( x) 的函数值，即已知 y f ( x) 在互异的三个点x0 ， x1 和 x2 处的函数值y0 ， y1 和 y2，可以构造过这三点的二次曲线 y L2 ( x) ，用 L2 (

5、) 作为准确值 f ( ) 的近似值。k 多项式；精选文库一般的，若已知 y f (x) 在互异的 n 1 个点 x0 ，x1 ， xn 处的函数值yn ，则可以考虑构造一个过这n 1 个点的次数不超过n 的多项式 Ln ( x)L n ( x) a0 x m 1m 1 am 1xam通过所有 n 1 个点，即满足a xL n ( xk ) yk ， k 0， 1， n然后用 Ln ( ) 作为准确值 f () 的近似值。这样构造出来的多项式L n ( x) 称为拉格朗日插值多项式或插值函数 。分段插值y0 ， y1 ，（ 1）（ 2）f ( x) 的 n 次多项式历来都被认为是最好的逼近工具

6、之一，它插值光滑， 但不具有收敛性， 会随着节点数目增多而次数升高，一般不宜采用高次多项式（如m 7）插值，否则逼近的效果往往是不理想的， 甚至发生 龙格振荡 （当节点数目 n 不断增大时， Ln ( x) 在区间中部趋于f ( x) ，但对于区间两端的x ， Ln ( x) 并不趋于 f (x) ，也称 龙格现象 ）。在插值范围较小， 用低次插值往往就能奏效。最直观的办法就是将各数据点用折线连接起来， 这种增加节点， 用分段低次多项式插值的化整为零的处理方法称作分段插值法，即不去寻求整个插值区间上的一个高次多项式，而是把区间划分为若干个小区间。如果a x0 x1 xn b（ 3）那么分段线性

7、插值公式为P(x) xxiyi1 xxi1yi ， xi1 x xi ， i 0， 1， n（ 4）xi 1xixixi1分段线性插值通常有较好的收敛性和稳定性，算法简单， 克服了龙格现象， 其缺点是不如拉格朗日插值多项式光滑。样条插值分段线性插值函数在节点的一阶导数一般不存在，且不光滑， 这就导致了样条插值函数的提出。在机械制造、航海、航空工业中，经常需要解决下列问题：已知一些数据点（x0 ，y0 ），（ x1 ， y1），（ xn ， yn ），如何全部通过这些数据点作一条比较光滑的曲线呢？绘图员解决了这一问题，首先把数据描绘在平面上，再把一根富有弹性的细直条（称为样条 ）弯曲，使其一边通

8、过这些数据点，用压铁固定其形状， 沿样条边绘出一条光滑的曲线，往往要用几根样条， 分段完成上述工作，同时也应让连接点处保持光滑。对绘图员用样条画出的曲线， 进行数学模拟， 就导出了样条函数的概念。如今已经成为了一个应用极为广泛的数学分支。现在数学上所说的样条，实质上指分段多项式的光滑连接。设有区间 a ， b 的一个划分如 （3）式，称分段函数S(x) 为 k 次样条函数 ，若它有：（ 1） S( x) 在每个小区间上的次数不超过（ 2） Si (x) yi（ 3） S( x) 在区间 a ， b 上有 k 1 阶连续的导数；用样条函数作出的插值称为样条插值 ，工程上广泛采用三次样条插值。2曲

9、线拟合的基本方法曲线拟合问题是指：已知平面上 n 个点（ xi ， yi ）， i 0， 1， n ， xi 互不相同，寻求函数 y f ( x) ，使 f (x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。线性最小二乘法是解决曲线拟合最常用的方法，其基本思路是，令f (x) a1r1 ( x) a 2 r2 (x) am rm (x)（ 5）其中 rk (x) 是事先选定的一组函数，系数ak （ k 0， 1， m , m < n ）待定。寻求 a k ，n( f (xi ) y i ) 2使得残差平方和Q （ 6）i 1-213精选文库达到最小。这里的建模原理实质上与实验七

10、中的回归分析是一致的。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 要单

11、调。 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 ）处的

12、函数值。XI ，YI 应是方向不同的向量，即一个是行向量，一个是列向量。method 的定义与interp1 命令中的定义是一致的；ZI = griddata( x , y , z , XI , YI , method )插值基点为散乱的节点，各参数定义与二维插值中一致，只不过向量 x， y 散乱无序，同时 method 方法中多了一种 MATLAB 提供的网格插值方法 V4 ；有关上述命令的详细内容可在MATLAB帮助文件中查阅。对于线性最小二乘拟合，用得较多的是多项式拟合，使用前面所讲的polyfit命令；若要寻求 f（ x）任意的非线性函数，则称为非线性最小二乘拟合，MATLAB提供了两

13、个求解命令： 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)求含参量非线性函数fu

14、n 中的参变量c，使得各数据点函数值fun 的平方和最小，命令中各参数定义与curvefit 命令中一致，非线性函数fun 的 M 函数文件定义方式为： fun( c , xdata , ydata) ，这里的 fun 已经与数据点的函数值向量 ydata 作差；【实验方法与步骤】1引例问题的分析流量是单位时间流出的水的体积， 由于水塔是圆柱形， 横截面积是常数， 在水泵不工作时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量拟合得到， 作为拟合的原始数据， 我们希望水泵不工作时段的流量越准确越好。我们可以考虑先用表中数据拟合

15、水位时间函数，-214精选文库然后对之求导即可得到各时段的流量。有了任意时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如由某一时段水位下降量乘以水塔的截面积（水塔截面积是常数S (17.4/2)2 237.8（ m2）就得到这一时段的用水量。这个数值可以还可以用来检验拟合效果。流量是时间的连续函数，只取决于水位差，与水位本身无关，与水泵是否工作无关。按照 Torricelli 定律从小孔流出的液体的速度正比于水面高度的平方根，题目给出水塔的最高和最低水位分别为 10.8 米和 8.2米（设出水口的水位为0），因为10.8 / 8.2 1.15，可以忽

16、略水位对流速的影响。简单起见， 计算中将流量定义为单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的）。水泵第 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 小时之后不久结束） 。水泵工作时单位时间的供水量大致为常数，这个常数应该大于单位时间的平均流量。首先考虑拟合水位时间函

17、数，从表1 测量记录看，一天有两个供水时段（以下称第1供水时段和第2 供水时段），和三个水泵不工作时段（简称第1时段 t 0 到 t 8.97，第 2时段 t 10.95到 t 20.84，第 3时段 t 23以后）。对第1、2时段的测量数据可直接分别作多项式拟合，得到水位函数。为使拟合曲线比较光滑，多项式次数不要太高，一般在36 次。由于第 3 时段只有3 个测量记录，无法对这一时段的水位作出较好的拟合。接着确定流量时间函数，对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量， 则用供水时段前后（水泵不工作时段）的流量拟合得到，并将拟合得到的第2 供水时段流量外推，将第3 时段

18、流量包含在第2 供水时段内。最后一天总用水量等于两个水泵不工作时段和两个供水时段（将第3 时段包含在第 2供水时段内）用水量之和，它们都可以由流量对时间的积分再乘以水塔截面积得到。2 MATLAB 命令求解拟合第 1、 2 时段的水位，并导出流量，t， h 为时刻和水位测量记录（水泵启动的4 个时刻不输入），程度代码如下：>> t=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.8317.93 19.04 19.96 20.84 23.88 24.99 25.

19、91;>> h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;>> c1=polyfit(t(1:10),h(1:10),3);%用 3 次多项式拟合第1 时段的水位>> a1=polyder(c1);%对拟合的多项式求导数得到第1 时段流量>> tp1=0:0.1:9;%对第 1 时段的时刻进行划分>> x1=abs(polyval(a1,tp1);% 计算第 1 时段各时刻的

20、流量类似地，可计算第 2 时段各时刻的流量。>> c2=polyfit(t(11:21),h(11:21),3);>> a2=polyder(c2);>> tp2=11:0.1:20.8;>> x2=abs(polyval(a2,tp2);在第 1 供水时段（ t 9 11）之前（即第1 时段）和之后（即第2 时段）各取几点，其流量已经得到，用它们拟合水泵第1 供水时段的流量。为使流量函数在t 9 和 t 11 连续，我们简单地只取4 个点，拟合3 次多项式（即曲线必过这4 个点），实现如下：>> xx1=abs(polyval(a1

21、,8 9);>> xx2=abs(polyval(a2,11 12);-215精选文库>> xx12=xx1,xx2;>> c12=polyfit(8 9 11 12,xx12,3);%拟合水泵供水时段的流量函数>> tp12=9:0.1:11;>> x12=polyval(c12,tp12); %计算第 1 供水时段各时刻的流量在第 2 供水时段之前取 t 20， 20.8 两点的流水量，第 3 时段仅有 3 个水位记录，我们用差分得到流量，然后用这 4 个数值拟合第 2 供水时段的流量：>> dt3=diff(t(22

22、:24);%最后 3 个时刻的两两之差>> dh3=diff(h(22:24);%最后 3 个水位的两两之差>> dht3=-dh3./dt3;%用差分计算t（ 22）和 t（ 23）的流量>> t3=20 20.8 t(22) t(23);%取第 2 时段 20， 20.8 两点和第 3 时段 23.88， 24.99 两点>> xx3=abs(polyval(a2,t3(1:2),dht3;取第 2 时段 20， 20.8 两点和第 3 时段 23.88，24.99两点的流量>> c3=polyfit(t3,xx3,3)%拟合出第

23、 2 水泵供水时段的流量函数>> tp3=20.8:0.1:24;>> x3=polyval(c3,tp3);%输出第 2 供水时段（外推到t 24）各时刻的流量求第 1、2 时段和第 1、2 供水时段流量的积分之和，就是一天总用水量。虽然诸时段的流量已表示为多项式函数，积分可以解析地算出，这里仍可用数值积分计算：>> y1=0.1*trapz(x1)%第 1 时段用水量， 0.1 为积分步长y1 =146.1815>> y2=0.1*trapz(x2) %第 2 时段用水量y2 =258.0441>> y12=0.1*trapz(x

24、12) %第 1 水泵供水时段用水量y12 =50.3990>> y3=0.1*trapz(x3) % 第 2 水泵供水时段用水量 y3 =74.9138总用水量为水位差乘以水塔截面积，0.01 是因为流量单位为厘米y =1.2592e+003【结果分析】计算出来的各时段用水量可以用测量记录来检验，y1 可用第1 时段水位测量下降高度为 968 822 146 来检验，类似地，y2 用 1082 822 260 来检验。供水时段流量的一种检验方法如下：供水时段用水量加上水位上升值260 是该时段泵入的水量，除以时间长度得到水泵的功率（单位时间泵入水量），而两个供水时段的功率应大致相等。第1、 2 时段水泵的功率计算如下：>> p1=(y12+260)/2p1 =155.1995>> tp2=20.8:0.1:23;-216精选文库>> xp2=polyval(c3,tp2);>> p2=(0.1*trapz(x3)+260)/2.2 p2 =152.2335可以看到，两次水泵泵水的功率差别不大。下面是水塔一天的流量曲线图：34323