第五章 水塔用水量的估计----插值
§5—1 水塔用水量问题本题中所使用的长度单位为E(=30.24cm)；容积单位为G(=3.785L(升)).
某些州的用水管理机构需估计公众的用水速度(单位是G/h)和每天总用水量的数据。许多地方没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位(误差不超过5%),当水箱水位低于某最低水位L时，水泵抽水，灌入水箱内直至水位达到某最高水位H为止。但是也无法测量水泵的流量，因此在水泵启动时无法立即将水箱中的水位和水量联系起来。水泵一天灌水1~2次，每次约2h,试估计在任意时刻(包括水泵灌水期间)t流出水箱的流量f(t)，并估计一天的总用水量。
下表给出了某镇中某一天的真实用水数据，表中测量时间以秒为单位，水位以0.01E为单位。例如在3316s时，水箱中的水深为31.10E,
时间(s)
水位(0.01E)
时间(s)
水位(0.01E)
0
3175
46636
3350
3316
3110
49953
3260
6635
3054
53936
3167
10619
2994
57254
3087
13937
2947
60574
3012
17921
2892
64554
2927
21240
2850
68535
2842
25223
2795
71854
2767
28543
2752
75021
2697
32284
2697
79254
泵水
35932
泵水
82649
泵水
39332
泵水
85968
3475
39435
3550
89953
3397
43318
3445
93270
3340
已知水箱是直径为57E、高为40E的正圆柱体。当水位落到27E以下时，水泵自动启动把水泵入水箱；当水位回升至35.5E时，水泵停止工作。
§5—3 水塔用水量的计算
5.3.1 问题分析
1．假设
28次采集数据，得到24对数据(时刻：秒。水位：英尺)，要估计水流出水箱的速度f(t)(即：任意时刻t时的水流量)。先作几条假设。
(1),只考虑通常情况下的正常用水，不考虑突发事变，如水管破裂、自然灾害等。
(2),全天24小时能保证用户的正常用水，不会发生某个用户想用水时打开水龙头而没水流出这样的情况。即：水箱中水的流出量f(t)完全由用户的需求量决定，与水箱自身的情况如水位高度、出口口径等因素无关。
(3),题中所给“水位约27.00E时，水泵启动；水位约35.50时，水泵停止；每次工作约2小时”这样的已知信息说的比较含糊，故这里将其作为一条假设，同时也假设水泵的灌水速度为常数、水泵不会损坏或不需要维修。
(4),水箱的出水流量f(t)，前面已假设是由用户决定的。单个用户的用水流量是时间t的突变函数(即：突然很大，突然为0)，当大量的用户随机用水时(本题的实际背景正是这样)，总的用水流量函数f(t)就可以近似当作是关于时间t连续的、光滑的，故假设水流量函数f(t)一阶连续可导。
本题目标是计算出水箱的出水流量函数f(t)、及一天24小时总的用水量。已知条件为24个时刻的水位，我们的想法是：分析确定水泵的两次启动、停止工作的时刻，用数值计算的方法确定测量点时刻的水流量，以这些测量点作为节点做三次样条插值，得出水流量函数f(t),
2,水泵的工作情况为了便于计算，将时间单位变成小时h，长度单位变成英尺E，见下表：
时间(h)
水位(E)
时间(h)
水位(E)
0
31.75
12.954
33.50
0.921
31.10
13.876
32.60
1.843
30.54
14.982
31.67
2.950
29.94
15.904
30.87
3.871
29.47
16.826
30.12
4.978
28.92
17.932
29.27
5.900
28.50
19.038
28.42
7.006
27.95
19.959
27.67
7.929
27.52
20.839
26.97
8.968
26.97
22.015
泵水
9.981
泵水
22.958
泵水
10.926
泵水
23.880
34.75
10.954
35.50
24.987
33.97
12.033
34.45
25.908
33.40
为了方便分析，从上述28对信息中提取出一小部分来：
时间(h)
水位(E)
时间(h)
水位(E)
8.968
26.97
20.839
26.97
9.981
泵水
22.015
泵水
10.926
泵水
22.958
泵水
10.954
35.50
23.880
34.75
根据已知，水位降低到约27E时，水泵开始启动。在t=8.968h时刻，水位是26.97E，低于27E，水泵应该工作了，但实际情况是此时刻水泵没工作。显然，过了t=8.968h后很短的时间内，水泵将开始工作，故近似确定，t=8.968h是水泵第一次工作的启动时刻。同理，近似确定t=20.839h是水泵第二次工作的启动时刻。
根据已知，水位上升到约35.5E时，水泵停止工作。在t=10.954h时刻，水位是35.50E，且水泵已停止工作。再考虑到已知水泵每次工作约2小时，而10.954-8.968=1.986<2，即水泵这次实际工作的时间长度不到2小时，所以就近似确定t=10.954h是水泵第一次工作的停止时刻。
设水泵第二次工作的停止时刻为T，则22.958<T<23.880,且T时刻的水位大约是35.5E,我们依据以下两个理由推测出h,
理由(1)：已知水泵每次工作约2小时，启动时刻为20.839，且22.958-20.839=2.119>2，实际情况是，在22.958时刻水泵还再工作，本次工作时长已经超过2小时了，故可推断在此刻之后的很短时间内水泵就要停止工作了。
理由(2)：（考察邻近时间段的用水情况） 从19.959至 20.839时段，对应水位从 27.67降到 26.97，平均降速为0.796E/h；从23.88至 24.987时段，对应水位从 34.75降到 33.97，平均降速为0.705E/h；实际中T>22.958，从T至 23.88时段，对应水位从 35.5降到 34.75，平均降速为
E/h
注意到T离22.958越远，水位降速会越比0.813大，而0.813已明显比0.796、0.705大，该时段内的水位降速不应该再大了，故推断T离22.958很近。
基于上述两个理由，我们近似确定t=22.958h是水泵第二次工作的停止时刻，且此刻水位是35.50E,(注：原来已知24对数据，现在又多了一对已知数据)
3．测量点时刻的水流量
t时刻的水流量为f(t)，设水箱中水的体积为V(t)，则，当水泵不工作的时候必有
,
把25个测量点时刻从小到大排列，记为 ,
时刻的水位记为
时刻水箱中水的体积记为，则有
时刻的水流量记为,
根据已知数据，我们用“数值微分”的方法来计算：
数学知识：数值微分问题1：y=f(x)连续可导，不知到表达式，
已知 ,
求 y=f(x) 在每个节点处的一阶导数。
分析：学过插值法（如：Lagrange插值、三次样条插值）和最小二乘法的人，很容易想到用这样的方法：先拟合出y=f(x)的近似表达式，再求导函数，再计算导函数在每个节点处的值。此法不好，应该直接用已知数据计算节点处的一阶导数。
数值微分法：节点处的一阶导数记为
i=1,2时,
时：
时,,
问题2：y=f(x)表达式已知，且可导，但求不出导函数，要求任意一点处的一阶导数值。
数值微分法：取小步长h
补充知识完毕.
将25个节点分成3个组第一组
第二组
第三组
每组的最前两点，用公式：
每组的最末两点，用公式：
每组的其余点，用公式：
5.3.2 模型建立与求解
前面已经得到测量点时刻的水流量
用三次样条插值可得水流量函数
对积分，可得任何一个时间段的用水总量。
下面程序计算并画出图象（程序中的数据t0,v0取自书P67之表5.4）：
clear
t0=[0,0.921,1.843,2.949,3.871,4.978,5.9,7.006,7.928,8.967,10.954,12.032,12.954,13.875,14.982,15.903,16.826,17.931,19.037,19.959,20.839,22.958,23.88,24.986,25.908];
v0=[2294,2247,2206,2163,2129,2089,2059,2020,1988,1948,2564,2489,2420,2355,2288,2230,2176,2114,2053,1999,1948,2564,2510,2454,2413];
for i=1:25
switch i
case {1,2,11,12,22,23}
f0(i)=abs((-3*v0(i)+4*v0(i+1)-v0(i+2))/(2*(t0(i+1)-t0(i))));
case {9,10,20,21,24,25}
f0(i)=abs((3*v0(i)-4*v0(i-1)+v0(i-2))/(2*(t0(i)-t0(i-1))));
otherwise
f0(i)=abs((-v0(i+2)+8*v0(i+1)-8*v0(i-1)+v0(i-2))/(12*(t0(i+1)-t0(i))));
end
end
t=0:0.05:26;f=spline(t0,f0,t);plot(t,f),grid
上面程序中，时间t以小步长0.05遍历计算得到了每个点处的f值，用“数值积分中的复化梯形法”计算区间[0,24]上定积分就得一天总用水量。程序如下：
（续上面程序内容）
zysl=0;
for j=1:480
zysl=zysl+0.05*(f(j)+f(j+1))/2;
end
zysl
计算结果为：zysl= 1.2582e+003
5.3.3 模型的检验
从几个侧面，分别对模型进行检验：
1．以不同时刻作为起点，计算总用水量
下面程序实现书P70表5.6
（续上面程序内容）
for qd=0.25:0.25:1.75
j0=qd/0.05;
zysl=0;
for j=j0+1:j0+480
zysl=zysl+0.05*(f(j)+f(j+1))/2;
end
qd,zysl
end
计算结果为：
qd =0.2500 zysl =1.2593e+003
qd =0.5000 zysl =1.2608e+003
qd =0.7500 zysl =1.2627e+003
qd =1 zysl =1.2649e+003
qd =1.2500 zysl =1.2671e+003
qd =1.5000 zysl =1.2688e+003
qd =1.7500 zysl =1.2698e+003
2．分三段，实际用水量与模型计算量比较。
略算，看书P70表5.7
3．估计水泵的功率略算，书中得：
水泵第一次工作：充水量=733.446，工作时长=1.987
水泵第一次工作：充水量=736.705，工作时长=2.119
第一次功率
第一次功率
两次功率的相对误差为
§5—3 二维插值的应用
（黑板：简略介绍“二维插值问题”）
已知节点及节点处的函数值，求二元函数，并且画曲面图形。
关键：interp2(…………)
下面程序实现书P72之图：
clear
yf=1:12;wd=5:10:85;
qx=[2.4,18.7,20.8,22.1,37.3,48.2,25.6,5.3,0.3;1.6,21.4,18.5,20.1,28.8,36.6,24.2,5.3,0;
2.4,16.2,18.2,20.5,27.8,35.5,25.5,5.4,0;3.2,9.2,16.6,25.1,37.2,40,24.6,4.9,0.3;
1.0,2.8,12.9,29.2,40.3,37.6,21.1,4.9,0;0.5,1.7,10.1,32.6,41.7,35.4,22.2,7.1,0;
0.4,1.4,8.3,33.0,46.2,35,20.2,5.3,0.1;0.2,2.4,11.2,31.0,39.9,34.7,21.2,7.3,0.2;
0.5,5.8,12.5,28.6,35.9,35.7,22.6,7,0.3;0.8,9.2,21.1,32.0,40.3,39.5,28.5,8.6,0;
2.4,10.3,23.9,28.1,38.2,40,25.3,6.3,0.1;3.6,16,25.5,25.6,43.4,41.9,24.3,6.6,0.3];
[x,y]=meshgrid(1:1:12,5:5:85);
z=interp2(yf,wd,qx',x,y,'cubic');
mesh(x,y,z)
§5—1 水塔用水量问题本题中所使用的长度单位为E(=30.24cm)；容积单位为G(=3.785L(升)).
某些州的用水管理机构需估计公众的用水速度(单位是G/h)和每天总用水量的数据。许多地方没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位(误差不超过5%),当水箱水位低于某最低水位L时，水泵抽水，灌入水箱内直至水位达到某最高水位H为止。但是也无法测量水泵的流量，因此在水泵启动时无法立即将水箱中的水位和水量联系起来。水泵一天灌水1~2次，每次约2h,试估计在任意时刻(包括水泵灌水期间)t流出水箱的流量f(t)，并估计一天的总用水量。
下表给出了某镇中某一天的真实用水数据，表中测量时间以秒为单位，水位以0.01E为单位。例如在3316s时，水箱中的水深为31.10E,
时间(s)
水位(0.01E)
时间(s)
水位(0.01E)
0
3175
46636
3350
3316
3110
49953
3260
6635
3054
53936
3167
10619
2994
57254
3087
13937
2947
60574
3012
17921
2892
64554
2927
21240
2850
68535
2842
25223
2795
71854
2767
28543
2752
75021
2697
32284
2697
79254
泵水
35932
泵水
82649
泵水
39332
泵水
85968
3475
39435
3550
89953
3397
43318
3445
93270
3340
已知水箱是直径为57E、高为40E的正圆柱体。当水位落到27E以下时，水泵自动启动把水泵入水箱；当水位回升至35.5E时，水泵停止工作。
§5—3 水塔用水量的计算
5.3.1 问题分析
1．假设
28次采集数据，得到24对数据(时刻：秒。水位：英尺)，要估计水流出水箱的速度f(t)(即：任意时刻t时的水流量)。先作几条假设。
(1),只考虑通常情况下的正常用水，不考虑突发事变，如水管破裂、自然灾害等。
(2),全天24小时能保证用户的正常用水，不会发生某个用户想用水时打开水龙头而没水流出这样的情况。即：水箱中水的流出量f(t)完全由用户的需求量决定，与水箱自身的情况如水位高度、出口口径等因素无关。
(3),题中所给“水位约27.00E时，水泵启动；水位约35.50时，水泵停止；每次工作约2小时”这样的已知信息说的比较含糊，故这里将其作为一条假设，同时也假设水泵的灌水速度为常数、水泵不会损坏或不需要维修。
(4),水箱的出水流量f(t)，前面已假设是由用户决定的。单个用户的用水流量是时间t的突变函数(即：突然很大，突然为0)，当大量的用户随机用水时(本题的实际背景正是这样)，总的用水流量函数f(t)就可以近似当作是关于时间t连续的、光滑的，故假设水流量函数f(t)一阶连续可导。
本题目标是计算出水箱的出水流量函数f(t)、及一天24小时总的用水量。已知条件为24个时刻的水位，我们的想法是：分析确定水泵的两次启动、停止工作的时刻，用数值计算的方法确定测量点时刻的水流量，以这些测量点作为节点做三次样条插值，得出水流量函数f(t),
2,水泵的工作情况为了便于计算，将时间单位变成小时h，长度单位变成英尺E，见下表：
时间(h)
水位(E)
时间(h)
水位(E)
0
31.75
12.954
33.50
0.921
31.10
13.876
32.60
1.843
30.54
14.982
31.67
2.950
29.94
15.904
30.87
3.871
29.47
16.826
30.12
4.978
28.92
17.932
29.27
5.900
28.50
19.038
28.42
7.006
27.95
19.959
27.67
7.929
27.52
20.839
26.97
8.968
26.97
22.015
泵水
9.981
泵水
22.958
泵水
10.926
泵水
23.880
34.75
10.954
35.50
24.987
33.97
12.033
34.45
25.908
33.40
为了方便分析，从上述28对信息中提取出一小部分来：
时间(h)
水位(E)
时间(h)
水位(E)
8.968
26.97
20.839
26.97
9.981
泵水
22.015
泵水
10.926
泵水
22.958
泵水
10.954
35.50
23.880
34.75
根据已知，水位降低到约27E时，水泵开始启动。在t=8.968h时刻，水位是26.97E，低于27E，水泵应该工作了，但实际情况是此时刻水泵没工作。显然，过了t=8.968h后很短的时间内，水泵将开始工作，故近似确定，t=8.968h是水泵第一次工作的启动时刻。同理，近似确定t=20.839h是水泵第二次工作的启动时刻。
根据已知，水位上升到约35.5E时，水泵停止工作。在t=10.954h时刻，水位是35.50E，且水泵已停止工作。再考虑到已知水泵每次工作约2小时，而10.954-8.968=1.986<2，即水泵这次实际工作的时间长度不到2小时，所以就近似确定t=10.954h是水泵第一次工作的停止时刻。
设水泵第二次工作的停止时刻为T，则22.958<T<23.880,且T时刻的水位大约是35.5E,我们依据以下两个理由推测出h,
理由(1)：已知水泵每次工作约2小时，启动时刻为20.839，且22.958-20.839=2.119>2，实际情况是，在22.958时刻水泵还再工作，本次工作时长已经超过2小时了，故可推断在此刻之后的很短时间内水泵就要停止工作了。
理由(2)：（考察邻近时间段的用水情况） 从19.959至 20.839时段，对应水位从 27.67降到 26.97，平均降速为0.796E/h；从23.88至 24.987时段，对应水位从 34.75降到 33.97，平均降速为0.705E/h；实际中T>22.958，从T至 23.88时段，对应水位从 35.5降到 34.75，平均降速为
E/h
注意到T离22.958越远，水位降速会越比0.813大，而0.813已明显比0.796、0.705大，该时段内的水位降速不应该再大了，故推断T离22.958很近。
基于上述两个理由，我们近似确定t=22.958h是水泵第二次工作的停止时刻，且此刻水位是35.50E,(注：原来已知24对数据，现在又多了一对已知数据)
3．测量点时刻的水流量
t时刻的水流量为f(t)，设水箱中水的体积为V(t)，则，当水泵不工作的时候必有
,
把25个测量点时刻从小到大排列，记为 ,
时刻的水位记为
时刻水箱中水的体积记为，则有
时刻的水流量记为,
根据已知数据，我们用“数值微分”的方法来计算：
数学知识：数值微分问题1：y=f(x)连续可导，不知到表达式，
已知 ,
求 y=f(x) 在每个节点处的一阶导数。
分析：学过插值法（如：Lagrange插值、三次样条插值）和最小二乘法的人，很容易想到用这样的方法：先拟合出y=f(x)的近似表达式，再求导函数，再计算导函数在每个节点处的值。此法不好，应该直接用已知数据计算节点处的一阶导数。
数值微分法：节点处的一阶导数记为
i=1,2时,
时：
时,,
问题2：y=f(x)表达式已知，且可导，但求不出导函数，要求任意一点处的一阶导数值。
数值微分法：取小步长h
补充知识完毕.
将25个节点分成3个组第一组
第二组
第三组
每组的最前两点，用公式：
每组的最末两点，用公式：
每组的其余点，用公式：
5.3.2 模型建立与求解
前面已经得到测量点时刻的水流量
用三次样条插值可得水流量函数
对积分，可得任何一个时间段的用水总量。
下面程序计算并画出图象（程序中的数据t0,v0取自书P67之表5.4）：
clear
t0=[0,0.921,1.843,2.949,3.871,4.978,5.9,7.006,7.928,8.967,10.954,12.032,12.954,13.875,14.982,15.903,16.826,17.931,19.037,19.959,20.839,22.958,23.88,24.986,25.908];
v0=[2294,2247,2206,2163,2129,2089,2059,2020,1988,1948,2564,2489,2420,2355,2288,2230,2176,2114,2053,1999,1948,2564,2510,2454,2413];
for i=1:25
switch i
case {1,2,11,12,22,23}
f0(i)=abs((-3*v0(i)+4*v0(i+1)-v0(i+2))/(2*(t0(i+1)-t0(i))));
case {9,10,20,21,24,25}
f0(i)=abs((3*v0(i)-4*v0(i-1)+v0(i-2))/(2*(t0(i)-t0(i-1))));
otherwise
f0(i)=abs((-v0(i+2)+8*v0(i+1)-8*v0(i-1)+v0(i-2))/(12*(t0(i+1)-t0(i))));
end
end
t=0:0.05:26;f=spline(t0,f0,t);plot(t,f),grid
上面程序中，时间t以小步长0.05遍历计算得到了每个点处的f值，用“数值积分中的复化梯形法”计算区间[0,24]上定积分就得一天总用水量。程序如下：
（续上面程序内容）
zysl=0;
for j=1:480
zysl=zysl+0.05*(f(j)+f(j+1))/2;
end
zysl
计算结果为：zysl= 1.2582e+003
5.3.3 模型的检验
从几个侧面，分别对模型进行检验：
1．以不同时刻作为起点，计算总用水量
下面程序实现书P70表5.6
（续上面程序内容）
for qd=0.25:0.25:1.75
j0=qd/0.05;
zysl=0;
for j=j0+1:j0+480
zysl=zysl+0.05*(f(j)+f(j+1))/2;
end
qd,zysl
end
计算结果为：
qd =0.2500 zysl =1.2593e+003
qd =0.5000 zysl =1.2608e+003
qd =0.7500 zysl =1.2627e+003
qd =1 zysl =1.2649e+003
qd =1.2500 zysl =1.2671e+003
qd =1.5000 zysl =1.2688e+003
qd =1.7500 zysl =1.2698e+003
2．分三段，实际用水量与模型计算量比较。
略算，看书P70表5.7
3．估计水泵的功率略算，书中得：
水泵第一次工作：充水量=733.446，工作时长=1.987
水泵第一次工作：充水量=736.705，工作时长=2.119
第一次功率
第一次功率
两次功率的相对误差为
§5—3 二维插值的应用
（黑板：简略介绍“二维插值问题”）
已知节点及节点处的函数值，求二元函数，并且画曲面图形。
关键：interp2(…………)
下面程序实现书P72之图：
clear
yf=1:12;wd=5:10:85;
qx=[2.4,18.7,20.8,22.1,37.3,48.2,25.6,5.3,0.3;1.6,21.4,18.5,20.1,28.8,36.6,24.2,5.3,0;
2.4,16.2,18.2,20.5,27.8,35.5,25.5,5.4,0;3.2,9.2,16.6,25.1,37.2,40,24.6,4.9,0.3;
1.0,2.8,12.9,29.2,40.3,37.6,21.1,4.9,0;0.5,1.7,10.1,32.6,41.7,35.4,22.2,7.1,0;
0.4,1.4,8.3,33.0,46.2,35,20.2,5.3,0.1;0.2,2.4,11.2,31.0,39.9,34.7,21.2,7.3,0.2;
0.5,5.8,12.5,28.6,35.9,35.7,22.6,7,0.3;0.8,9.2,21.1,32.0,40.3,39.5,28.5,8.6,0;
2.4,10.3,23.9,28.1,38.2,40,25.3,6.3,0.1;3.6,16,25.5,25.6,43.4,41.9,24.3,6.6,0.3];
[x,y]=meshgrid(1:1:12,5:5:85);
z=interp2(yf,wd,qx',x,y,'cubic');
mesh(x,y,z)
