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

《数值分析课程设计三次样条插值.doc》由会员分享，可在线阅读，更多相关《数值分析课程设计三次样条插值.doc（26页珍藏版）》请在三一办公上搜索。

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

2、的估计问题进行求解是重要的解题方法。整个模型求解程序均是在MATLAB2010中运行的。关键字：城镇化、水流量、曲线插值、曲线拟合、MATLAB2010目录摘 要21理论基础41.1问题重述42 算法分析52.1 功能分析52.2 算法分析62.2.1 拉格朗日插值法62.2.2 三次样条插值72.2.3 最小二乘法的计算方法113 程序设计123.1 选单和主窗口设计123.2 模块设计224 总结225 参考文献241 理论基础1.1问题重述许多供水单位由于没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位。试通过测得的某时刻水箱中的水位的数据，估计在任意时刻（包括水泵灌水期间）流出

3、水箱的流量。给出原始数据表（如表1），其中长度单位为E （1E=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)0317544636335033163110499533260663530545393631

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

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

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

7、作时不维修，也不中途停止工作。水泵冲水的水流量远大于水塔的水流量。假设四：表1中水位数据取得的时间准确在1秒之内。假设五：水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小。2.2 算法分析通过上面各项功能的分析、分类、综合，按照模块化程序设计的要求，得到模块结构（作为后面选单设计的主要依据2.2.1 拉格朗日插值法假设取区间上的n+1个结点，并且已知函数f(x)在此点的函数值xx0x1x2xnf(x)f(x0)f(x1)f(x2)f(xn)现在求一个次数不超过n的多项式使得满足条件这种差值方法称为n次多项式插值(或称代数插值)，利用拉格朗日插值方法可得，其中上述多项式称为n次

8、拉格朗日(lagrange)插值多项式，函数称为拉格朗日插值基函数。特别地，当n1，2时，n次拉格朗日(lagrange)插值多项式即为线性差值多项式和抛物插值多项式2.2.2 三次样条插值分段线性插值，具有良好的稳定性和收敛性，但光滑性较差。在数学上若函数（曲线）的k阶导数存在且连续，则称该曲线具有k阶光滑性。易见，分段线性插值不光滑，这影响了它在某些工程技术实际问题中的应用。例如：在船体、飞机等外形曲线的设计中，不仅要求曲线连续而且还要求曲线的曲率连续，这就要求插值函数具有连续的二阶导数。为解决这一类问题，就产生了三次样条插值。所谓样条（Spline），本来是指一种绘图工具，它是一种富有弹

9、性的细长木条，在飞机或轮船制造过程中，被用于描绘光滑的外形曲线。使用时，用压铁将其固定在一些给定的节点上，在其他地方任其自然弯曲，然后依样画下的光滑曲线，就称为样条曲线。它实际上是由分段三次曲线拼接而成，在连续点即节点上，不仅函数自身是连续的，而且它的一阶和二阶导数也是连续的。从数学上加以概括，可得到样条函数的定义如下：三次样条函数记作，满足：在每个小区间是三次多项式。 在每个内节点上具有二次连续导数。 由三次样条函数中的条件知，有个待定系数。由条件知，在个内节点上具有二阶连续导数，即满足条件： （1）共有个条件。由条件，知，共有个条件。因此，要确定一个三次样条，还需要外加个条件.最常用的三次

10、样条函数的边界条件有两类：第一类边界条件：第二类边界条件：特别地，称为自然边界条件。第三类边界条件： 称为周期边界条件。三次样条插值不仅光滑性好，而且稳定性和收敛性都有保证，具有良好的逼近性质。样条插值函数的建立。构造满足条件的三次样条插值函数的表达式可以有多种方法。下面我们利用的二阶导数值表达，由于 在区间 上是三次多项式，故在上是线性函数，可表示为 （2）其中对积分两次并利用及，可定出积分常数，于是得三次样条表达式 （3）上式中是未知的，为确定，对求导得 （4） 由此可得。类似地可求出在区间上的表达式，从而得利用 可得 （5） 其中 （6）对第一类边界条件，可导出两个方程 （7）如果令，则

11、式（5）及其（7）可写出矩阵 （8）通过求解上述三对角矩阵可求得。对于第二类边界条件，直接得端点方程 （9）如果令，则式（5）及式（9）也可以写成矩阵（8）的形式。对于第三类边界条件，可得 （10）其中， 则式（5）及式（10）可以写成矩阵形式 求解上述矩阵可得 。2.2.3 最小二乘法的计算方法设为个线性无关的函数，对给定的数据，求，使 最小。利用极值的必要条件。 得到关于的线性方程组则方程组可表示为，其中，由于线性无关，所以是列满秩，是可逆矩阵，方程组的解存在且唯一，并且。取，得多项式拟合。3 程序设计3.1 选单和主窗口设计流量估计方法首先依照表1所给数据，用MATLAB作出时间水位散点

12、图（图1）。下面来计算水箱流量与时间的关系。 根据图1，一种简单的处理方法为，将表1中的数据分为三段，然后对每一段的数据做如下处理： 设某段数据 ，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位左端点的水位）/区间长度）： 每段数据首尾点的流速用下面的公式计算： 用以上公式求得时间与流速之间的数据如表2。表2时间/h流速/cmh-1时间/h流速/ cmh-1029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.42515.3217.3823.645.

13、4413.0418.4923.426.4515.4519.5025.007.46513.9820.4023.868.4516.3520.8422.178.9719.2922.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表2作出时间流速散点图如图2。插值法由表2，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由

14、于最后一段水泵不工作时段数据太少，我们将它忽略，只对水泵工作时段2进行插值处理。我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。 调用程序3t = 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);l

15、glrjf = 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结果 运行结果：lglrjf = 145.6250 fdxxjf =147.1469 sancytjf =145.6885图中曲线lglr、f

16、dxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。 由表1知，第1未供水时段的总用水高度为146(=968822)，可见上述三种插值方法计算的结果与实际值（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;t0 = 8.97 : 0.1 : 10.95; lglr = l

17、glrcz(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运行结果：lglrjf =56.4426 fdxxjf =49.6051 sancytjf =53.5903第

18、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.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 = i

19、nterp1(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运行结果：lglrjf = 258.8664 fdxxjf = 258.9697 sancytjf = 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、 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)gtext(lglr)gtext(fdxx)gtext(scyt)图6运行结果：lglrjf =104.1526fdxx

21、jf =73.7929sancytjf =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.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.

22、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.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运行结果：fdxxjf

23、 = 534.4311 sancytjf = 541.8148表3 各时段及一天的总用水量（用水高度(cm)）第1未供水段第2未供水段第1供水段第2供水段第3未供水段全天（除3未供）拉格朗日插值法145.625258.866456.4426104.152640.7825565.0866分段线性插值法147.1465258.969749.605173.792940.7516529.5142三次样条插值法145.6885258.654753.590379.817240.7825537.7507表4 各时段及一天的总用水量（用水体积(m3)）第1未供水段第2未供水段第1供水段第2供水段第3未供水段全

24、天（除3未供）拉格朗日插值法339.8123604.0583131.7074243.037595.16491318.616分段线性插值法343.3627604.2993115.7523172.193995.09281235.608三次样条插值法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.5

25、49029.4173注：拉格朗日插值法分段线性插值法三次样条插值法3.2.2 实验模型检验1、水泵未工作时的检验首先计算水泵未工作的三段时间内的实际用水量第一时间段，即t0，8.97时的实际用水量为：233.3475(9.68-8.22)=340.69第二时间段，即t10.95，20.84时的实际用水量为：233.3475(10.82-8.22)= 606.70第三时间段，即t23.88，25.91时的实际用水量为：233.3475(10.59-10.18)=95.67其次用三次样条插值模型得到的函数分别在水泵未工作的三段时间内作数值积分，执行程序后，可得水泵未工作的三段时间内，用模型计算出的

26、用水量（见表4-2）。分三段的实际用水量与模型用水量的比较可见，用三次样条插值模型计算的用水量与实际用水量吻合的很好。2、水泵工作时的检验水泵充水量=充水后的水量+充水期间的流出量-充水前的水量。首先利用模型计算两次充水期间的流出量，执行后可得：第一次充水期间的用水量为117.95，第二次充水期间的用水量为179.06。进一步可计算第一次水泵充水量=233.347510.82+117.95-233.34758.22=724.6535第二次水泵充水量=233.347510.59+179.06-233.34758.22=732.0936易见，水泵两次注水量相差不大。4 总结4.1模型的优缺点优点：

27、(1)这模型很灵活，具有可移植性，输入的数据可为任何规则逼近时间区间的水位，即时间分布可以是随机的。(2)模型用的数学知识简单易懂，曲线插值和曲线拟合的内容现在已普遍为人所知且在实际中广泛使用。(3)模型使用的软件操作简单，计算速度快。(4)模型不仅提供了水流量及一天用水量的较基本准确的估计，还可以估计任何时刻的水流量，包括水泵工作时的数流量。缺点：(1)用光滑曲线拟合的方法无法模拟真实水流量曲线的微小变化，除非时间取得足够小。(2)最小二乘法拟合的误差不能估计，这样做的结果可能对这一组的特定数据拟合得好坏有一定的影响。 4.2 实验心得体会 通过该课程设计，增加了自身对处理大型问题的能力，能

28、够认真的分析、尝试、写代码、最后生成论文。虽然过程中出现了许多困难，但凭借着一颗执着的心最终顺利完成了本实验报告。5 参考文献1 数值分析，李庆扬，王能超，易大义，2001，清华大学出版社（第四版）。2 数值方法，关治，陆金甫，2006，清华大学出版社。3数值分析与实验学习指导，蔡大用，2001，清华大学出版社。4数值分析与实验，薛毅，2005，北京工业大学出版社。5杨桂元等.数学模型应用实例.合肥工业大学.2007-066叶其孝.数学建模教育与国际数学建模工科数学专辑M .合肥:工科数学杂志社 ,1995课程设计成绩评定表评定项目内 容满分评分总分学习态度学习认真，态度端正，遵守纪律。10答疑和设计情况认真查阅资料，勤学好问，提出的问题有一定的深度，分析解决问题的能力较强。40说明书质量设计方案正确、表达清楚；设计思路、实验（论证）方法科学合理；达到课程设计任务书规定的要求；图、表、文字表达准确规范，上交及时。40回答问题情况回答问题准确，基本概念清楚，有理有据，有一定深度。10总成绩采用五级分制：优、良、中、及格、不及格指导教师评语： 签名： 年 月 日