估计水塔的水流量

估计水塔的水流量

美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率

以及每天所用的总水量.许多社区没有测量流入或流出当地水塔的水量的装置,他们只能代之以每小时测量水塔中的水位,其精度不超过5%,更重要的是,当水塔中的水位下降最低水位L 时水泵就启动向水塔输水直到最高水位H,但也不能测量水泵的供水量.因此,当水泵正在输水时不容易建立水塔中水位和水泵工作时用水量之间的关系.水泵每两天输水一次或两次,每次约二小时.

试估计任何时刻(包括水泵正在输水的017921 时间内)从水塔流出的流量f(t),并估计

一天的总用水量.附表给出了某各小镇一天中真实的数据.

附表给出了从第一次测量开始的以秒为单位的时刻.以及该时刻的高度单位

为百分之一英尺的水位测量值.例如,3316 秒后,水塔中水位达到31.10 英尺.水塔是一个高为40 英尺,直径为57 英尺的正圆柱.通常当水塔水位降至约27.00 英尺的水泵开始工作,当水位升到35.50 英尺时水泵停止工作.

问题分析与数据处理

由问题的要求,关键在于确定用水率函数,即单位时间内用水体积,记为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英尺 得到不同时刻水塔中水的体积如下

三,.水流速度的估算

水流速度应该是水塔中水的体积对时间的导数(微商).由于没有水的体积关于时间的函数表达式,而只有一个离散的函数值表,因此考虑用差商代替微商, 我们已经得到了水塔中水的体积V i 与时间t i 的关系，由于f(t)=|

dt

)t (dV |，要得到f(t)曲线，我

们用差分公式得到f i ~t i 关系，当然由于没有水泵工作期间的V i ，由差分得到的f

i 也没有水泵工作期间的数据。

可以看出，V i 很像均匀隔开的数据点，并因水泵工作而分成三组，我们这样来处理这些数据：对每一组数据点，为了减少误差，我们采用中心差分公式： fi:=f(ti)=|)

t t (12V V 8V 8V i 1i 2i 1i 1i 2i -+-+-+--++| 及

F i :=f(ti)=|)

t t (2V V 4V 3i 1i 2i 1i i --+-+++|,

F i :=f:(ti)=|

)

t t (2V V 4V 3i 1i 2

i 1i i -+-+++| 得到水流量值的表格如下:

F r e q u e n c y

四, 用三次样条拟合f i

上面所得到的f i 是相当粗糙的，而且还不包括水泵工作期间的数据。我们采用三次样条函数插值表数据，得到光滑的水流量曲线f(t)。但有一个问题：水泵工作期间的水流量如何拟合？根据假设，水流量只依靠于公众对水的需求，是一种自然的规律，它本身是一条相当光滑的曲线，有水泵工作时的数据当然最好，现在不知道，我们只能依据连续性，领先充水前后的数据来拟合曲线，为了得到水泵工作时的水流量，我们忽略水泵工作期间的数据，直接对充水前后的数据用三次样条插值来拟合。

PS ：三次插值样条函数

定义 设在区间［a,b ］上给定一个分割∏：a=x 0

①在每个小区间［x i-1,x i ］(ｉ=1,2, …，n)内S(x)是三次多项式；

②在整个区间［a,b ］上，S(x)为二阶连续可导函数，也就是说，在每个节点x i (i=1,2,…,n-1)处， S (k)

(x i -0)=S (k)

(x i +0),k=0,1,2

(2)

则称S(x)为三次样条函数．

对定义在区间［a,b ］上的函数f(x),如果存在三次样条函数S(x)，使得在节点处还满足S(x i )=f(x i )(i=0,1, …，n),就称S(x)为插值于f(x)的三次样条函数．

对给定的一组有序数组y i (i=0,1, …，n),如果三次样条函数S(x)满足S(x i )=y i (i=0,1, …，n),就称S(x)为插值于｛y i ｝的三次样条函数．

现在，如果对函数f(x)，我们并不知道其解析表达式，而只知道其在节点处的值f i =f(x i ) (i=0,1, …，n),如何估计f(x)?一个很自然的方法就是求插值于｛f i ｝的三次样条函数S(x),以S(x)作为对f(x)的逼近．那么，如何求出S(x)?我们将利用f i 及一阶、二阶导数来建立求S(x)的表示式及连续性方程． (１)Ｍ连续性方程与S(x)的表示式

记S(x)在节点x i 处的函数值、一阶导数和二阶导数分别为 S(x i )=f i ,S ′(x i )=m i ,,S"(x i )=M i , (i=0,1, …，n) (3)

由于S(x)是分片三次多项式，在每个小区间[x i-1,x i ]上，S(x)的二阶导数是线性函数，记h i =x i -x i-1表示小

区间长度，有 S 〃(x)=M i-1

i

1

i i

i i h x x M h x x --+-, (x i-1≤x ≤x i ) (4)

将（4）式积分一次，得

S ＇(x)=-M i-1

i 1i

2

1i i 2i C h )x x (Mi h 2)x x (+-+-- , (x i-1≤x ≤x i ) (5)

再将（5）式积分一次，有 S(x)=M i-1

,C x C h 63

)x x (Mi h 63)x x (i 2i 1i

1i i i ++-+-- (x i-1≤x ≤x i )

(6)

由插值条件（3），S （x i ）=f i ,S(x i-1)=f i-1,代入（6）式，有

⎪⎪⎩

⎪

⎪⎨⎧

+-+-=---=-----1

i i i i i i 1i i i 1i i 21i i i 1i f i 1x )6M h h f (x )6M h h f (C 6)M M (h hi f C i 而由（5）式，有

⎪⎪⎩

⎪

⎪⎨⎧----=++---=-+++++--2h M 6)M M (h h f f )0x ('S 2h M 6)M M (h hi f f )0x ('S 1i i

i 1i 1i 1i i 1i i i i 1i i i 1i i i (7)

但由一阶导数连续，S ＇(x i -0)=S’(x i +0)(i=1,…,n-1),由（7）式就得到n-1个等式 μi M i-1+2M i +λi M i+1=d i , (i=1,…,n-1) (8)

其中 λi=1i i 1i h h h +++,μi =i

1i i

h h h ++

di=

)h f f h f f (h h 6

i

1i i 1i i 1i 1i i -+++---+ (i=1,…,n-1)。

(9)

方程组（8）称为S(x)的M 连续性方程。

我们来看（6）式，假定f i (i=0,1,…,n)都已知，如果再知道M 0,…,M n ,则（6）式就决定了S(x)的具体表达式。 实际上方程组（8）具有和有直观的力学意义，那就是λi 与μi 表示相邻区间[x i-1,x i ]与[x i ,x i+1]的长度比，而d i 为插值数据在x i 处的二阶中心差商的3倍，那么，（8）式说明了插值函数的二阶导数在x i-1,x i ,x i+1三点处的加权平均（权因子分别为μi /3,2/3,λi /3）为被插数据在x i 处的二阶中心差商，这就是力学上的“三弯矩”