估计水塔的水流量

估计水塔的水流量美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天所用的总水量.许多社区没有测量流入或流出当地水塔的水量的装置,他们只能代之以每小时测量水塔中的水位,其精度不超过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 ，由差分得到的fi 也没有水泵工作期间的数据。

可以看出，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 2i 1i i -+-+++| 得到水流量值的表格如下:F r e q u e n c y四, 用三次样条拟合f i上面所得到的f i 是相当粗糙的，而且还不包括水泵工作期间的数据。

我们采用三次样条函数插值表数据，得到光滑的水流量曲线f(t)。

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

PS ：三次插值样条函数定义 设在区间［a,b ］上给定一个分割∏：a=x 0<x 1<…<x n-1<x n =b,定义在［a,b ］上的一个函数S(x)如果满足下列条件：①在每个小区间［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-1i1i ii i h x x M h x x --+-, (x i-1≤x ≤x i ) (4)将（4）式积分一次，得S ＇(x)=-M i-1i 1i21i 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 1i1i i i ++-+-- (x i-1≤x ≤x i )(6)由插值条件（3），S （x i ）=f i ,S(x i-1)=f i-1,代入（6）式，有⎪⎪⎩⎪⎪⎨⎧+-+-=---=-----1i 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 ii 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 =i1i ih h h ++di=)h f f h f f (h h 6i1i 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 处的二阶中心差商，这就是力学上的“三弯矩”关系。

设有三次样条函数S(t)通过表的数据点，具有端点t o =0和t 24=25.9083，设在每一个小区间[t i ,t i+1]上，S(t)是一个三次方程，为 S i (t)=a 0i +a 1i (t-t i )+a 2i (t-t i )2+a 3i (t-t i )3在每一个节点t i 处，S i-1(t)与S i (t)的函数值、一阶导数和二阶导数相等，都反映了水流量曲线f(t)的光滑性质，这时就有限制方程 fi= a 0i +a 1i (t-t i )+a 2i (t-t i )2+a 3i (t-t i )3 =a 0I-1+a 1I-1(t i -t i-1)+a 2I-1(t i -t i-1)2+a 3I-1(t i -t i-1)3, f’i=a 1I +2a 2I (t-t i )+3a 3I (t-t i )2=a 1I-1+2a 2I-1(t i -t i-1)+3a 3I-1(t i -t i-12) f”i=2a 2I +6a 3I (t-t i )=2a 2I-1+6a 3I-1(t i -t i-1)但我们利用M 连续性方程来求解。

首先为了得到完整的连续性方程，设在端点t 0和t 24处被拟合曲线f(t)的斜率为⎪⎪⎩⎪⎪⎨⎧-≈-+-=-+-=-≈--+-=-++-=56.5789)9869.249083.25(2152645483228899)t t (2f f 4f 3m 53.4645)09211.0(2100634472043215)t t (2f f 4f 3m 232422232424012100此时有h i ,λi ,μi 和d i (I=0,1,…,24)，求解M i (I=0,1,,24),进而得到S(t)的表示式，算得q i ,p i ,从而得到M i 最后，有⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧--==+----==++++++)t t (6M M a M 21a 6)M M 2)(t t (t t f f a f a i 1i i 1i i 3ii 21i i i 1i i 1i i 1i i 1i i 0 （*）把f i 与M i (i=0,1,…,24)代入（*）就有S i (t)的系数，这样就得到了-S(t)的分段表示式，我们把S(t)作为对水流量的拟合曲线f(t)， 表 λ和d 的值表Si(x)的系数a的值因为样条曲线S(t)的定义域区间并起来覆盖了整个区间[t0,t27],这样拟合的f(t)就对任何时刻t 都有了定义。