估计水塔的水流量.doc

估计水塔的水流量
美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率
以及每天所用的总水量.许多社区没有测量流入或流出当地水塔的水量的装置,他
们只能代之以每小时测量水塔中的水位,其精度不超过5%,更重要的是,当水塔中
的水位下降最低水位L 时水泵就启动向水塔输水直到最高水位H,但也不能测量
水泵的供水量.因此,当水泵正在输水时不容易建立水塔中水位和水泵工作时用水
量之间的关系.水泵每两天输水一次或两次,每次约二小时.
试估计任何时刻(包括水泵正在输水的017921 时间内)从水塔流出的流量f(t),并估计
一天的总用水量.附表给出了某各小镇一天中真实的数据.
附表给出了从第一次测量开始的以秒为单位的时刻.以及该时刻的高度单位
为百分之一英尺的水位测量值.例如,3316 秒后,水塔中水位达到31.10 英尺.水塔
是一个高为40 英尺,直径为57 英尺的正圆柱.通常当水塔水位降至约27.00 英尺
的水泵开始工作,当水位升到35.50 英尺时水泵停止工作.

时间（秒）
水位(0.01英尺)
时间（秒）
水位(0.01英尺)
时间（秒）
水位(0.01英尺)
0
3175
35932
水泵工作
68535
2842
3316
3110
39332
水泵工作
71854
2767
6635
3054
39435
3550
75021
2697
10619
2994
43318
3445
79154
水泵工作
13937
2947
46636
3350
82649
水泵工作
17921
2892
49953
3260
85968
3475
21240
2850
53936
3167
89953
3397
25223
2797
57254
3087
93270
3340
28543
2752
60574
3012

32284
2697
64554
2927

问题分析与数据处理
由问题的要求,关键在于确定用水率函数,即单位时间内用水体积,记为f(t),又称水流速度.如果能够通过测量数据,产生若干个时刻的用水率,也就是f(t)在若干个点的函数值,则f(t)的计算问题就可以转化为插值或拟合问题

一,问题假设
1)水塔中水流量是时间的连续光滑函数,与水泵工作与否无关,并忽略水位高度对水流速度的影响.
2)水泵工作与否完全取决于水塔内水位的高度,且每次加水的工作时间为2小时.
3)水塔为标准圆柱体.
4)水泵第一次供水时间为[32284, 39435],第二次供水时间段为[75021,85948].
5)为了方便计算我们把表格中的秒转化成小时.
6）我们规定以下符号：
h:水塔中水位的高度，是时间的函数，单位为英尺；
v:水塔中水的体积，是时间的函数，单位为加仑；
t:时间，单位为小时；
f:模型估计的水塔水流量，是时间的函数，单位为加仑/小时
p:水泵工作时的充水水流量，也是时间的函数，单位为加仑/小时。

二.体积计算
水塔是一个圆柱体,体积为v=pi*D^2*h,其中D=57英尺 得到不同时刻水塔中水的体积如下
水塔中水的体积
时间
（小时）
水体积
（加仑）
时间
（小时）
水体积
（加仑）
时间
（小时）
水体积
（加仑）
0
0．9211
1．8431
2．9497
3．8714
4．9781
5．9000
7．0064
7．9286
8．9678
606125
593716
583026
571571
562599
552099
544081
533963
525372
514872
9．98101
10.9256
10．9542
12．0328
12．9544
13．8758
14．9822
15．9039
16．8261
17．9317

水泵工作
水泵工作
677715
657670
639534
622352
604598
589325
575008
558781

19．0375
19．9594
20．8392
22．015．
22．9581
23．8800
24．9869
25．9．83

542554
528236
514872
水泵工作
水泵工作
663397
648506
637625

三,.水流速度的估算
水流速度应该是水塔中水的体积对时间的导数(微商).由于没有水的体积关于时间的函数表达式,而只有一个离散的函数值表,因此考虑用差商代替微商,
我们已经得到了水塔中水的体积Vi与时间ti的关系，由于f(t)=||，要得到f(t)曲线，我们用差分公式得到fi~ti关系，当然由于没有水泵工作期间的Vi，由差分得到的fi也没有水泵工作期间的数据。
可以看出，Vi很像均匀隔开的数据点，并因水泵工作而分成三组，我们这样来处理这些数据：对每一组数据点，为了减少误差，我们采用中心差分公式：
fi:=f(ti)=|| 及
Fi:=f(ti)=||,
Fi:=f:(ti)=| | 得到水流量值的表格如下:
时间
（小时）
水体积
（加仑）
时间
（小时）
水体积
（加仑）
时
（小时）
水体积
（加仑）
0
0．9211
1．8431
2．9497
3．8714
4．9781
5．9000
7．0064
7．9286
8．9678
14405
11180
10063
11012
8797
9992
8124
10160
8488
11018

9．9811
109256
10．9542
12．0328
12．9544
13．8758
14．9822
15．9039
16．8261
17．9317

水泵工作
水泵工作
19469
20196
18941
15903
18055
15646
13741
14962

19．0375
19．9594
20．8392
22．0150
22．9581
23．8800
24．9869
25．9083

16653
14496
14648
水泵工作
15225
15264
13708
9633
用SPSS大致画出函数图象如下：

四, 用三次样条拟合fi
上面所得到的fi是相当粗糙的，而且还不包括水泵工作期间的数据。我们采用三次样条函数插值表数据，得到光滑的水流量曲线f(t)。但有一个问题：水泵工作期间的水流量如何拟合？根据假设，水流量只依靠于公众对水的需求，是一种自然的规律，它本身是一条相当光滑的曲线，有水泵工作时的数据当然最好，现在不知道，我们只能依据连续性，领先充水前后的数据来拟合曲线，为了得到水泵工作时的水流量，我们忽略水泵工作期间的数据，直接对充水前后的数据用三次样条插值来拟合。
PS：三次插值样条函数
定义 设在区间［a,b］上给定一个分割∏：a=x0在每个节点xi(i=1,2,…,n-1)处，
S(k)(xi-0)=S(k)(xi+0),k=0,1,2 (2)
则称S(x)为三次样条函数．
对定义在区间［a,b］上的函数f(x),如果存在三次样条函数S(x)，使得在节点处还满足S(xi)=f(xi)(i=0,1, …，n),就称S(x)为插值于f(x)的三次样条函数．
对给定的一组有序数组yi(i=0,1, …，n),如果三次样条函数S(x)满足S(xi)=yi(i=0,1, …，n),就称S(x)为插值于｛yi｝的三次样条函数．
现在，如果对函数f(x)，我们并不知道其解析表达式，而只知道其在节点处的值fi=f(xi) (i=0,1, …，n),如何估计f(x)?一个很自然的方法就是求插值于｛fi｝的三次样条函数S(x),以S(x)作为对f(x)的逼近．那么，如何求出S(x)?我们将利用fi及一阶、二阶导数来建立求S(x)的表示式及连续性方程．
(１)Ｍ连续性方程与S(x)的表示式
记S(x)在节点xi处的函数值、一阶导数和二阶导数分别为
S(xi)=fi,S′(xi)=mi,,S"(xi)=Mi, (i=0,1, …，n) (3)
由于S(x)是分片三次多项式，在每个小区间[xi-1,xi]上，S(x)的二阶导数是线性函数，记hi=xi-xi-1表示小区间长度，有
S〃(x)=Mi-1, (xi-1≤x≤xi) (4)
将（4）式积分一次，得
S＇(x)=-Mi-1 , (xi-1≤x≤xi) (5)
再将（5）式积分一次，有
S(x)=Mi-1 (xi-1≤x≤xi) (6)
由插值条件（3），S（xi）=fi,S(xi-1)=fi-1,代入（6）式，有

而由（5）式，有
(7)
但由一阶导数连续，S＇(xi-0)=S’(xi+0)(i=1,…,n-1),由（7）式就得到n-1个等式
μiMi-1+2Mi+λiMi+1=di, (i=1,…,n-1) (8)
其中
λi=,μi=
di= (i1,…,n-1)。 (9)
方程组（8）称为S(x)的M连续性方程。
我们来看（6）式，假定fi(i=0,1,…,n)都已知，如果再知道M0,…,Mn,则（6）式就决定了S(x)的具体表达式。
实际上方程组（8）具有和有直观的力学意义，那就是λi与μi表示相邻区间[xi-1,xi]与[xi,xi+1]的长度比，而di为插值数据在xi处的二阶中心差商的3倍，那么，（8）式说明了插值函数的二阶导数在xi-1,xi,xi+1三点处的加权平均（权因子分别为μi/3,2/3,λi/3）为被插数据在xi处的二阶中心差商，这就是力学上的“三弯矩”关系。

设有三次样条函数S(t)通过表的数据点，具有端点to=0和t24=25.9083，设在每一个小区间[ti,ti+1]上，S(t)是一个三次方程，为
Si(t)=a0i+a1i(t-ti)+a2i(t-ti)2+a3i(t-ti)3
在每一个节点ti处，Si-1(t)与Si(t)的函数值、一阶导数和二阶导数相等，都反映了水流量曲线f(t)的光滑性质，这时就有限制方程
fi= a0i+a1i(t-ti)+a2i(t-ti)2+a3i(t-ti)3
= a0I-1+a1I-1(ti-ti-1)+a2I-1(ti-ti-1)2+a3I-1(ti-ti-1)3,
f’i=a1I+2a2I(t-ti)+3a3I(t-ti)2
=a1I-1+2a2I-1(ti-ti-1)+3a3I-1(ti-ti-12)
f”i=2a2I+6a3I(t-ti)=2a2I-1+6a3I-1(ti-ti-1)
但我们利用M连续性方程来求解。首先为了得到完整的连续性方程，设在端点t0和t24处被拟合曲线f(t)的斜率为

此时有hi,λi,μi和di(I=0,1,…,24)，求解Mi(I=0,1,,24),进而得到S(t)的表示式，算得qi,pi,从而得到Mi最后，有
（*）
把fi与Mi(i=0,1,…,24)代入（*）就有Si(t)的系数，这样就得到了-S(t)的分段表示式，我们把S(t)作为对水流量的拟合曲线f(t)，
表 λiμi和di的值
i
hi
λi
μI
di
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

0.9211
0.9220
1.1066
0.9217
1.1067
0.9219
1.1064
0.9222
1.0392
1.9864
1.0786
0.9216
0.9214
1.1064
0.9217
0. 9222
1.1056
1.1058
0.9219
0.8798
2.1189
0.9219
1.1069
0.9214
1
0.5002
0.5455
0.4544
0.5456
0.4544
0.5455
0.4546
0.5298
0.6565
0.3519
0.4607
0.4999
0.5456
0.4545
0.5001
0.5452
0.5000
0.4546
0.4883
0.7066
0.3032
0.5456
0.4543
0
0
0.4998
0.4545
0.5456
0.4544
0.5456
0.4545
0.5454
0.4702
0.3435
0.6481
0.5393
0.5001
0.4544
0.5455
0.4999
0.4548
0.5000
0.5454
0.5117
0.2934
0.6968
0.4544
0.5457
1
7453.78
7454.02
6119.69
-9645.76
10302.58
-9186.75
11437.51
-10805.26
12993.63
3608.92
-7008.96
-6106.69
-6300.49
15511.03
-13496.61
1782.98
9379.89
1152.65
-11448.26
8371.74
199.17
-453.85
-4282.39
-8924.42
-11164.92

表 qi,pi与Mi的值
I
qi
pI
Mi
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
-0.5000
-0.2858
-0.2917
-0.2468
-0.2890
-0.2467
-0.2890
-0.2468
-0.2812
0.3449
-0.1981
-0.2434
-0.2662
-0.2904
-0.2468
-0.2665
-0.2902
-0.2696
-0.2453
-0.2605
-0.3673
-0.1739
-0.2840
-0.2462
0
3726.89
3194.86
2495.92
-5901.07
6948.20
-6957.01
7812.85
-8076.74
8989.61
276.90
-3840.76
-2159.41
-2796.42
8983.59
-9857.11
3594.72
4146.11
-493.46
-5990.27
6126.83
-844.74
71.73
-2306.31
-4108.96
-3527.98
2832.16
1789.46
4917.46
-8301.38
9725.71
-9610.77
10757.05
-10187.53
8552.68
1553.81
-3702.27
-699.11
-599.61
12033.09
-10501.04
2609.13
3698.26
1543.23
-7554.50
6376.81
-959.62
312.76
-1386.04
-3240.37
-3527.98

表Si(x)的系数a0i,a1i,a2i,a3i的值
i
区间[ti,ti+1]
a0i,
a1
a2i
A3i
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
[0,0.9211]
[0. 9211,1.8431]
[1.8431,2.9497]
[2.9497,3.8714]
[3.8714,4.9781]
[4.9781,5.9000]
[5.9000,7.0064]
[7.0064,7.9286]
[7.9286,8.9678]
[8.9678,10.9542]
[10.9542,12.0328]
[12.0328,12.9544]
[12.9544,13.8758]
[13.8758,14.9822]
[14.9822,15.9039]
[15.9039,16.8261]
[16.8261,17.9317]
[17.9317,19.0375]
[19.0375,19.9594]
[19.9594,20.8392]
[20.8392,22.8581]
[22.9581,23.8800]
[23.8800,24.9869]
[24.9869,25.9083]
14405
11180
10063
11012
8797
9992

8124
10160
8488
11018
19469
20196
18941
15903
18055
15646
13741
14962
16653
14296
14648
15225
15264
13708
-4646
-2517
575
-1347
-735
-726
-248
4
-797
4451
2131
-225
-3302
-556
212
-3436
-543
2353
-998
-1557
840
159
-297
-2886
1416
895
2459
-4151
4863
-4805
5379
-5094
4276
777
-1851
-350
-3000
6017
-5251
1305
1849
772
-3777
3188
-480
156
-693
-1620
-189
565
-1991
3260
-2912
3682
-3155
3387
-1122
-441
464
-959
3262
-3395
2371
197
-325
-1371
2519
-1390
100
-307
-279
-52

因为样条曲线S(t)的定义域区间并起来覆盖了整个区间[t0,t27],这样拟合的f(t)就对任何时刻t都有了定义。
PS：由三次样条插值计算得到用水率函数f(t)
MATLAB程序:
x0=t;y0=r;
[1,n]=size(x0); d1=x0(n)-x0(1);
x=x0(1):1/3600:x0(n); %被插值点
ys=interp1(x0,y0,x,`spline`); %样条插值输出
plot(x,ys);
xlabel(`时间(小时)`)；ylabel(`流速(立方米/小时)`）
title(`样条插值下的流速图`）
五’对水泵充水的两段时间的处理
水泵充水的水流量，单从我们已知的数据是无法准确得到的。但我们可以分析一下，由假设，水泵抽水与水塔的水流量是无关,其充水水流量只与水泵本身的性能有关，至多与水塔的高度有关，而与水位也无关，我们可以估计一下水泵充水的平均水流量，即使水泵充水水流量不是常数，它也不应该是变化不定的，对第一次充水，水塔的水体积之差为：△V1=677715-514872=162843加仑，
充水时间：
△ t1=10.9542-8.9678=1.9864小时
则第一次充水水泵平均水流量为

=81979+15597
=97576（加仑/小时）。
对第二次充水，水塔的水体积之差△V2=△V1，充水时间△t2=22.9581-20.8392=2.1189小时，所以第二次充水水泵平均水流量为

=76853+15057
=91910（加仑/小时）。
两次充水的平均流量就为
p=(p1+p2)=94743（加仑/小时）.
得一天总用水量
我们已得到水流量曲线f(t),要计算一天的总用水量就简单了：对f(t)在24小时区间积分即可，利用f(t)的表达式，对[0，24]区间进行积分，有
(加仑).
即该社区一天的总用水量约为333189加仑。