插值和拟合在水流流量计算中的运用Word文档下载推荐.doc

插值和拟合在水流流量计算中的运用Word文档下载推荐.doc

《插值和拟合在水流流量计算中的运用Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《插值和拟合在水流流量计算中的运用Word文档下载推荐.doc（22页珍藏版）》请在冰点文库上搜索。

1.2分段线性插值法

假设区间[a,b]的连续函数g（x）在n+1个节点上的函数值。

则得到xy平面上的n+1个数据点。

连接相邻数据点，得到n条线段，它们组成一条折线。

把区间[a,b]上这n条折线段表示的函数称为被插值函数g（x）关于这n+1个数据点的分段线性插值函数，记作I（x）具有如下性质：

1）I（x）可以分段表示，在每个小区间上，它是线性函数，即

1）

2）

3）在[a,b]上连续。

若构造插值函数

因则

当g（x）在[a,b]上连续时，分段线性插值函数I（x）具有良好的收敛性，即而且当g（x）在[a,b]上二阶导数连续时，对于任意有

其中用计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n+1无关。

但n越大，分段越多，插值误差越小。

实际上用数据点作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。

+

1.3三次样条插值

分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。

例如：

在船体、飞机等外形曲线的设计中，不仅要求曲线连续而且还要求曲线的曲率连续，这就要求插值函数具有连续的二阶导数。

为解决这一类问题，就产生了三次样条插值。

从数学上加以概括，可得到样条函数的定义如下：

三次样条函数记作，，满足：

①在每个小区间是三次多项式。

②在每个内节点上具有二次连续导数。

③

由三次样条函数中的条件①知，有个待定系数。

由条件②知，在个内节点上具有二阶连续导数，即满足条件：

共有个条件。

由条件③，知，共有个条件。

因此，要确定一个三次样条，还需要外加个条件，最常用的三次样条函数的边界条件有两类：

第一类边界条件：

第二类边界条件：

特别地，,称为自然边界条件。

第三类边界条件：

称为周期边界条件。

三次样条插值不仅光滑性好，而且稳定性和收敛性都有保证，具有良好的逼近性质。

样条插值函数的建立。

构造满足条件的三次样条插值函数的表达式可以有多种方法。

下面我们利用的二阶导数值表达，由于在区间上是三次多项式，故在上是线性函数，可表示为（5）

其中对积分两次并利用及，可定出积分常数，于是得三次样条表达式

（6）

上式中是未知的，为确定，对求导得

（7）

由此可得。

同样求出在区间上的表达式，从而得

利用可得（8）

其中；

；

（9）

对第一类边界条件，可导出两个方程

（10）

如果令，则式（8）及其（10）可写出矩阵（11）

通过求解上述三对角矩阵可求得。

对于第二类边界条件，直接得端点方程（12）

如果令，则式（8）及式（12）也可以写成矩阵（11）的形式。

对于第三类边界条件，可得（13）

其中，

则式（8）及式（13）可以写成矩阵形式

求解上述矩阵可得。

1.4一维插值法总结

拉格朗日插值函数在整个插值区间上有统一的解析表达式，其形式关于节点对称，光滑性好。

但缺点同样明显，这主要体现在高次插值收敛性差（龙格现象）；

增加节点时前期计算作废，导致计算量大；

一个节点函数值的微小变化（观测误差存在）将导致整个区间上插值函数都发生改变，因而稳定性差等几个方面。

因此拉格朗日插值法多用于理论分析，在采用拉格朗日插值方法进行插值计算时通常选取n<

7。

分段线性插值函数（仅连续）与三次样条插值函数（二阶导数连续）虽然光滑性差，但他们都克服了拉格朗日插值函数的缺点，不仅收敛性、稳定性强，而且方法简单实用，计算量小。

因而应用十分广泛。

分段线性插值，具有良好的稳定性和收敛性，但光滑性较差。

在数学上若函数（曲线）的k阶导数存在且连续，则称该曲线具有k阶光滑性。

易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。

二.曲线拟合

已知一组数据，即平面上的n个点，互不相同，寻求一个函数（曲线）y=f（x），使f（x）在某种准则下与所有数据点最为接近，即曲线拟合得最好。

2.1最小二乘法拟合法

设为个线性无关的函数，对给定的数据，

求使最小。

利用极值的必要条件。

得到关于的线性方程组

则方程组可表示为，

其中，，

由于线性无关，所以是列满秩，是可逆矩阵，方程组的解存在且唯一，并且。

取，

得多项式拟合方程。

三．数值实验

为了更准确明白插值分析与拟合方法的实用性，采用以下示例进行测试.

3.1实例测试

许多社区没有测量流入或流出水塔的水量装置，他们只能代之以每小时测量水塔中的水位，其误差不超过0.5%。

更重要的是，当水塔中的水位下降到最低水位L时水泵就启动向水塔输水，直到最高水位H，但也不能测量水泵的供水量。

因此，当水泵在输水时不容易建立水塔中水位和水泵工作时用水量之间的关系。

水泵每天输水一次或两次，每次约二小时。

试估计任意时刻（包括水泵在输水工作的时候）从水塔流出的流量f（t），并估计一天的总水量。

表1给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺的水位测量值。

表1

时间/h

水位/m

9.68

12.95

10.21

0.92

9.45

13.88

9.94

1.84

9.31

14.98

9.65

2.95

9.13

15.90

9.41

3.87

8.98

16.83

9.18

4.98

8.81

17.93

8.92

5.90

8.69

19.04

8.66

7.00

8.52

19.96

8.43

7.93

8.39

20.84

8.22

8.97

22.02

水泵开动

9.98

22.96

10.93

23.88

10.59

10.95

10.82

24.99

10.35

12.03

10.50

25.91

10.18

（1）影响水箱流量的唯一因素是该区公众对水的普通需要；

（2）水泵的灌水速度为常数；

（3）从水箱中流出水的最大流速小于水泵的灌水速度；

（4）每天的用水量分布都是相似的；

（5）水箱的流水速度可用光滑曲线来近似；

（6）当水箱的水容量达到514×

103g时，开始泵水；

达到677.6×

103g时，便停止泵水

3.3实例解答

首先依照表1所给数据，用MATLAB作出时间—水位散点图（图1）。

下面来计算水箱流量与时间的关系。

根据图1一种简单的处理方法为将表1中的数据分为三段，然后对每一段的数据做如下处理：

设某段数据，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位－左端点的水位）/区间长度）：

每段数据首尾点的流速用下面的公式计算：

用以上公式求得时间与流速之间的数据如表2。

表2

流速/cm·

h-1

流速/cm·

29.89

12.49

31.52

0.46

21.74

13.42

29.03

1.38

18.48

14.43

26.36

2.395

16.22

15.44

26.09

3.41

16.30

16.37

24.73

4.425

15.32

17.38

23.64

5.44

13.04

18.49

23.42

6.45

15.45

19.50

25.00

7.465

13.98

20.40

23.86

8.45

16.35

22.17

19.29

27.09

33.50

24.43

21.62

11.49

29.63

25.45

13.30

由表2作出时间—流速散点图如图2。

3.3.1插值法分析

由表2，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。

我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；

对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它忽略，只对水泵工作时段2进行插值处理。

我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。

调用附件程序1实现图3结果

图3

运行结果：

lglrjf=145.6250fdxxjf=147.1469sancytjf=145.6885

图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。

由表1知，第1未供水时段的总用水高度为146（=968－822），可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。

考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。

其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图像.

第一供水段时间—流速示意图（程序代码见附件程序2）如下

图4

lglrjf=56.4426fdxxjf=49.6051sancytjf=53.5903

第2未供水段时间—流速示意图（程序代码见附件程序3）如下

图5

lglrjf=258.8664fdxxjf=258.9697sancytjf=258.6547

第2供水段时间—流速示意图（程序代码见附件程序4）如下

图6

lglrjf=104.1526 fdxxjf=73.7929 sancytjf=79.8172

下图是用分段线性及三次样条插值方法得到的整个过程的时间—流速函数示意图。

（程序代码见附件程序5）

图7

fdxxjf=534.4311sancytjf=541.8148

表3-1各时段及一天的总用水量（用水高度（cm））

第1未供水段

第2未供水段

第1供水段

第2供水段

第3未供水段

全天（除3未供）

拉格朗日插值法

145.625

258.8664

56.4426

104.1526

40.7825

565.0866

分段线性插值法

147.1465

258.9697

49.6051

73.7929

40.7516

529.5142

三次样条插值法

145.6885

258.6547

53.5903

79.8172

537.7507

表3-2各时段及一天的总用水量（用水体积（m3）

339.8123

604.0583

131.7074

243.0375

95.1649

1318.616

343.3627

604.2993

115.7523

172.1939

95.0928

1235.608

339.9605

603.5643

125.0516

186.2514

1254.828

表4是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值.

时间（h）

6.88

10.88

15.88

22.88

①

15.9826

33.7426

25.5662

34.7099

②

14.8272

32.9976

25.4465

25.4715

③

15.0527

33.7089

25.5490

29.4173

注：

①拉格朗日插值法②分段线性插值法③三次样条插值法

3.3.2拟合法分析

（1）拟合水位—时间函数

从表1中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以3～6次为宜）。

对第3未供水时段来说，数据过少不能得到很好的拟合。

设t，h分别为已输入的时刻和水位测量记录（由表1提供，水泵启动的4个时刻不输入），这样第1未供水时段各时刻的水位可由MATLAB程序完成（程序代码见附件程序6）：

图8

（2）确定流量—时间函数

对于第1，2未供水时段的流量可直接对水位函数求导（程序代码见附件程序7）：

图9

wgsysl1=145.6657 wgsysl2=260.6561

下图为5次多项式拟合，相比较显然较三次拟合的效果好。

（程序代码见附件程序8）：

图10

而第1供水时段的流量则用前后时期的流量进行拟合得到。

为使流量函数在t=9和t=11连续，我们只取4个点，用三次多项式拟合得到第1供水时段的时间—流量图，可以看到与总时间—流速函数示意图（图7）中的相应部分比较吻合。

图11（程序代码见附件程序9）

gsysl1=49.8215。

在第2供水时段之前取t=19.96,20.84两点的流量，用第3未供水时段的3个记录做差分得到两个流量数据21.62，18.48，然后用这4个数据做三次多项式拟合得到第2供水时段与第3未供水时段的时间—流量图（图7），可以看到与总时间—流速函数示意图中的相应部分也比较吻合。

图12（程序代码见附件程序10）

gsysl2=73.9635

（3）一天总用水量的估计

分别对供水的两个时段和不供水的两个时段积分（流量对时间）并求和得到一天的总用水量约为530.1068（此数据是总用水高度，单位为cm）。

表6列出各段用水量，与插值法算得的表4相比，二者较为吻合。

表6

时段

全天用水

用水高度（cm）

145.6657

260.6561

49.8215

73.9635

530.1068

（4）结果分析

由表3可以看出，使用三次样条插值法得到的结果（145.6885，258.6547）与表1中记录的下降高度146cm，260cm相差不大，说明插值结果与原始数据比较吻合。

由表6可以得全天的用水量约为526.7148*233.3475*10=1229075.82升

（5）流量及总用水量的检验

计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。

例如，算得第1未供水段的用水量高度是145.67cm，而实际记录的水位下降高度为968－822=146cm，两者是吻合的；

算得第2未供水段的用水量高度是260.66cm，而实际记录的水位下降高度为1082－822=260cm，两者也是吻合的。

从算法设计和分析可知，计算结果与各时段所用的拟合多项式的次数有关。

表7给出的是对第1，2未供水时段用五多项式拟合后得到的用水量结果。

表7

146.5150

257.7605

46.1317

76.3076

526.7148

四、模型检验

1、分三段（水泵未工作）的实际用水量与由模型推算的用水量的差异很小，见表5，最

大的不超过4%，三段总计不超过0.5%。

其中由模型推算的用水量为水流量函数f（t）的积分；

实际用水量为水塔内水的体积之差。

表5分三段的实际用水量与模型计算出的24小时用水量比较（单位：

立方米）

实际用水量模型用水量绝对误差相对误差

第一段[0，8.967]345.3682338.05117.31712.1%

第二段[10.954，20.839]616.3161620.83414.51800.7%

第三段[22.958，25.908]151.7308156.60004.86923.2%

三段总计1113.41511115.48522.07010.2%

2、两次充水期间，水泵注入量的差异大约为3个立方米，不到0.5%。

水泵充水量=充水后的水量+充水期间的流出量—充水前的水量。

第一次水泵充水量=2564+117.4463-1948=733.4463立方米；

第二次水泵充水量=2564+120.7052-1948=736.7052立方米。

由于水泵功率=充水量/充水时间，由上述计算知，水泵的功率大约为367立方米/小时，而且两次冲水期间计算出的功率大致相等，这与实际问题是一致的。

3、用水高峰的比较

实际与模型之间相差无几。

实际用水高峰可近似地用差商最大值点表示为t=11，即上午11点钟左右；

模型得到的用水高峰可由模型得到的水流量函数依次在单位区间上积分或者找出水流量函数的最大值点得出，它们都在11点左右。

五、总结

1、优点

（1）模型灵活性好、稳定性强，可用于那些拥有地方性的竖直圆柱型水塔的小城镇和乡镇。

模型中的输入数据可以是任何近似均匀的时间间隔时的水位，时间间隔大约2小时。

（2）模型中的数学概念简单，并且容易理解。

只用到数值计算知识。

（3）模型容易实现，且给出了一天里水流速度和用水量的精确估计。

2、缺点

（1）本模型受水塔的几何形状限制。

（2）光滑曲线的逼近方法不能模拟真实流动中流速的微小变化，实际流动中流速可能会有一种高程度的“噪音”，即激波出现。

3、改进与推广

（1）我们可以在模型中用一个参数来限定不同几何形状的水塔。

（2）可以通过对流速数据进行回归分析等一系列处理，以便得到一些随机变化的特征。

参考文献

[1]亨塞尔曼.美.精通Matlab.清华大学出版社.2006-05

[2]马修斯.美数值方法（MATLAB版）（第4版）.电子工业出版社.2005-12

[3]司守奎.数学建模算法与程序.海军航空工程学院

[4]陆元鸿.数理统计方法[M].上海:

华东理工大学出版社.2005

[5]杨桂元等.数学模型应用实例.合肥工业大学.2007-06

[6]陈桂明，戚红雨，潘伟编著，Matlab数理统计（6.X），北京：

科学出版社，2002

[7]谢云荪，张志让.数学实验.北京：

科学出版社.2000

[8]飞思科技产品研发中心编著.MATLAB6.5辅助优化计算与设计.北京：

电子工

业出版社，2003。

附件

程序1：

t=[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.29];

t0=0:

0.1:

8.97;

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,'

:

'

t0,fdxx,'

-.'

t0,scyt,'

b'

）

gtext（'

lglr'

fdxx'

scyt'

程序2：

t=[8.45,8.97,10.95,11.49,12.49];

v=[16.35,19.29,33.50,29.63,31.52];

t0=8.97:

10.95;

spl