数值分析课程设计三次样条插值

郑州轻工业学院

《数 值 分 析》 课 程 设 计 报 告

题 目： 姓 名： 院（系）： 专业班级： 学 号： 指导教师：

水箱水流量问题

赵威威

数学与信息科学学院

信科11-01班 541110010159 汪远征

时 间： 2013年12月30日至2014年1月4日

摘 要

（摘要是论文内容的简短陈述，一般不超过200字。）

随着城镇化的发展，水资源的分配与管理成为城市发展的核心问题，高效合理地分配水资源成为评价水厂的重要指标。对自来水厂在尽可能短的时间内准确对水塔水流量的预测在高规格的用水管理机构中越来越成为至关紧要的一个环节，对城镇化的发展也具有极其重要的意义。分析该问题运用曲线插值和曲线拟合两种方法建立数学模型，对水流量的估计问题进行求解是重要的解题方法。

整个模型求解程序均是在MATLAB2010中运行的。

关键字：城镇化、水流量、曲线插值、曲线拟合、MATLAB2010

目录

摘 要 .......................................................................................................................... II 1 理论基础 .................................................................................................................. 1

1.1问题重述 ............................................................................................................ 1 2 算法分析 .................................................................................................................... 2

2.1 功能分析 ......................................................................................................... 2 2.2 算法分析 ......................................................................................................... 3

2.2.1 拉格朗日插值法 ................................................................................... 3 2.2.2 三次样条插值 ....................................................................................... 4 2.2.3 最小二乘法的计算方法 ....................................................................... 8

3 程序设计 .................................................................................................................... 9

3.1 选单和主窗口设计 .......................................................................................... 9 3.2 模块设计 ....................................................................................................... 13 4 总结 .......................................................................................................................... 21 5 参考文献 .................................................................................................................. 22

I

1 理论基础

1.1问题重述

许多供水单位由于没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位。试通过测得的某时刻水箱中的水位的数据，估计在任意时刻（包括水泵灌水期间）t流出水箱的流量f(t)。给出原始数据表（如表1），其中长度单位为E （1E=30.24cm），水箱为圆柱体，其直径为57 E。 假设：

(1) 影响水箱流量的唯一因素是该区公众对水的普通需要； (2) 水泵的灌水速度为常数；

(3) 从水箱中流出水的最大流速小于水泵的灌水速度； (4) 每天的用水量分布都是相似的； (5) 水箱的流水速度可用光滑曲线来近似；

(6) 当水箱的水容量达到514.8 Kg时，开始泵水；达到677.6 Kg时，停止泵水。

表1

时间 (s) 0 3316 6635 10619 13937 17921 21240 25223 28543 32284 水位(10 –2 E) 3175 3110 3054 2994 2947 2892 2850 2795 2752 2697 时间 (s) 44636 49953 53936 57254 60574 64554 68535 71854 75021 79254 水位(10 –2 E) 3350 3260 3167 3087 3012 2927 2842 2767 2697 泵水 1

35932 39332 39435 43318 泵水 泵水 3550 3445 82649 85968 89953 93270 泵水 3475 3397 3340 1.2 需求分析

通过对水流速进行估计，越是正确的估计将节省各方面额外的开销，在相同的资源配置下获取最大的利润。

2 算法分析

2.1 功能分析

本文问题所指流量可视为单位时间内流出水的体积，一天的时间为0到23.88小时止（即忽略第3未供水时段）。由于水塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段的流量。

水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这此流量大体上可由两种方法计算，一是直接对表1中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续时间的流量。

有了任何时刻的流量，就不难计算一天的总用水量。其中，水泵不工作时段的用水量可以由测量记录直接得到，由表1中下降水位乘以水搭的截面积就是这一时段的用水量。这个数值可以用来检验数据插值或拟合的结果。

符号说明

t：测量的时刻；

2

h：水位的高度； v：水塔中水的体积；

模型假设

假设一：流量只取决于水位差，与水位本身无关，故由物理学中

Torriceli定律：小孔流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是8.1648m（设出

口的水位为零）因为sqrt(10.73528.1648)1.1467，约为1，

所以可忽略水位对流速的影响。

假设二：将流量看作时间的连续光滑函数，为计算简单，不妨将流量定

义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），水塔截面积为

S(57*0.3024)2*4233.3475m2,得到结果后乘以

s即可.

假设三：水泵工作起止时间由水塔的水位决定。水泵工作时不维修，也不中途停止工作。水泵冲水的水流量远大于水塔的水流量。

假设四：表1中水位数据取得的时间准确在1秒之内。

假设五：水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少

水流量的大小。

2.2 算法分析

通过上面各项功能的分析、分类、综合，按照模块化程序设计的要求，得到模块结构（作为后面选单设计的主要依据

2.2.1 拉格朗日插值法

假设取区间a,b上的n+1个结点x0,x1,在此点的函数值 x x0 x1 3

,xn，并且已知函数f(x)

… xn x2 f(x) f(x0) f(x1) f(x2) … f(xn) 现在求一个次数不超过n的多项式

2Pn(x)a0a1xa2xanxn

使得满足条件

Pn(xi)f(xi),i0,1,,n

这种差值方法称为n次多项式插值(或称代数插值)，利用拉格朗日插值方法可得

Pn(x)l0,n(x)f(x0)l1,n(x)f(x1)l2,n(x)f(x2)其中lkn(x)j0jknln,n(x)f(xn)，

xxjxkxjk0,1,,n

上述多项式称为n次拉格朗日(lagrange)插值多项式，函数

lk,n(x),k0,1,,n称为拉格朗日插值基函数。

特别地，当n＝1，2时，n次拉格朗日(lagrange)插值多项式即为线性差值多项式和抛物插值多项式

2.2.2 三次样条插值

分段线性插值，具有良好的稳定性和收敛性，但光滑性较差。在数学上若函数（曲线）的k阶导数存在且连续，则称该曲线具有k阶光滑性。易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。例如：在船体、飞机等外形曲线的设计中，不仅要求曲线连续而且还要求曲线的曲率连续，这就要求插值函数具有连续的二阶导数。为解决这一类问题，就产生了三次样条插值。

所谓样条（Spline），本来是指一种绘图工具，它是一种富有弹性的细长木条，在飞机或轮船制造过程中，被用于描绘光滑的外形曲线。使用时，用压铁将其固定在一些给定的节点上，在其他地方任其自然弯曲，然后依样画下的光滑曲线，就称为样条曲线。它实际上是由分段三次曲线拼接而成，在连续点即节点上，不仅函数自身是连续的，而且它

4

的一阶和二阶导数也是连续的。从数学上加以概括，可得到样条函数的定义如下：

三次样条函数记作S(x)，axb，满足： ①在每个小区间[xi,xi1](i0,1, ②在每个内节点xi(i1, ③S(xi)yi,i0,1,,n1)是三次多项式。

,n1)上具有二次连续导数。

,n

由三次样条函数中的条件①知，S(x)有4n个待定系数。由条件②知，

S(x)在n1个内节点上具有二阶连续导数，即满足条件：S(xi0)S(xi0)S'(xi0)S'(xi0) (i0,1,2,S(x0)S\"(x0)ii,n1) （1）

共有3(n1)个条件。由条件③，知S(xi)yi,i0,1,,n，共有n1个

条件。因此，要确定一个三次样条，还需要外加4n3(n1)(n1)2个条件.

最常用的三次样条函数S(x)的边界条件有两类：

第一类边界条件：S'(x0)f0,S'(xn)fn 第二类边界条件：S''(x0)f0,S''(xn)fn 特别地，S''(x0)0,S''(xn)0，称为自然边界条件。

第三类边界条件： S(x0)S(xn),S'(x0)S'(xn),S(x0)S(xn) 称为周期边界条件。

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

5

构造满足条件的三次样条插值函数S(x)的表达式可以有多种方法。下面我们利用S(x)的二阶导数值S\"(xj)Mj(j0,1,n,表达)由于 S(x)在区间[xj,xj1] 上是三次多项式，故S\"(x)在[xj,xj1]S(x)，

上是线性函数，可表示为 S\"(x)Mjxj1xhjMj1xxjhj （2）

其中hjxj1xj对S\"(x)积分两次并利用S(xj)yj及S(xj1)yj1，可定出积分常数，于是得三次样条表达式

S(x)Mj(xj1x)36hjMj1(xxj)36hjMjh2xj1xMj1h2xxjjj(yj)(yj1) (j0,1,6hj6hj,n1) （3）

上式中Mj(j0,1,导得

S'(x)Mj(xj1x)22hjMj1(xxj)22hjyj1yjhjMj1Mj6hj （4）

,n)是未知的，为确定Mj,(j0,1,,n)，对S(x)求

由此可得S'(xj0)hj3Mjhj6Mj1yj1yj6。

类似地可求出S(x)在区间[xj1,xj]上的表达式，从而得

S'(xj0)hj16Mj1hj13Mjyjyj1hj1

利用 S'(xj0)S'(xj0)可得

6

jMj12MjjMj1dj(j1,2,n1) （5）

其中 jhj1hj1hj,jhjhj1hj,j0,1,,n

dj6f[xj,xj1]f[xj1,xj]hj1hjf(xj1)f(xj)hj6f[xj1,xj,xj1] （6）

f[xj,xj1]

f[xj1,xj,xj1]f[xj1,xj1]f[xj1,xj]hj

对第一类边界条件S'(x0)f'0,S'(xn)f'n，可导出两个方程

  （7）

6Mn12Mn(f'nf[xn1,xn])hn12M0M1如果令01,d066(f[x0,x1]f'0),n1，dn(f'nf[xn1,xn])则h0hn16(f[x0,x1]f'0)h0式（5）及其（7）可写出矩阵

02M0d0Md21111 （8） 2Mdn1n1n1n1n2Mndn通过求解上述三对角矩阵可求得M0,M1,对于第二类边界条件，直接得端点方程

,Mn。

M0f\"0,Mnf\"n （9）

7

如果令000,d02f\"0,dn2f\"n，则式（5）及式（9）也可以写成矩阵（8）的形式。

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

M0Mn nM1nMn12Mndn （10） 其中nh0hn1f[x0,x1]f[xn1,xn]，dn6 ,n1nhn1h0hn1h0h0hn1 则式（5）及式（10）可以写成矩阵形式

22n122n12nMd22 n1Mn1dn12Mndn1M1d1求解上述矩阵可得 M1,M2,,Mn1,Mn。

2.2.3 最小二乘法的计算方法

设1(x),2(x),,m(x)为m个线性无关的函数，对给定的数据

(xi,yi),(i1,2,Q(a1,a2,,n)，求P(x)a11x(a)2x2()ammx，(使)

n2in,am)[yiP(xi)]2 最小。利用极值的必要条件

i1i1Q0(k1,2,ak,m)。 得到关于a1,a2,,am的线性方程组

8

mn1(xi)(yiakk(xi))0k1i1mn2(xi)(yiakk(xi))0 i1k1mnm(xi)(yiakk(xi))0k1i1则方程组可表示为GTGaGTY，其中a(a1,a2,1(x1)2(x1)1(x2)2(x2),yn)T，G(x)(x)2n1n,am)T，

Y(y1,y2,m(x1)m(x2)m(xn)nm

由于{1(x),,m(x)}线性无关，所以G是列满秩，GTG是可逆矩阵，

方程组的解存在且唯一，并且a(GTG)1GTY。 取

1(x)1,2(x)x,,m(x)xmamxm。

，得多项式拟合

P(x)a1a2x

3 程序设计

3.1 选单和主窗口设计

流量估计方法

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

9

图11110.5109.598.58051015202530

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

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

{(x , y 设某段数据 y 0 ),( x 1 ,( x n , y n )} ，相邻数据中点的平均流速用下面0,1),的公式（流速=（右端点的水位－左端点的水位）/区间长度）：

v(xi1xiyyi1 )i2xi1xi每段数据首尾点的流速用下面的公式计算：

v(x0)(3y04y1y2)(x2x0)

v(xn)(3yn4yn1yn2)(xnxn2)

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

表2

时间/h 0 0.46 1.38 2.395 流速/cm·h-1 29.89 21.74 18.48 16.22 10

时间/h 12.49 13.42 14.43 15.44 流速/ cm·h-1 31.52 29.03 26.36 26.09 3.41 4.425 5.44 6.45 7.465 8.45 8.97 9.98 10.93 10.95 11.49

16.30 15.32 13.04 15.45 13.98 16.35 19.29 水泵开动 水泵开动 33.50 29.63 16.37 17.38 18.49 19.50 20.40 20.84 22.02 22.96 23.88 24.43 25.45 25.91 24.73 23.64 23.42 25.00 23.86 22.17 水泵开动 水泵开动 27.09 21.62 18.48 13.30 由表2作出时间—流速散点图如图2。

图2353025201510051015202530

插值法

由表2，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的

11

流速，从而可以知道任意时刻的流量。

我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它忽略，只对水泵工作时段2进行插值处理。 我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。 调用程序3

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') gtext('fdxx') gtext('scyt')

实现图3结果

12

第一未供水时间-流速图30282624222018161412fdxxscytlglr0123456789

运行结果：lglrjf = 145.6250 fdxxjf =147.1469 sancytjf =145.6885

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

由表1知，第1未供水时段的总用水高度为146(=968－822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。

其他三段的处理方法与第1未供水时段的处理方法类似.

3.2 模块设计

3.2.1 四个主要阶段t-v代码

第一供水段时间—流速示意图 程序4

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

13

v = [16.35,19.29,33.50,29.63,31.52]; t0 = 8.97 : 0.1 : 10.95;

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') gtext('fdxx') gtext('scyt')

图4

第一供水段时间流速图363432302826242220181688.599.51010.51111.51212.5lglrscytfdxx

运行结果：lglrjf =56.4426 fdxxjf =49.6051 sancytjf =53.5903

14

第2未供水段时间—流速示意图 程序5

t=[10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84];

v=[33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17];

t0 = 10.95 : 0.1 : 20.84;

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') gtext('fdxx') gtext('scyt')

图5

15

第2未供水段时间—流速示意图3432302826lglrscytfdxx24222010121416182022

运行结果：lglrjf = 258.8664 fdxxjf = 258.9697 sancytjf = 258.6547

第2供水段时间—流速示意图 程序6

t = [20.40,20.84,23.88,24.43,25.45,25.91]; v = [23.86,22.17,27.09,21.62,18.48,13.3]; t0 = 20.84 : 0.1 : 23.88;

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

16

gtext('lglr') gtext('fdxx') gtext('scyt')

图6

第2供水段时间—流速示意图45lglr403530scytfdxx2520151020212223242526

运行结果：lglrjf =104.1526 fdxxjf =73.7929 sancytjf =79.8172

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

t =

[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97,10.95,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,20.84,23.88,24.43,25.45,25.91]; v =

[29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29,33.50,29.63,31.52,29.03,26.36,26.09,24.73,23.64,23.42,25.00,23.86,22.17,27.09,21.

17

62,18.48,13.30]; t0 = 0 : 0.1 : 23.88;

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, fdxx ,'-.', t0, scyt, 'b')

gtext('fdxx') gtext('scyt')

图7

总时间—流速函数示意图3530scyt25fdxx201510051015202530

运行结果：fdxxjf = 534.4311 sancytjf = 541.8148 表3 各时段及一天的总用水量（用水高度(cm)）

18

第1未第2未第1供供水段 供水段 水段 145.625 147.1465 145.6885 258.8664 258.9697 258.6547 56.4426 49.6051 53.5903 第2供水段 104.1526 73.7929 79.8172 全天第3未（除3供水段 未供） 40.7825 40.7516 40.7825 565.0866 529.5142 537.7507 拉格朗日插值法 分段线性插值法 三次样条插值法

表4 各时段及一天的总用水量（用水体积(m3)）

第1未第2未第1供供水段 供水段 水段 339.8123 604.0583 131.7074 第2供水段 243.0375 全天第3未（除3供水段 未供） 95.1649 1318.616 拉格朗日插值法 分段线性插值法 三次样条插值法

343.3627 604.2993 115.7523 172.1939 95.0928 1235.608 339.9605 603.5643 125.0516 186.2514 95.1649 1254.828 表5是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值。 时间(h) ① ② 6.88 15.9826 14.8272 10.88 33.7426 32.9976 19

15.88 25.5662 25.4465 22.88 34.7099 25.4715 ③ 15.0527 33.7089 25.5490 29.4173 注：①拉格朗日插值法②分段线性插值法③三次样条插值法

3.2.2 实验模型检验

1、水泵未工作时的检验

首先计算水泵未工作的三段时间内的实际用水量

第一时间段，即t∈[0，8.97]时的实际用水量为：233.3475×(9.68-8.22)=340.69

第二时间段，即t∈[10.95，20.84]时的实际用水量为：233.3475×(10.82-8.22)= 606.70

第三时间段，即t∈[23.88，25.91]时的实际用水量为：233.3475×(10.59-10.18)=95.67

其次用三次样条插值模型得到的函数分别在水泵未工作的三段时间内作数值积分，执行程序后，可得水泵未工作的三段时间内，用模型计算出的用水量（见表4-2）。分三段的实际用水量与模型用水量的比较可见，用三次样条插值模型计算的用水量与实际用水量吻合的很好。

2、水泵工作时的检验

水泵充水量=充水后的水量+充水期间的流出量-充水前的水量。首先利用模型计算两次充水期间的流出量，执行后可得：第一次充水期间的用水量为117.95，第二次充水期间的用水量为179.06。 进一步可计算

第一次水泵充水量=233.3475×10.82+117.95-233.3475×8.22=724.6535

第二次水泵充水量=233.3475×10.59+179.06-233.3475×8.22=732.0936

易见，水泵两次注水量相差不大。

20

4 总结

4.1模型的优缺点

优点：

(1)这模型很灵活，具有可移植性，输入的数据可为任何规则逼近时间区间的水位，即时间分布可以是随机的。

(2)模型用的数学知识简单易懂，曲线插值和曲线拟合的内容现在已普遍为人所知且在实际中广泛使用。

(3)模型使用的软件操作简单，计算速度快。

(4)模型不仅提供了水流量及一天用水量的较基本准确的估计，还可以估计任何时刻的水流量，包括水泵工作时的数流量。

缺点：

(1)用光滑曲线拟合的方法无法模拟真实水流量曲线的微小变化，除非时间取得足够小。

(2)最小二乘法拟合的误差不能估计，这样做的结果可能对这一组的特定数据拟合得好坏有一定的影响。

4.2 实验心得体会

通过该课程设计，增加了自身对处理大型问题的能力，能够认真的分析、尝试、写代码、最后生成论文。虽然过程中出现了许多困难，但凭借着一颗执着的心最终顺利完成了本实验报告。

21

5 参考文献

[1] 数值分析，李庆扬，王能超，易大义，2001，清华大学出版社（第四版）。 [2] 数值方法，关治，陆金甫，2006，清华大学出版社。 [3]数值分析与实验学习指导，蔡大用，2001，清华大学出版社。 [4]数值分析与实验，薛毅，2005，北京工业大学出版社。 [5]杨桂元等.数学模型应用实例.合肥工业大学.2007-06

[6]叶其孝.数学建模教育与国际数学建模《工科数学》专辑[M] .合肥:《工科数学》杂志社 ,1995

22

课程设计成绩评定表

评定项目 学习态度 答疑和设计情况 内 容 学习认真，态度端正，遵守纪律。 认真查阅资料，勤学好问，提出的问题有一定的深度，分析解决问题的能力较强。 设计方案正确、表达清楚；设计思路、说明书质量 实验（论证）方法科学合理；达到课程设计任务书规定的要求；图、表、文字表达准确规范，上交及时。 回答问题情况 总成绩 回答问题准确，基本概念清楚，有理有据，有一定深度。 采用五级分制：优、良、中、及格、不及格 10 40 40 满分 10 评分 总分 指导教师评语： 签名： 年 月 日

23