数学建模_水塔流量问题
实验十四水塔流量问题
【实验目的】
1．了解有关数据处理的基本概念和原理。
2．初步了解处理数据插值与拟合的基本方法，如样条插值、分段插值等。
3．学习掌握用MATLAB命令处理数据插值与拟合问题。
【实验内容】
某居民区有一供居民用水的圆形水塔，一般可以通过测量其水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间是无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约两小时。水塔是一个高12.2米、直径17.4米的正圆柱。按照设计，水塔水位降到约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作。
某一天的水位测量记录如表1所示，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量。
【实验准备】
在生产实践和科学研究中，常常遇到这样的问题：由实验或测量得到的一批离散样点，需要确定满足特定要求的曲线或曲面（即变量之间的函数关系或预测样点之外的数据）。如果要求曲线（面）通过所给的所有数据点（即确定一个初等函数通过已知各数据，一般用多项式或分段多项式），这就是数据插值。在数据较少的情况下，这样做能够取得好的效果。但是，如果数据较多，那么插值函数是一个次数很高的函数，比较复杂。如果不要求曲线（面）通过所有的数据点，而是要求它反映对象整体的变化趋势，可得到更简单实用的近似函数，这就是数据拟合。函数插值和曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。
1．数据插值的基本方法
拉格朗日插值
若知道函数y ＝)(x f 在互异的两个点0x 和1x 处的函数值0y 和1y ，而想估计该函数在另一点ξ处的函数值，最自然的想法是作过点（0x ，0y ）和点（1x ，1y ）的直线y ＝)(1x L ，用)(1ξL 作为准确值的近似值，如果得到的结果误差太大，还可增加一点)(x f 的函数值，即已知y ＝)(x f 在互异的三个点0x ，1x 和2x 处的函数值0y ，1y 和2y ，可以构造过这三点的二次曲线y ＝)(2x L ，用)(2ξL 作为准确值)(ξf 的近似值。
一般的，若已知y ＝)(x f 在互异的n ＋1个点0x ，1x ，…，n x 处的函数值0y ，1y ，…，
n y ，则可以考虑构造一个过这n ＋1个点的次数不超过n 的多项式)(x L n
)(x L n ＝m
x a 0＋1
1-m x
a ＋…＋x a m 1-＋m a （1）
通过所有n ＋1个点，即满足
)(k n x L ＝k y ，k ＝0，1，…，n （2） 然后用)(ξn L 作为准确值)(ξf 的近似值。这样构造出来的多项式)(x L n 称为)(x f 的n 次拉格朗日插值多项式或插值函数。 分段插值
多项式历来都被认为是最好的逼近工具之一，它插值光滑，但不具有收敛性，会随着节点数目增多而次数升高，一般不宜采用高次多项式（如m ＞7）插值，否则逼近的效果往往是不理想的，甚至发生龙格振荡（当节点数目n 不断增大时，)(x L n 在区间中部趋于)(x f ，但对于区间两端的x ，)(x L n 并不趋于)(x f ，也称龙格现象）。
在插值范围较小，用低次插值往往就能奏效。最直观的办法就是将各数据点用折线连接起来，这种增加节点，用分段低次多项式插值的化整为零的处理方法称作分段插值法，即不去寻求整个插值区间上的一个高次多项式，而是把区间划分为若干个小区间。如果 a ＝0x ＜1x ＜…＜n x ＝b （3） 那么分段线性插值公式为 )(x P ＝
11----i i i i y x x x x ＋i i i i y x x x x 1
1
----，1-i x ＜x ≤i x ，i ＝0，1，…，n （4）
分段线性插值通常有较好的收敛性和稳定性，算法简单，克服了龙格现象，其缺点是不如拉
格朗日插值多项式光滑。 样条插值
分段线性插值函数在节点的一阶导数一般不存在，且不光滑，这就导致了样条插值函数的提出。在机械制造、航海、航空工业中，经常需要解决下列问题：已知一些数据点（0x ，
0y ），（1x ，1y ），（n x ，n y ），如何全部通过这些数据点作一条比较光滑的曲线呢？
绘图员解决了这一问题，首先把数据描绘在平面上，再把一根富有弹性的细直条（称为样条）弯曲，使其一边通过这些数据点，用压铁固定其形状，沿样条边绘出一条光滑的曲线，往往要用几根样条，分段完成上述工作，同时也应让连接点处保持光滑。对绘图员用样条画出的曲线，进行数学模拟，就导出了样条函数的概念。如今已经成为了一个应用极为广泛的数学分支。现在数学上所说的样条，实质上指分段多项式的光滑连接。
设有区间［a ，b ］的一个划分如（3）式，称分段函数)(x S 为k 次样条函数，若它有：
（1）)(x S 在每个小区间上的次数不超过k 多项式； （2）)(x S i ＝i y
（3）)(x S 在区间［a ，b ］上有k －1阶连续的导数；
用样条函数作出的插值称为样条插值，工程上广泛采用三次样条插值。
2．曲线拟合的基本方法
曲线拟合问题是指：已知平面上n 个点（i x ，i y ），i ＝0，1，…，n ，i x 互不相同，寻求函数y ＝)(x f ，使)(x f 在某种准则下与所有数据点最为接近，即曲线拟合得最好。 线性最小二乘法是解决曲线拟合最常用的方法，其基本思路是，令
)(x f ＝)(11x r a ＋)(22x r a ＋…＋)(x r a m m （5） 其中)(x r k 是事先选定的一组函数，系数k a （k ＝0，1，…，m ,m ∑=-n i i x f 1 2 i ) y )(( （6） 达到最小。这里的建模原理实质上与实验七中的回归分析是一致的。 3．数据插值与拟合的MATLAB 命令 对于多项式插值和拟合，有一个方便的方法 对于一维和二维插值，其命令格式如下 对于线性最小二乘拟合，用得较多的是多项式拟合，使用前面所讲的polyfit 命令；若要寻求f （x ）任意的非线性函数，则称为非线性最小二乘拟合，MATLAB 提供了两个求解命令：curfit 和leastsq 。二者都要事先定义M －函数文件，但定义方式稍有不同： 【实验方法与步骤】 1．引例问题的分析 流量是单位时间流出的水的体积，由于水塔是圆柱形，横截面积是常数，在水泵不工作时段，流量很容易从水位对时间的变化率算出，问题是如何估计水泵供水时段的流量。 水泵供水时段的流量只能靠供水时段前后的流量拟合得到，作为拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。我们可以考虑先用表中数据拟合水位～时间函数，然后对之求导即可得到各时段的流量。有了任意时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，如由某一时段水位下降量乘以水塔的截面积（水塔截面积是常数S＝(17.4/2)2 ＝237.8（m2））就得到这一时段的用水量。这个数值可以还可以用来检验拟合效果。 流量是时间的连续函数，只取决于水位差，与水位本身无关，与水泵是否工作无关。按照Torricelli定律从小孔流出的液体的速度正比于水面高度的平方根，题目给出水塔的最高 10＝1.15，可和最低水位分别为10.8米和8.2米（设出水口的水位为0），因为2.8/8. 以忽略水位对流速的影响。简单起见，计算中将流量定义为单位时间流出的水的高度，即水位对时间变化率的绝对值（水位是下降的）。 水泵第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=[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.83 17.93 19.04 19.96 20.84 23.88 24.99 25.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时段各时刻的流量 类似地，可计算第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,[8 9])); >> xx2=abs(polyval(a2,[11 12])); >> 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: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)% 拟合出第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(x12) % 第1水泵供水时段用水量 y12 = 50.3990 >> y3=0.1*trapz(x3) % 第2水泵供水时段用水量 y3 = 74.9138 >> y=(y1+y2+y12+y3)*237.8*0.01% 总用水量为水位差乘以水塔截面积，0.01是因为流量单位为厘米 y = 1.2592e+003 【结果分析】 计算出来的各时段用水量可以用测量记录来检验，y1可用第1时段水位测量下降高度为968－822＝146来检验，类似地，y2用1082－822＝260来检验。 供水时段流量的一种检验方法如下：供水时段用水量加上水位上升值260是该时段泵入的水量，除以时间长度得到水泵的功率（单位时间泵入水量），而两个供水时段的功率应大致相等。第1、2时段水泵的功率计算如下： >> p1=(y12+260)/2 p1 = 155.1995 >> tp2=20.8:0.1:23; >> xp2=polyval(c3,tp2); >> p2=(0.1*trapz(x3)+260)/2.2 p2 = 152.2335 可以看到，两次水泵泵水的功率差别不大。下面是水塔一天的流量曲线图： 图14.1 当取三次多项式拟合的流量曲线图 由图14.1我们可以看到，流量曲线与原始记录基本上相吻合，但在第1时段和第1泵水时段的交接处曲线不太光滑，这说明我们采用3次曲线通过4点的做法不够好，应该多取几点进行拟合。0点到10点很流量很低，10点到下午3点即中午时间段是用水高峰期。【练习与思考】 1．假定某天的气温变化见下表，试找出这一天的气温变化规律： 2．在化工生产中常常需要知道丙烷在各种温度T和压力P下的导热系数K。下面是实验得到的一组数据： 试求T＝99（℃）和P＝10.3（1032 kN）下的导热系数。 /m 3．下表给出了某一海域以码为单位的直角坐标为X、Y的水面一点处以英尺为单位的水深Z，水深数据是在低潮时测得的。船的吃水深度为5英尺，问在矩形区域（85200）×（－40150）里，哪些地方船要避免进入。