水塔流量估计

1、水塔流量估计水塔流量估计摘摘要要 .1关键字关键字 .2一一.问题重述问题重述 .2二二.问题分析问题分析 .2三三.模型假设模型假设 .3四问题求解四问题求解 .3（1）数据差值法 .6（2）数值拟合法 .11五五.结果简单分析结果简单分析 .18六模型检验六模型检验 .18七模型改进七模型改进 .18八八.实验小结实验小结 .18九参考文献九参考文献 .18附录附录 .181摘要针对水塔流量估计问题，由于水塔流量与水位变化密切相关，所以通过水位变化来确定其流量；用流速和时间的关系来确定水位的变化，即水位差；由物理学中 Torricelli 定律确定流速，继而解决具体问题。我们对某一天水塔水

2、位统计数据进行分析，将其分为水泵供水时段和水泵不工作时段，然后通过插值与拟合模型分别进行计算。计算结果如下：图 1.1 插值法算得各时段句全天的用水量（用高度计）第一未供水时段第二未供水时段第一供水时段混合时段全天拉格朗日插值法145.6231258.866454.0689127.8244586.3828分段线性插值法147.1430258.969749.6051115.0578570.7756三次样条插值法145.6870258.654753.3334120.0178577.6929图 1.2 数据拟合法算得各个时段及全天的用水量（用高度计）时段第一未供水时段第二未供水时段第一供水时段混合时

3、段全天用水高度146.5150257.760546.131776.3076526.7148关键字：流量 插值与拟合 torricelli 定律 一.问题重述:某一地区的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用水的总量。许多社区没有测量流入或流出水塔的水量装置，他们只能通过测量水塔每小时的水位高度来代替，其误差不超过 0.05.但是，当水塔中的水位下降到最低水位 L 时水泵就自动启动向水塔输水直到最高水位 H，在此期间无法测量水泵的供水量。这样，当水泵正在输水时就不容易建立水塔中水位和用水量之间的关系，水泵每天输水一次或两次，每次约两小时。当水位降至约 27.00ft（英

4、尺）时，水泵开始启动向水塔供水，当水位升到35.50ft 时，水泵自动停止工作。已知某一小镇的水塔高度为 40ft，直径为 57ft 的正圆柱，下表记录的是某一天的水塔水位的真实数据（注：1ft=0.3024m） 。表 1.1 某一天水塔水位统计表时间/s水位/0.01ft时间/s水位/0.01ft203175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275

5、021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开动859683475394353550899533397433183445932703340试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量和一天的用水总量。二.问题分析：流量是单位时间内流出水的体积，由于水塔是正圆柱形，横截面积时常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化率算出。问题的难点在于如何估计水泵供水时段的流量。水泵供水时段的流量只能依据供水时段前后的流量经差值或拟合得到。作为用于差值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。事实上，水泵不工

6、作时段的用水量可以由测量记录直接得到吗，由表中记录的下降水位乘以水塔的横截面积就该时段的用水量。这个数值可以用来检验数据差值或拟合的结果。为了表示方便，我们将表中的数据全部化为国际标准单位（如下表） ，时间用小时（h） ，高度用米（m）。表 2.1 一天内水塔水位记录时间（h）水位（m）时间（h）水位（m）09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.8

7、48.228.978.2222.02水泵开动39.98水泵开动22.96水泵开动10.93水泵开动23.8810.5910.9510.8224.9910.3512.0310.5025.9110.18三.模型假设（1）假设流量只取决于水位差，与水位本身无关。由水塔的最低和最高水位分别是 8.1648m（27x0.3024）和 10.7352m（35.50 x0.3024） ，设出口的水位为零，从小孔流出的流速正比于水面高度的平方根，可忽略水位对流速的影响。（2）假设流量是时间的连续函数。为计算简单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位时下降的） ，得到结果后在乘

8、以水塔横截面积即可。（3）已知水塔横截面积为 S=（57*0.3048）2*/4=237.8（平方米） 。四问题求解首先根据表 2.1 的数据，用 MATLAB 软件作出水位时间散点图 4.1。程序见附录 1图 4.1 水位时间散点图下面计算流量与时间的关系。根据数据散点图 1.0,一种简单的处理方法是先将表 2.0 中的数据分为三段，然后对每一段数据做如下处理：设某段数句为（x0，y0） ， （x1,y1）,（xn，yn），相邻数据中点的平均流速采用公式：流速=（左端点的水位-右端点的水位）/区间长度，即：v（）=(yi-y(i+1)/(x（i+1）-xi)2) 1(xiix4每段数据首尾点

9、的柳树采用下面的公式计算：v(x0)=(3y0-4y1+y2)/(x2-x0),v(xn)=(-3yn+4y(n-1)-y(n-2).用以上公式算得流速与时间之间的数据表 4.0 如下：表 4.1 流速与时间关系数据表时间（h）流速（cm/h）时间（h）流速（cm/h）029.8913.4229.030.4621.7414.4326.361.3818.4815.4426.092.39516.2216.3724.733.4116.3017.3823.644.42515.3218.4923.425.4413.0419.5025.006.4515.4520.4023.867.46513.9820.8

10、422.178.4516.3522.02水泵开动8.9719.2922.96水泵开动9.98水泵开动23.8827.0910.93水泵开动24.4321.6210.9533.5025.4518.4811.4929.6325.9113.3012.4931.52由表 4.1 作出流速时间散点图 4.2，程序见附录 2：图 4.2 流速时间关系散点图5下面用数据差值和数据拟合两种方法来估计水塔水流量：（1）数据差值法 由表 4.1，对水泵不工作时段 1 和水泵不工作时段 2 采用差值方法。可以得到任意时刻的流速，从而可以知道任意时刻的流量。对于水泵工作时段数据太少，我们将它与水泵工作时段 2 合并一

11、同进行差值。由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段 2 合并一同进行差值处理（该时段一下简称混合时段） 。这样，总共需要对四段数据（第 1，2 未供水时段，第 1供水时段，混合时段）进行差值处理。A.第 1 未供水时段数据差值求解，分别用拉格朗日插值法，分段线性差值及三次样条差值三种方法算出流量函数和用水量（用水高度） 。图 4.3 中曲线 lglr，fdxx，scyt 分别表示用拉格朗日差值法，分段线性差值法及是三次样条差值法对第 1 未供水时段数据差值得到的流速曲线，程序见附录 3图 4.3 第 1 未供水时段流速曲线B. 第 2 未供水时段数据差值求解，分别用拉格朗日

12、插值法，分段线性差值及三次样条差值三种方法算出流量函数和用水量（用水高度） 。图 4.4 中曲线 lglr，fdxx，scyt 分别表示用拉格朗日差值法，分段线性差值法及是那次样条差值法对第 2 未供水时段数据差值得到的流速曲线， 程序代码见附录 4：6图 4.4 第 2 未供水时段流速曲线C.第 1 供水时段数据差值求解，分别用拉格朗日插值法，分段线性差值及三次样条差值三种方法算出流量函数和用水量（用水高度） 。图 4.5 中曲线 lglr，fdxx，scyt 分别表示用拉格朗日差值法，分段线性差值法及是那次样条差值法对第 2 未供水时段数据差值得到的流速曲线，程序代码见附录 5：图 4.5

13、 第 2 未供水时段流速曲线7D.混合时段数据差值求解，分别用拉格朗日插值法，分段线性差值及三次样条差值三种方法算出流量函数和用水量（用水高度） 。图 4.2 中曲线 lglr，fdxx，scyt 分别表示用拉格朗日差值法，分段线性差值法及是那次样条差值法对第 2 未供水时段数据差值得到的流速曲线:程序执行后显示结果，程序代码见附录 6：图 4.6 混合时段流速曲线综合得到：表 4.2 差值法算得个时段及全天的用水总量（以高度计）时段及方法第 1 未供水时段第 2 未供水时段第 1 供水时段混合时段全天拉格朗日插值法145.6231258.866454.0689127.8244586.3828

14、分段线性差值法147.1430258.969749.6051115.0578570.7756三次样条差值法145.6870258.654753.3334120.0178577.6929图 4.6 是用分段线性及三次样条差值方法算得全天水塔水流速与时间的关系曲线，图中曲线 scyt 表示三次样条差值曲线，曲线 fdxx 表示分段线性差值的曲线，程序代码见附录 7：8图 4.7 全天水塔流速与时间的关系曲线(2)数值拟合法1）拟合水位时间函数根据表 4.1 中的数据，分别对第 1、2 未供水时段的测量数据直接做三次多项式拟合，得到水位与实践的关系函数。程序中，变量 x1 存放了以 0.1 为时间步

15、长第 1 未供水时段各个时刻的水位高度。图 4.7 绘出的是第 1 未供水时段水位与时间的三次多项式拟合曲线，程序代码见附录 8：图 4.8 第 1 未供水时段水位与时间的三次多项式拟合曲线9下图4.9，Matlab程序对第2未供水时段的数据进行拟合,得到的第2未供水时段水位与时间三次多项式拟合曲线，Matlab程序见附录9.图 4.9 第 2 未供水时段水位与时间的三次多项式拟合曲线2）确定流量时间函数对第 1,2 未供水时段的水位求导可得流量，用三次多项式拟合第 1,2 未供水时段的流速与时间关系曲线，程序代码见附录 10：图 4.10 第 1,2 未供水时段的流速与时间关系曲线（三次多项

16、式拟合）10与图 4.6 中相应部分段对比发现，用三次多项式拟合效果不是很好，若改用 5 次多项式拟合则得效果较好的第 1,2 未供水时段的流速与时间关系曲线，见图 4.9(程序代码见附录 11):图 4.11 第 1,2 未供水时段的流速与时间关系曲线（五次多项式拟合）第 1 供水时段的流速则用前后时刻的流速拟合得到。为使流速函数在 t=9和 t=11 处连续，只取 4 个点，用 3 次多项式拟合得第 1 供水时段流速与时间关系曲线图 4.12，同图 4.1 中相应部分段较为吻合，程序代码见附录 12：图 4.12 第一供水时段流速与时间的三次多项式似合曲线11对于混合时段，在第 2 供水时

17、段之前取 t=19.96,20.84 两点的流速，用第三末供水时段的 3 个记录做差分得到两个流速数据 21.62,18.48，然后用这 4 个数据做 3 次多项式似合得到混合时段的流速与时间的关系曲线如图 4.13 所示，同数据插值法得到的曲线图 4.1 中相应曲线段较为吻合，程序代码见附录 13：图 4. 11 混合时段流速与时间的三次多项式拟合曲线3）估计一天的用水总量分别对供水的两个时段和不供水的两个时段积分（流量对时间积分）并求和得到一天的用水总量约为 526.8935（此数据是用水总高度 cm）.表 4.3.1 列出各段用水量，同数据插值法算得的表 4.1 数据相比，较为吻合。表

18、4.3 数据拟合算得各个时段的用水高度（一）时段第 1 未供水时段第 2 未供水时段第 1 供水时段混合时段全天用水高度145.67260.6646.6073.9635526.89354）流量及用水总量的检验计算出各时刻的流量可以用水位记录数据的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。例如，算得第 1 未供水时段的用水高度是 145.67，而实际记录的水位下降高度为 968-822=146，两者是吻合的同样地，算得第 2 未供水时段的用水高度是 260.66，而实际记录的水位下降高度为 1082-822=260，两者也是吻合的。从算法设计和分析可知，计算结果与各时段

19、所用的拟合多项式的次数有关。表 4.4 给出的是对第 1，2 未供水时段分别用 5，6 次多项式拟合得到的用水高度。12表 4.4 数据拟合法算得各个时段的用水高度（二）时段第 1 未供水时段第 2 未供水时段第 1 供水时段混合时段全天用水高度146.5150257.760546.131776.3076526.7148五.结果简单分析对于数据差值法，由表 4.1 可以看出，使用三次样条差值法得到的第 1,2 未供水时段的用水高度结果 145.6870 和 258.6547 与表中记录的下降高度 146 和260 相差不大，说明差值结果与原始数据比较吻合。用三次样条差值法估计全天的用水量约为

20、539.4450 x237.8x10=1282800.21（升） 。对于数据拟合法，由表 4.3.2 可得全天的用水总量约为526.7148x237.8x10=1252527.7944（升） ，同数据差值法得到的结果也很接近。六模型检验七模型改进八.实验小结本次实验主要进行水塔水量的估算。第一种估计为数据差值方法，我们用了三种不同的差值方法进行估计。在求解的过程中，可以熟悉数据差值的理论和方法。第二种估计方法为数据拟合法， 用多项式进行拟合，得到水塔水流量的估计。九参考文献附录：下面程序是对上面问题的补充，均由 matlab 软件编写：附录附录 1：附录附录 2：x=xlsread(x=xls

21、read(stgj.xlsstgj.xls, ,sheet1sheet1, ,A2:A29A2:A29););y=xlsread(y=xlsread(stgj.xlsstgj.xls, ,sheet1sheet1, ,B2:B29B2:B29););plot(x,y,plot(x,y,b*b*) )xlabel(xlabel(xx轴轴 );ylabel();ylabel(yy轴轴 ););x=xlsread(x=xlsread(stgjls.xlsstgjls.xls, ,sheet1sheet1, ,A2:A32A2:A32););y=xlsread(y=xlsread(stgjls.xls

22、stgjls.xls, ,sheet1sheet1, ,B2:B32B2:B32););plot(x,y,plot(x,y,b*b*) )xlabel(xlabel(xx轴轴 );ylabel();ylabel(yy轴轴 ););13附录补充：拉格朗日函数文件附录附录 3：functionunction y=lglrcz(x0,y0,x)y=lglrcz(x0,y0,x)n=length(x0);n=length(x0);m=length(x);m=length(x);forfor i=1:mi=1:m z=x(i);z=x(i); s=0.0;s=0.0; forfor k=1:nk=1:n

23、 p=1.0;p=1.0; forfor j=1:nj=1:n ifif j=kj=k p=p*(z-x0(j)/(x0(k)-x0(j);p=p*(z-x0(j)/(x0(k)-x0(j); endend endend s=p*y0(k)+s;s=p*y0(k)+s; endend y(i)=s;y(i)=s;endendt=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.9

24、7;lglr=lglrcz(t,v,t0);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)14附录附录 4t=10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.4t=10.95,11.49

25、,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84;0,20.84;v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.8v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17;6,22.17;t0=10.95:0.1:20.84;t0=10.95:0.1:20.84;lglr=lglrcz(t,v,t0);lglr=lglrcz(t,v,t0);lglrj

26、f=0.1*trapz(lglr)lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,spline);scyt=interp1(t,v,t0,spline);sancytjf=0.1*trapz(scyt)sancytjf=0.1*trapz(scyt)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)plot(t,v,*,t0,lglr,r,t0,fdxx,g

27、,t0,scyt,b)gtext(lglr)gtext(lglr)gtext(fdxx)gtext(fdxx)gtext(scyt)gtext(scyt)附录附录 5：t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84;49,12.49,13.42,14.43,

28、15.44,16.37,17.38,18.49,19.50,20.40,20.84;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.2v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.869,33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.

29、00,23.86,22.17;,22.17;t0=8.97:0.1:10.95;t0=8.97:0.1:10.95;lglr=lglrcz(t,v,t0);lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,spline);scyt=interp1(t,v,t0,spline);sancytjf=0.1*tra

30、pz(scyt)sancytjf=0.1*trapz(scyt)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(lglr)gtext(fdxx)gtext(fdxx)gtext(scyt)gtext(scyt)15附录附录 6：t=10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20t=10.95,11.49,12.49,13.42,14.43,15.44,16.

31、37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91;.84,23.88,24.43,25.45,25.91;v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22v=33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.62,18.48,13.30;.17,27.09,21.62,18.48,13.30;t0=20.84:0.1:22.96;

32、t0=20.84:0.1:22.96;lglr=lglrcz(t,v,t0);lglr=lglrcz(t,v,t0);lglrjf=0.1*trapz(lglr)lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,spline);scyt=interp1(t,v,t0,spline);sancytjf=0.1*trapz(scyt)sancytjf=0.1*trapz(scyt)plot(t,

33、v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(lglr)gtext(fdxx)gtext(scyt)附录附录 7：t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,1

34、8.49,19.50,20.40,20.84,23.88,24,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91;.43,25.45,25.91;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63,31.52,29.03,26.36,26

35、.09,24.73,23.64,23.42,25.00,23.86,22.33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.62,18.48,13.30;17,27.09,21.62,18.48,13.30;t0=0:0.1:25.91;t0=0:0.1:25.91;fdxx=interp1(t,v,t0);fdxx=interp1(t,v,t0);fdxxjf=trapz(fdxx)fdxxjf=trapz(fdxx)scyt=interp1(t,v,t0,scyt=interp1

36、(t,v,t0,splinespline););sancytjf=trapz(scyt)sancytjf=trapz(scyt)plot(t,v,plot(t,v,*,t0,fdxx,t0,fdxx,gg,t0,scyt,t0,scyt,bb) )gtext(gtext(fdxxfdxx) )gtext(gtext(scytscyt) )附录附录 8：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.t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.9

37、5,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.6688,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66h=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,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

38、,8.92,8.66,8.43,8.22,10.59,10.35,10.18;9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),3);c1=polyfit(t(1:10),h(1:10),3); % %输入项输入项 3 3 表示拟合多项式的次数表示拟合多项式的次数tp1=0:0.1:8.9;tp1=0:0.1:8.9;x1=polyval(c1,tp1);x1=polyval(c1,tp1);plot(tp1,x1)plot(tp1,x1)16附录附录 9：t=0,0.92,

39、1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,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.6613.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66h=9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.3

40、9,8.22,10.82,10.50,10.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;21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18;c1=polyfit(t(1:10),h(1:10),5);c1=polyfit(t(1:10),h(1:10),5); % %输入项输入项 5 5 表示拟合多项式的次数表示拟合

41、多项式的次数tp1=0:0.1:8.9;tp1=0:0.1:8.9;x1=polyval(c1,tp1);x1=polyval(c1,tp1);plot(tp1,x1)plot(tp1,x1)附录附录 10：c1=polyfit(t(1:10),h(1:10),3);c1=polyfit(t(1:10),h(1:10),3); % %用三次多项式拟合用三次多项式拟合c2=polyfit(t(11:21),h(11:21),3);c2=polyfit(t(11:21),h(11:21),3);a1=polyder(c1);a1=polyder(c1);a2=polyder(c2);a2=poly

42、der(c2);tp1=0:0.01:8.97;tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x13=-polyval(a1,tp1);x113=-polyval(a1,0:0.01:8.97);x113=-polyval(a1,0:0.01:8.97);wgsysl1=100*trapz(tp1,x113);wgsysl1=100*trapz(tp1,x113); % %计算第计算第 1 1 未供水时段的总用水量未供水时段的总用水量x14=-polyval(a1,7.93,8.97

43、);x14=-polyval(a1,7.93,8.97); % %为后面的程序准备数据为后面的程序准备数据x23=-polyval(a2,tp2);x23=-polyval(a2,tp2);x114=-polyval(a2,10.95:0.01:20.84);x114=-polyval(a2,10.95:0.01:20.84);wgsysl2=100*trapz(tp2,x114);wgsysl2=100*trapz(tp2,x114); % %计算第计算第 2 2 未供水时段的总用水量未供水时段的总用水量x24=-polyval(a2,10.95,12.03);x24=-polyval(a2

44、,10.95,12.03); % %为后面的程序准备数据为后面的程序准备数据x25=-polyval(a2,19.96,20.84);x25=-polyval(a2,19.96,20.84); % %为后面的程序准备数据为后面的程序准备数据subplot(1,2,1)subplot(1,2,1)plot(tp1,x13*100)plot(tp1,x13*100)subplot(1,2,2)subplot(1,2,2)plot(tp2,x23*100)plot(tp2,x23*100)17附录附录 11：c1=polyfit(t(1:10),h(1:10),3);c1=polyfit(t(1:1

45、0),h(1:10),3); % %用五次多项式拟合用五次多项式拟合c2=polyfit(t(11:21),h(11:21),3);c2=polyfit(t(11:21),h(11:21),3);a1=polyder(c1);a1=polyder(c1);a2=polyder(c2);a2=polyder(c2);tp1=0:0.01:8.97;tp1=0:0.01:8.97;tp2=10.95:0.01:20.84;tp2=10.95:0.01:20.84;x13=-polyval(a1,tp1);x13=-polyval(a1,tp1);x113=-polyval(a1,0:0.01:8.97);x113=-polyval(a1,0:0.01:8.97);wgsysl1=100*trapz(tp1,x113);wgsysl1=100*trapz(tp1,x113); % %计算第计算第 1 1 未供水时段的总用水量未供水时段的总用水量x14=