數(shù)值分析課程設(shè)計(jì)三次樣條插值

1、鄭州輕工業(yè)學(xué)院數(shù) 值 分 析課 程 設(shè) 計(jì) 報(bào) 告題 目：水箱水流量問(wèn)題姓 名：趙威威院（系）：數(shù)學(xué)與信息科學(xué)學(xué)院專業(yè)班級(jí)：信科11-01班學(xué) 號(hào)：541110010159指導(dǎo)教師：汪遠(yuǎn)征時(shí) 間：2013年12月30日至2014年1月4日摘 要（摘要是論文內(nèi)容的簡(jiǎn)短陳述，一般不超過(guò)200字。）隨著城鎮(zhèn)化的發(fā)展，水資源的分配與管理成為城市發(fā)展的核心問(wèn)題，高效合理地分配水資源成為評(píng)價(jià)水廠的重要指標(biāo)。對(duì)自來(lái)水廠在盡可能短的時(shí)間內(nèi)準(zhǔn)確對(duì)水塔水流量的預(yù)測(cè)在高規(guī)格的用水管理機(jī)構(gòu)中越來(lái)越成為至關(guān)緊要的一個(gè)環(huán)節(jié)，對(duì)城鎮(zhèn)化的發(fā)展也具有極其重要的意義。分析該問(wèn)題運(yùn)用曲線插值和曲線擬合兩種方法建立數(shù)學(xué)模型，對(duì)水流量

2、的估計(jì)問(wèn)題進(jìn)行求解是重要的解題方法。整個(gè)模型求解程序均是在matlab2010中運(yùn)行的。關(guān)鍵字：城鎮(zhèn)化、水流量、曲線插值、曲線擬合、matlab2010目錄摘 要21理論基礎(chǔ)41.1問(wèn)題重述42 算法分析52.1 功能分析52.2 算法分析62.2.1 拉格朗日插值法62.2.2 三次樣條插值72.2.3 最小二乘法的計(jì)算方法113 程序設(shè)計(jì)123.1 選單和主窗口設(shè)計(jì)123.2 模塊設(shè)計(jì)224 總結(jié)225 參考文獻(xiàn)24i1 理論基礎(chǔ)1.1問(wèn)題重述許多供水單位由于沒(méi)有測(cè)量流入或流出水箱流量的設(shè)備，而只能測(cè)量水箱中的水位。試通過(guò)測(cè)得的某時(shí)刻水箱中的水位的數(shù)據(jù)，估計(jì)在任意時(shí)刻（包括水泵灌水期間）流

3、出水箱的流量。給出原始數(shù)據(jù)表（如表1），其中長(zhǎng)度單位為e （1e=30.24cm），水箱為圓柱體，其直徑為57 e。 假設(shè)： (1) 影響水箱流量的唯一因素是該區(qū)公眾對(duì)水的普通需要； (2) 水泵的灌水速度為常數(shù)； (3) 從水箱中流出水的最大流速小于水泵的灌水速度； (4) 每天的用水量分布都是相似的； (5) 水箱的流水速度可用光滑曲線來(lái)近似； (6) 當(dāng)水箱的水容量達(dá)到514.8 kg時(shí)，開始泵水；達(dá)到677.6 kg時(shí)，停止泵水。表1時(shí)間 (s)水位(10 2 e)時(shí)間 (s)水位(10 2 e)031754463633503316311049953326066353054539363

4、16710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254泵水35932泵水82649泵水39332泵水8596834753943535508995333974331834459327033401.2 需求分析通過(guò)對(duì)水流速進(jìn)行估計(jì)，越是正確的估計(jì)將節(jié)省各方面額外的開銷，在相同的資源配置下獲取最大的利潤(rùn)。2 算法分析2.1 功能分析本文問(wèn)題所指流量可視為單位時(shí)間內(nèi)流出水的體積，一天的時(shí)間為0到23.88小時(shí)

5、止（即忽略第3未供水時(shí)段）。由于水塔是正圓柱形，橫截面積是常數(shù)，所以在水泵不工作時(shí)段，流量很容易根據(jù)水位相對(duì)時(shí)間的變化算出。問(wèn)題的難點(diǎn)在于如何估計(jì)水泵供水時(shí)段的流量。水泵供水時(shí)段的流量只能靠供水時(shí)段前后的流量經(jīng)插值或擬合得到。作為用于插值或擬合的原始數(shù)據(jù)，我們希望水泵不工作時(shí)段的流量越準(zhǔn)確越好。這此流量大體上可由兩種方法計(jì)算，一是直接對(duì)表1中的水量用數(shù)值微分算出各時(shí)段的流量，用它們擬合其它時(shí)刻或連續(xù)時(shí)間的流量；二是先用表中數(shù)據(jù)擬合水位一時(shí)間函數(shù)，求導(dǎo)數(shù)即可得到連續(xù)時(shí)間的流量。有了任何時(shí)刻的流量，就不難計(jì)算一天的總用水量。其中，水泵不工作時(shí)段的用水量可以由測(cè)量記錄直接得到，由表1中下降水位乘以水

6、搭的截面積就是這一時(shí)段的用水量。這個(gè)數(shù)值可以用來(lái)檢驗(yàn)數(shù)據(jù)插值或擬合的結(jié)果。 符號(hào)說(shuō)明t：測(cè)量的時(shí)刻；h：水位的高度；v：水塔中水的體積；模型假設(shè)假設(shè)一：流量只取決于水位差，與水位本身無(wú)關(guān)，故由物理學(xué)中torriceli定律：小孔流出的液體的流速正比于水面高度的平方根。題目給出水塔的最低和最高水位分別是8.1648m（設(shè)出口的水位為零）因?yàn)閟qrt，約為1，所以可忽略水位對(duì)流速的影響。假設(shè)二：將流量看作時(shí)間的連續(xù)光滑函數(shù)，為計(jì)算簡(jiǎn)單，不妨將流量定義成單位時(shí)間流出水的高度，即水位對(duì)時(shí)間變化率的絕對(duì)值（水位是下降的），水塔截面積為,得到結(jié)果后乘以s即可.假設(shè)三：水泵工作起止時(shí)間由水塔的水位決定。水泵

7、工作時(shí)不維修，也不中途停止工作。水泵沖水的水流量遠(yuǎn)大于水塔的水流量。假設(shè)四：表1中水位數(shù)據(jù)取得的時(shí)間準(zhǔn)確在1秒之內(nèi)。假設(shè)五：水塔的水流量與水泵狀態(tài)獨(dú)立，并不因水泵工作而增加或減少水流量的大小。2.2 算法分析通過(guò)上面各項(xiàng)功能的分析、分類、綜合，按照模塊化程序設(shè)計(jì)的要求，得到模塊結(jié)構(gòu)（作為后面選單設(shè)計(jì)的主要依據(jù)2.2.1 拉格朗日插值法假設(shè)取區(qū)間上的n+1個(gè)結(jié)點(diǎn)，并且已知函數(shù)f(x)在此點(diǎn)的函數(shù)值xx0x1x2xnf(x)f(x0)f(x1)f(x2)f(xn)現(xiàn)在求一個(gè)次數(shù)不超過(guò)n的多項(xiàng)式使得滿足條件這種差值方法稱為n次多項(xiàng)式插值(或稱代數(shù)插值)，利用拉格朗日插值方法可得，其中上述多項(xiàng)式稱為n

8、次拉格朗日(lagrange)插值多項(xiàng)式，函數(shù)稱為拉格朗日插值基函數(shù)。特別地，當(dāng)n1，2時(shí)，n次拉格朗日(lagrange)插值多項(xiàng)式即為線性差值多項(xiàng)式和拋物插值多項(xiàng)式2.2.2 三次樣條插值分段線性插值，具有良好的穩(wěn)定性和收斂性，但光滑性較差。在數(shù)學(xué)上若函數(shù)（曲線）的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。易見，分段線性插值不光滑，這影響了它在某些工程技術(shù)實(shí)際問(wèn)題中的應(yīng)用。例如：在船體、飛機(jī)等外形曲線的設(shè)計(jì)中，不僅要求曲線連續(xù)而且還要求曲線的曲率連續(xù)，這就要求插值函數(shù)具有連續(xù)的二階導(dǎo)數(shù)。為解決這一類問(wèn)題，就產(chǎn)生了三次樣條插值。所謂樣條（spline），本來(lái)是指一種繪圖工具，它是一種富有

9、彈性的細(xì)長(zhǎng)木條，在飛機(jī)或輪船制造過(guò)程中，被用于描繪光滑的外形曲線。使用時(shí)，用壓鐵將其固定在一些給定的節(jié)點(diǎn)上，在其他地方任其自然彎曲，然后依樣畫下的光滑曲線，就稱為樣條曲線。它實(shí)際上是由分段三次曲線拼接而成，在連續(xù)點(diǎn)即節(jié)點(diǎn)上，不僅函數(shù)自身是連續(xù)的，而且它的一階和二階導(dǎo)數(shù)也是連續(xù)的。從數(shù)學(xué)上加以概括，可得到樣條函數(shù)的定義如下：三次樣條函數(shù)記作，滿足：在每個(gè)小區(qū)間是三次多項(xiàng)式。 在每個(gè)內(nèi)節(jié)點(diǎn)上具有二次連續(xù)導(dǎo)數(shù)。 由三次樣條函數(shù)中的條件知，有個(gè)待定系數(shù)。由條件知，在個(gè)內(nèi)節(jié)點(diǎn)上具有二階連續(xù)導(dǎo)數(shù)，即滿足條件： （1）共有個(gè)條件。由條件，知，共有個(gè)條件。因此，要確定一個(gè)三次樣條，還需要外加個(gè)條件.最常用的三

10、次樣條函數(shù)的邊界條件有兩類：第一類邊界條件：第二類邊界條件：特別地，稱為自然邊界條件。第三類邊界條件： 稱為周期邊界條件。三次樣條插值不僅光滑性好，而且穩(wěn)定性和收斂性都有保證，具有良好的逼近性質(zhì)。樣條插值函數(shù)的建立。構(gòu)造滿足條件的三次樣條插值函數(shù)的表達(dá)式可以有多種方法。下面我們利用的二階導(dǎo)數(shù)值表達(dá)，由于 在區(qū)間 上是三次多項(xiàng)式，故在上是線性函數(shù)，可表示為 （2）其中對(duì)積分兩次并利用及，可定出積分常數(shù)，于是得三次樣條表達(dá)式 （3）上式中是未知的，為確定，對(duì)求導(dǎo)得 （4） 由此可得。類似地可求出在區(qū)間上的表達(dá)式，從而得利用 可得 （5） 其中 （6）對(duì)第一類邊界條件，可導(dǎo)出兩個(gè)方程 （7）如果令，

11、則式（5）及其（7）可寫出矩陣 （8）通過(guò)求解上述三對(duì)角矩陣可求得。對(duì)于第二類邊界條件，直接得端點(diǎn)方程 （9）如果令，則式（5）及式（9）也可以寫成矩陣（8）的形式。對(duì)于第三類邊界條件，可得 （10）其中， 則式（5）及式（10）可以寫成矩陣形式 求解上述矩陣可得 。2.2.3 最小二乘法的計(jì)算方法設(shè)為個(gè)線性無(wú)關(guān)的函數(shù)，對(duì)給定的數(shù)據(jù)，求，使 最小。利用極值的必要條件。 得到關(guān)于的線性方程組則方程組可表示為，其中，由于線性無(wú)關(guān)，所以是列滿秩，是可逆矩陣，方程組的解存在且唯一，并且。取，得多項(xiàng)式擬合。3 程序設(shè)計(jì)3.1 選單和主窗口設(shè)計(jì)流量估計(jì)方法首先依照表1所給數(shù)據(jù)，用matlab作出時(shí)間水位散

12、點(diǎn)圖（圖1）。下面來(lái)計(jì)算水箱流量與時(shí)間的關(guān)系。 根據(jù)圖1，一種簡(jiǎn)單的處理方法為，將表1中的數(shù)據(jù)分為三段，然后對(duì)每一段的數(shù)據(jù)做如下處理： 設(shè)某段數(shù)據(jù) ，相鄰數(shù)據(jù)中點(diǎn)的平均流速用下面的公式（流速=（右端點(diǎn)的水位左端點(diǎn)的水位）/區(qū)間長(zhǎng)度）： 每段數(shù)據(jù)首尾點(diǎn)的流速用下面的公式計(jì)算： 用以上公式求得時(shí)間與流速之間的數(shù)據(jù)如表2。表2時(shí)間/h流速/cmh-1時(shí)間/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水泵開動(dòng)9.98水泵開動(dòng)22.96水泵開動(dòng)10.93水泵開動(dòng)23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表2作出時(shí)間流速散點(diǎn)圖如圖2。插值法由表2，對(duì)水泵不工作時(shí)段1，2采取插值方法，可以得到任意時(shí)刻的流速，從而可以知道任意時(shí)刻的流量。我們分別采取拉格朗日插值法，分段線性插值法及三次樣條插值法；對(duì)于水泵工作時(shí)段1應(yīng)用前后時(shí)期的流速進(jìn)行插值，

14、由于最后一段水泵不工作時(shí)段數(shù)據(jù)太少，我們將它忽略，只對(duì)水泵工作時(shí)段2進(jìn)行插值處理。我們總共需要對(duì)四段數(shù)據(jù)（第1，2未供水時(shí)段，第1供水時(shí)段，混合時(shí)段）進(jìn)行插值處理，下面以第1未供水時(shí)段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量（用水高度）。 調(diào)用程序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);

15、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)實(shí)現(xiàn)圖3結(jié)果 運(yùn)行結(jié)果：lglrjf = 145.6250 fdxxjf =147.1469 sancytjf =145.6885圖中曲線lglr、

16、fdxx和scyt分別表示用拉格朗日插值法，分段線性插值法及三次樣條插值法得到的曲線。 由表1知，第1未供水時(shí)段的總用水高度為146(=968822)，可見上述三種插值方法計(jì)算的結(jié)果與實(shí)際值（146）相比都比較接近?？紤]到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。其他三段的處理方法與第1未供水時(shí)段的處理方法類似.3.2 模塊設(shè)計(jì)3.2.1 四個(gè)主要階段t-v代碼第一供水段時(shí)間流速示意圖程序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 =

17、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運(yùn)行結(jié)果：lglrjf =56.4426 fdxxjf =49.6051 sancytjf =53.5903

18、第2未供水段時(shí)間流速示意圖程序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 =

19、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運(yùn)行結(jié)果：lglrjf = 258.8664 fdxxjf = 258.9697 sancytjf = 258.6547第2供水段時(shí)間流速示意圖程序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運(yùn)行結(jié)果：lglrjf =104.1526fdx

21、xjf =73.7929sancytjf =79.8172下圖是用分段線性及三次樣條插值方法得到的整個(gè)過(guò)程的時(shí)間流速函數(shù)示意圖。程序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運(yùn)行結(jié)果：fdxxj

23、f = 534.4311 sancytjf = 541.8148表3 各時(shí)段及一天的總用水量（用水高度(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 各時(shí)段及一天的總用水量（用水體積(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是對(duì)一天中任取的4個(gè)時(shí)刻分別用3種方法得到的水塔水流量近似值。時(shí)間(h)6.8810.8815.8822.8815.982633.742625.566234.709914.827232.997625.446525.471515.052733.708925.

25、549029.4173注：拉格朗日插值法分段線性插值法三次樣條插值法3.2.2 實(shí)驗(yàn)?zāi)Ｐ蜋z驗(yàn)1、水泵未工作時(shí)的檢驗(yàn)首先計(jì)算水泵未工作的三段時(shí)間內(nèi)的實(shí)際用水量第一時(shí)間段，即t0，8.97時(shí)的實(shí)際用水量為：233.3475(9.68-8.22)=340.69第二時(shí)間段，即t10.95，20.84時(shí)的實(shí)際用水量為：233.3475(10.82-8.22)= 606.70第三時(shí)間段，即t23.88，25.91時(shí)的實(shí)際用水量為：233.3475(10.59-10.18)=95.67其次用三次樣條插值模型得到的函數(shù)分別在水泵未工作的三段時(shí)間內(nèi)作數(shù)值積分，執(zhí)行程序后，可得水泵未工作的三段時(shí)間內(nèi)，用模型計(jì)算出的用水量（見表4-2）。分三段的實(shí)際用水量與模型用水量的比較可見，用三次樣條插值模型計(jì)算的用水量與實(shí)際用水量吻合的很好。2、水泵工作時(shí)的檢驗(yàn)水泵充水量=充水后的水量+充水期間的流出量-充水前的水量。首先利用模型計(jì)算兩次充水期間的流出量，執(zhí)行后可得：第一次充水期間的用水量為117.95，第二次充水期間的用水量為179.06。進(jìn)一步可計(jì)算第一次水泵充水量=233.347510.82+117.95-233.34758.22=724.6535第二次水泵充水量=233