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

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

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

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

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

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

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

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

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

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

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

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

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

14、由于最后一段水泵不工作時段數(shù)據(jù)太少，我們將它忽略，只對水泵工作時段2進(jìn)行插值處理。我們總共需要對四段數(shù)據(jù)（第1，2未供水時段，第1供水時段，混合時段）進(jìn)行插值處理，下面以第1未供水時段數(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未供水時段的總用水高度為146(=968822)，可見上述三種插值方法計算的結(jié)果與實(shí)際值（146）相比都比較接近。考慮到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。其他三段的處理方法與第1未供水時段的處理方法類似.3.2 模塊設(shè)計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 =

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未供水段時間流速示意圖程序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供水段時間流速示意圖程序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下圖是用分段線性及三次樣條插值方法得到的整個過程的時間流速函數(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 各時段及一天的總用水量（用水高度(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.

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