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

1、题 目:姓 名:院（系）:专业班级:学 号: 指导教师: 时 间:部州轻工池曇協数值分析课程设计报告水箱水流量问题赵威威数学与信息科学学院信科11-01班541110010159汪远征2013年12月30日至2014年1月4日（摘要是论文内容的简短陈述，一般不超过200字。）随着城镇化的发展，水资源的分配与管理成为城市发展的核心问 题，高效合理地分配水资源成为评价水厂的重要指标。对自来水厂在尽 可能短的时间内准确对水塔水流量的预测在高规格的用水管理机构中 越来越成为至关紧要的一个环节，对城镇化的发展也具冇极英重要的意 义。分析该问题运用曲线插值和曲线拟合两种方法建立数学模型，对水 流量的估计问

2、题进行求解是重要的解题方法。整个模型求解程序均是在matlab2010中运行的。关键字：城镇化、水流量、曲线插值、曲线拟合、matlab2010目录摘 要ii1理论基础11问题重述12算法分析22功能分析22.2算法分析32.2.1拉格朗fi插值法32. 2.2三次样条插值42.2.3最小二乘法的计算方法83程序设计93.1选单和主窗口设计93.2模块设计134 总结215 参考文献221理论基础1.1问题重述许多供水单位由于没有测量流入或流出水箱流量的设备，而只能测量水 箱中的水位。试通过测得的某吋刻水箱中的水位的数据，估计在任意吋刻(包 括水泵灌水期间)(流出水箱的流fi/(r)o给出原始

3、数据表(如表1),其中 长度单位为e (le=30.24cm),水箱为圆柱体，其直径为57 e。假设：(1) 影响水箱流量的唯一因索是该区公众对水的普通需要；(2) 水泵的灌水速度为常数；(3) 从水箱屮流出水的最人流速小于水泵的灌水速度；(4) 每天的川水量分布都是相似的；(5) 水箱的流水速度可用光滑曲线來近似；(6) 当水箱的水容最达到514.8 kg时，开始泵水；达到677.6 kg时，停 止泵水。表1时间(s)水位(10 -2 e)吋间(s)水位(10 -2 e)0317544636335033163110499533260663530545393631671061929945725

4、4308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254泵水35932泵水82649泵水39332泵水8596834753943535508995333974331834459327033401.2需求分析通过对水流速进行估计，越是正确的估计将节省各方面额外的开销, 在相同的资源配置下获取最大的利润。2算法分析2.1功能分析木文问题所指流量可视为单位吋间内流出水的体积，-天的吋间为 0到23.88小时止（即忽略第3未供水时段）。由于水

5、塔是正圆柱形, 横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时 间的变化算出。问题的难点在于如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得 到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量 越准确越好。这此流量大体上可由两种方法计算，一是直接对表1中的 水量用数值微分算出各吋段的流量，用它们拟合其它吋刻或连续吋间的 流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续时 间的流量。有了任何时刻的流量，就不难计算一天的总用水量。其中，水泵不 工作时段的用水量可以由测量记录直接得到，由表1屮下降水位乘以水 搭的截而积就是这

6、一时段的用水量。这个数值可以用来检验数据插值或 拟合的结果。符号说明t：测量的吋刻；h：水位的咼度；r：水塔屮水的体积；模型假设假设一：流量只取决于水位差，与水位本身无关，故由物理学中 torriceli定律：小孔流出的液体的流速正比于水面高度的平 方根。题目给出水塔的最低和最高水位分别是& 1648m （设出 口 的水位为零）因为 sqrt （10.7352/8.1648） « 1.1467 ,约为 1, 所以可忽略水位对流速的影响。假设二：将流量看作时间的连续光滑函数，为计算简单，不妨将流量定 义成单位吋间流出水的高度，即水位对吋间变化率的绝对值 （水位是下降的），水塔截

7、面积为 s = （57*0.3024）2 *刃4 = 233.3475（），得到结果后乘以s即可.假设三：水泵工作起止时间由水塔的水位决定。水泵工作时不维修，也 不中途停止工作。水泵冲水的水流量远大于水塔的水流量。假设四：表1屮水位数据取得的吋间准确在1秒之内。假设五：水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少 水流量的大小。2.2算法分析通过上面各项功能的分析、分类、综合，按照模块化程序设计的要 求，得到模块结构（作为后面选单设计的主要依据2. 2.1拉格朗日插值法假设取区间“上的n+1个结点并且已知函数f（x）在此点的函数值xxox1x-2 xnf（x）i f（xo） |f（x

8、j | f（x2）| |f（xj现在求一个次数不超过n的多项式pn （x） = a0 + axx + a2x + + ©/"使得满足条件代（兀）二/（兀），心0,1,这种差值方法称为n次多项式插值（或称代数插值）,利用拉格册日 插值方法可得代= d（qfg）+仁/（和+厶.“/（兀2）+仁（兀）/（£），" r y其中人=口丄虫 “0,1,'戶0耳一旺上述多项式称为n次拉格朗r （lagrange）插值多项式，函数 仁（兀），2 =0丄，兀称为拉格朗日插值基函数。特别地，当n=l, 2时，n次拉格朗h （lagrange）插值多项式即为 线性差值

9、多项式和抛物插值多项式2. 2. 2三次样条插值分段线性插值，具冇良好的稳定性和收敛性，但光滑性较差。在数 学上若函数（曲线）的k阶导数存在口连续，则称该曲线具有k阶光滑 性。易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题 中的应用。例如：在船体、飞机等外形曲线的设计中，不仅耍求曲线连 续而且还耍求曲线的曲率连续，这就要求插值函数具有连续的二阶导 数。为解决这一类问题，就产生了三次样条插值。所谓样条（spline）,本来是指一种绘图工具，它是一种富有弹性 的细长木条，在飞机或轮船制造过程屮，被用于描绘光滑的外形曲线。 使用时，用压铁将其固定在一些给定的节点上，在其他地方任其自然弯

10、曲，然后依样画下的光滑曲线，就称为样条曲线。它实际上是由分段三 次曲线拼接而成，在连续点即节点上，不仅函数门身是连续的，而且它 的一阶和二阶导数也是连续的。从数学上加以概扌舌，可得到样条函数的 定义如下：三次样条函数记作s(x) , a<x<b ,满足： 在每个小区间x,.,和(心0丄必-1)是三次多项式。 在每个内节点耳0 = 1,/-1)上具有二次连续导数。 soj = x，i = 0 丄/由三次样条函数屮的条件知，s(x)有加个待定系数。由条件知，5(%)在t个内节点上具有二阶连续导数，即满足条件：s(兀一 o) = s(xj+o)< sa-o) = sa+o) q =

11、 0,l,2,/-l)(1)s"ao)=s”y+o)共有35-1)个条件。由条件，知s(e = x，心0丄曲，共有+ 1个条件。因此，要确定一个三次样条，还需要外加4h-3(h-1)-(/7 + 1) = 2 个条件.最常用的三次样条函数scx)的边界条件有两类：第一类边界条件:sxq) = f； 9sxn) = f；第二类边界条件：sxxq) = f；ysxxn) = f；特别地，sn(xo) = o,svj = o,称为自然边界条件。第三类边界条件:s(x0) = s(£),s 仇)二 s xn), sff(x0) = sg 称为周期边界条件。三次样条插值不仅光滑性好，

12、而且稳定性和收敛性都有保证，貝有 良好的逼近性质。样条插值函数的建立。构造满足条件的三次样条插值函数s(q的表达式可以有多种方法。下面我们利用s(x)的二阶导数值s”(®) = m2 = 0丄加表达s(x),由于s在区间xrxj+l上是三次多项式,故s”(兀)在xpxj+上是线性函数，可表示为x xhjx-x.jhj(2)其中i* = xy+1 -xj对sx)积分两次并利用5(xy)=州及s(s)=儿+i,可定出积分常数，于是得三次样条表达式6花s(x) = mi(x-兀)6hj0=0 丄一1)上式中= 丄/)是未知的，为确定mj,(j = 0,l,加，对s(q求 导得=二空+ m旧

13、上n +畑二比一仏也勺(4)72/2.7+12/1.h.67jjj由此可得s'(厂+ 0) = _zm zm田+_儿。73767+16类似地可求出s(x)在区间is，叩上的表达式，从而得sg 0)=丛+丛側,+片_八767-13 j hm 利用sq+o)= sqo)可得= 6fxj_l,xj,xj+l其小d 二§门®'b+j_门®t'® 厂s + %fxj_vxrxj+hjfxj,xj+i =对第一类边界条件s go) = f，0,5,(xn) = f，可导出两个方程、2m°+m严/心卸-几)%*mn_2mn=-(f-f

14、xn.9xn)hn-l如果令g心訓“-几亡(几丁“|)则'2入、w。、仏、“2 a d “一12a,-1<儿2 >& ja丿式(5)及其(7)可写出矩阵(8)通过求解上述三对角矩阵可求得m(), a/】,m”。对于第二类边界条件，直接得端点方程(9)如果令入=他=0竝=2厂0，心=2厂”,则式(5)及式(9)也可以写 成矩阵(8)的形式。对于第三类边界条件，可得anml+pnmn_+2mtl=dn(10)其心*7七心兽"i,2-1 'n-则式(5)及式(10)可以写成矩阵形式2“2人 2入 如 2a, ym22八mn求解上述矩阵可得m|,m2，,m

15、t，m”。2. 2. 3最小二乘法的计算方法设 cl),© (兀), (兀)为加个线性无关的函数，对给定的数据a，x),(zt,2,m),求 p(x) =(x) + a2(p2 (x) + 4- am(pm (x),使q(q"2,)=工皆=/=1龙兀-") 最小。利用极值的必要条件 /=!学= 0(" 1,2,加)。得到关于4宀，4的线性方程组“加s叭（无）（x 一 工ak探（） = 0 /=!*=1s02a）（x 一工条久 a）= 01=1r=1工（无）（）1-工色5（兀）= 0.；=1k=则方程组可表示为g1 ga = gjy ,其屮a =（坷宀,，

16、so） ®（兀j（西）丫 = 0，2,，儿 丫，g =（勺）0（兀2）”（兀2） 4（暫）©（£）久（£）丿”枷由于%,（兀）线性无关，所以g是列满秩，gt是可逆矩阵,方程组的解存在且唯一，并ra = （gtgygty o取 =1,°2（兀）=兀，,0加=， 得 多项式 拟 合p(x) = a + a2x + + amxm。3程序设计3.1选单和主窗口设计流量估计方法首先依照表1所给数据，用matlab作出吋间一水位散点图（图1）。1110.5*9.5、*98.58051015202530下面来计算水箱流量与时间的关系。根据图1, 一种简单的处

17、理方法为，将表1小的数据分为三段，然 后对每一段的数据做如下处理：设某段数据（兀0，儿），（初）,，（£，儿）,相邻数据中点的平均流速用下面 的公式（流速二（右端点的水位一左端点的水位）/区间长度）:1（兀+1+兀）=x -冶】2和一无每段数据首尾点的流速用下而的公式计算："（兀0） = （3% - 4） + 儿）/（兀2 - 兀0）"（耳）=（一3儿 + 4儿“一儿一2）/（£ 一 £-2）用以上公式求得时间与流速z间的数据如表2。 表2时间/h流速/cm h-1时间/h流速/ cm h-1029. 8912. 4931. 520. 4621

18、.7413. 4229. 031.3818.4814. 4326. 362. 39516. 2215. 4426. 093.4116. 3016. 3724. 734. 42515. 3217. 3823. 645. 4413. 041& 4923. 426. 4515. 4519. 5025. 007. 46513. 9820. 4023. 86& 4516. 3520. 8422. 178. 9719. 2922.02水泵开动9.98水泵开动22. 96水泵开动10. 93水泵开动23.8827.0910. 9533. 5024. 4321.6211.4929. 6325.

19、 4518. 4825.9113. 3035100图225201551025303i1520由表2作出时间一流速散点图如图2o插值法由表2,对水泵不工作时段1, 2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。我们分别采取拉格朗h插值法，分段线性插值法及三次样条插值法；对 于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不 工作时段数据太少，我们将它忽略，只对水泵工作时段2进行插值处理。 我们总共需要对四段数据（第1, 2未供水时段，第1供水时段，混合 时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法 算出流量函数和用水量（用水高度）。调用程序3t

20、 = 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 :& 97;lglr = lglrcz (t, v, to);lglrjf = 0. 1 * trapz(lglr)fdxx = interpl (t, v, to);fdxxjf 二 0. 1 * trapz (fdxx)s

21、cyt = inte rpl (t, v, to, 'spline');sancytjf 二 0 1 * trapz (scyt)plot (t, v, ', to, lglr,to, fdxx，to, scyt, ' b')gtext c lglr') gtext (' fdxx') gtext (' scyt，)实现图3结果运行结果：lglrjf 二 145. 6250 fdxxjf =147. 1469sancytjf=145. 6885图中曲线lglr. fdxx和scyt分别表示用拉格朗口插值法，分段线 性插值法

22、及三次样条插值法得到的曲线。由表1知，第1未供水时段的总用水高度为146(=968-822),可见上 述三种插值方法计算的结果与实际值(146)相比都比较接近。考虑到 三次样条插值方法具有更加良好的性质，建议采取该方法。其他三段的处理方法与第1未供水时段的处理方法类似.3.2模块设计3. 2. 1四个主要阶段t-v代码第一供水段时间一流速示意图程序4t = 8. 45, 8. 97, 10. 95, 11.49,12. 49;v = 16. 35, 19. 29, 33. 50, 29. 63, 31. 52;to 二 8.97 : 0. 1 : 10. 95;lglr = lglrcz(t,

23、 v, to);lglrjf = 0. 1 * trapz(lglr)fdxx = in terpl (t, v, to):fdxxjf = 0. 1 * trapz (fdxx)scyt = in tcrpl (t, v, to, 'spline');sancyt jf = 0. 1 * trapz (scyt)plot (t, v, ', to, lglr, j :, to, fdxx,'-.', to, scyt, ' b)gtext( lgl?)gtext (' fdxx')gtext (' scyl')第

24、一供水段时间流速图fdxxjf =49.6051运行结果：lglrjf 二56. 4426sancytjf =53.5903第2未供水段时间一流速示意图程序5t二10. 95, 11. 49, 12. 49, 13. 42, 14. 43, 15. 44, 16. 37, 17. 38, 18. 49, 19. 50, 20. 40, 20. 849v=33. 50, 29. 63, 31. 52, 29. 03, 26. 36, 26. 09, 24. 73, 23. 64, 23. 42, 25. 00, 23. 86, 22. 179to 二 10. 95 : 0. 1 : 20.84;

25、lglr = lglrcz(t, v, t0);lglrjf = 0. 1 * trapz(lglr)fdxx = interpl (t, v, to);fdxxjf = 0. 1 * trapz(fdxx)scyt = interpl (t, v, to, 'spline');sancytjf 二 0. 1 * trapz(scyt)plot(t, v, ', to, lglr,to, fdxx，'一.'，to, scyt, ' b')gtext(' lglr')gtext (? fdxx')gtext (

26、9; scyt')图5第2未供水段时间一流速示意图运行结果：lglrjf = 258. 8664fdxxjf = 25& 9697sancytjf = 258.6547第2供水段吋间一流速示意图程序6t = 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;iglr = lglrcz(t, v, 10);iglrjf = 0. 1 * trapz(lglr)fdxx = interp

27、l (t, v, to);fdxxjf = 0. 1 * trapz(fdxx)scyt = in terpl (t, v, to, ' spl in e');sancytjf = 0. 1 * trapz(scyt)plot(t, v, ', to, lglr,to, fdxx，'一'，to, scyt, ' b)gtext (/ lgir*)gtext c fdxx，)gtext c scyt，)第2供水段时间一流速示意图图6iglr454035302520151020212223242526运行结果：lglrjf =104. 1526 fd

28、xxjf =73. 7929 sancytjf =79.8172下图是用分段线性及三次样条插值方法得到的整个过程的时间一流速函 数示意图。程序7t =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. 3& 18. 49, 19. 50, 20. 40, 20. 84, 23. 8& 24. 43, 25. 45, 25 .91;29. 89, 21. 74,48

29、, 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.62, 18. 4& 13. 30;t0 = 0 : 0. 1 : 23. 88;fdxx = interpl (t, v, to);fdxxjf = 0. 1 * trapz(fdxx)scyt = interpl (t, v, to, '

30、;spline');sancytjf 二 0. 1 * trapz(scyt)plot (t, v, ',10, fdxx, 10, scyt, ' b) gtext (? fdxx') gtext (' scyt，)运行结果：fdxxjf = 534. 4311sancyt jf = 541. 8148表3各时段及一天的总用水量(用水高度(cm)第1未 供水段第2未 供水段第1供 水段第2供 水段第3未 供水段全天（除3 未供）拉格朗 日插值 法145. 62525& 866456. 4426104. 152640. 7825565. 0866

31、分段线 性插值 法147. 146525& 969749. 605173. 792940. 7516529.5142三次样 条插値 法145. 6885258. 654753. 590379.817240. 7825537. 7507表4各吋段及一天的总用水量（用水体积（h?）第1未 供水段第2未 供水段第1供 水段第2供 水段第3未 供水段全天（除3 未供）拉格 朗f1 插值 法339.8123604.0583131.7074243. 037595. 16491318.616善法343. 3627604. 2993115. 7523172. 193995. 09281235.608證

32、龌法339. 9605603. 5643125. 0516186. 251495. 16491254.828表5是对一天中任取的4个时刻分別用3种方法得到的水塔水流量近似 值。时间（h）6. 8810. 8815. 8822. 8815. 982633. 742625. 566234.709914.827232.997625. 446525. 471515. 052733. 708925. 549029.4173注：拉格朗日插值法分段线性插值法三次样条插值法322实验模型检验1、水泵未工作时的检验首先计算水泵未工作的三段时间内的实际用水量第一吋间段，即t e 0, 8.97吋的实际用水量为：233. 3475 x(9. 68-8.22)=340. 69第二时间段，即t e10. 95, 20. 84时的实际用水量为：233. 3475 x(10. 82-8. 22)= 606. 70第三时间段,即te 23. 88, 25.91时的实际用水量为：233. 3475 x(10. 59-10. 18)=95. 67其次用三次样条插值模型得到的函数分别在水泵未工作的三段时间内 作数值积分，执行程序后，可得水泵未工作的三段时间内，用模型计算 出的用水量(见表4-2)。分三段的实际用水量与模型用水量