數(shù)學(xué)建模里MATLAB 插值與擬合

1、1數(shù)學(xué)建模與數(shù)學(xué)實驗數(shù)學(xué)建模與數(shù)學(xué)實驗后勤工程學(xué)院數(shù)學(xué)教研室 插插 值值2實驗?zāi)康膶嶒災(zāi)康膶嶒瀮?nèi)容實驗內(nèi)容2、掌握用數(shù)學(xué)軟件包求解插值問題。、掌握用數(shù)學(xué)軟件包求解插值問題。1、了解插值的基本內(nèi)容。、了解插值的基本內(nèi)容。11一維插值一維插值22二維插值二維插值33實驗作業(yè)實驗作業(yè)3一一 維維 插插 值值一、插值的定義一、插值的定義二、插值的方法二、插值的方法三、用三、用Matlab解插值問題解插值問題返回返回4返回返回二維插值二維插值一、二維插值定義一、二維插值定義二、網(wǎng)格節(jié)點插值法二、網(wǎng)格節(jié)點插值法三、用三、用MatlabMatlab解插值問題解插值問題最鄰近插值最鄰近插值分片線性插值分片線性

2、插值雙線性插值雙線性插值網(wǎng)格節(jié)點數(shù)據(jù)的插值網(wǎng)格節(jié)點數(shù)據(jù)的插值散點數(shù)據(jù)的插值散點數(shù)據(jù)的插值5一維插值的定義一維插值的定義已知已知 n+1個節(jié)點個節(jié)點, 1 , 0(),(njyxjj其中其中jx互不相同，不妨設(shè)互不相同，不妨設(shè)),10bxxxan求任一插值點求任一插值點)(*jxx 處的插值處的插值.*y0 x1xnx0y1y節(jié)點可視為由節(jié)點可視為由)(xgy 產(chǎn)生產(chǎn)生,，g表達(dá)式復(fù)雜表達(dá)式復(fù)雜,，或無封閉形式或無封閉形式,，或未知或未知.。*x*y6 構(gòu)造一個構(gòu)造一個(相對簡單的相對簡單的)函數(shù)函數(shù)),(xfy 通過全部節(jié)點通過全部節(jié)點, 即即), 1 ,0()(njyxfjj再用再用)(xf

3、計算插值，即計算插值，即).(*xfy 0 x1xnx0y1y*x*y返回返回7 稱為拉格朗日插值基函數(shù)拉格朗日插值基函數(shù)。n0iiiny)x(L)x(P 已知函數(shù)f(x)在n+1個點x0,x1,xn處的函數(shù)值為 y0,y1,yn 。求一n次多項式函數(shù)Pn(x)，使其滿足： Pn(xi)=yi,i=0,1,n. 解決此問題的拉格朗日插值多項式公式如下其中Li(x) 為n次多項式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i拉格朗日拉格朗日(Lagrange)插值插值8拉格朗日拉格朗日(Lagrange)插值插值

4、特別地特別地:兩點一次兩點一次(線性線性)插值多項式插值多項式: 101001011yxxxxyxxxxxL三點二次三點二次(拋物拋物)插值多項式插值多項式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL .,滿足插值條件直接驗證可知xLn9 拉格朗日多項式插值的這種振蕩現(xiàn)象叫 Runge現(xiàn)象現(xiàn)象55,11)(2xxxg 采用拉格朗日多項式插值：選取不同插值節(jié)點個數(shù)n+1，其中n為插值多項式的次數(shù)，當(dāng)n分別取2,4,6,8,10時，繪出插值結(jié)果圖形.例例返回返回To MatlabTo Matlablch(larg1)lch(larg1)

5、10分段線性插值分段線性插值其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL計算量與n無關(guān);n越大，誤差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy11To MATLABxch11，xch12，xch13，xch14返回返回66,11)(2xxxg例例用分段線性插值法求插值用分段線性插值法求插值,并觀察插值誤差并觀察插值誤差.1.在在-6,6中平均選取中平均選取5個點作插值個點作插值(xch11)4.在在-6,6中平均選取中平均選取41個點作插值個點作插值(xch14)2.在在-6,6中平均選取中

6、平均選取11個點作插值個點作插值(xch12)3.在在-6,6中平均選取中平均選取21個點作插值個點作插值(xch13)12比分段線性插值更光滑。比分段線性插值更光滑。xyxi-1 xiab 在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。 光滑性的階次越高，則越光滑。是否存在較低次的分段多項式達(dá)到較高階光滑性的方法？三次樣條插值就是一個很好的例子。三次樣條插值三次樣條插值13 三次樣條插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs

7、) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然邊界條件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii)()(limxgxSng g( (x x) )為被插值函數(shù)為被插值函數(shù)。14例例66,11)(2xxxg用三次樣條插值選取用三次樣條插值選取11個基點計算插值個基點計算插值(ych)返回返回To MATLABych(larg1)15用用MATLABMATLAB作插值計算作插值計算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點被插值點插值節(jié)點插值節(jié)點xixi

8、處的插處的插值結(jié)果值結(jié)果nearest ：最鄰近插值：最鄰近插值linear ： 線性插值；線性插值；spline ： 三次樣條插三次樣條插值；值；cubic ： 立方插值。立方插值。缺省時：缺省時： 分段線性插值。分段線性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是單調(diào)的，并且是單調(diào)的，并且xi不不能夠超過能夠超過x的范圍。的范圍。16 例：在例：在1-121-12的的1111小時內(nèi)，每隔小時內(nèi)，每隔1 1小時測量一次小時測量一次溫度，測得的溫度依次為：溫度，測得的溫度依次為：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525

9、，2727，2424。試估計每隔。試估計每隔1/101/10小時的小時的溫度值。溫度值。To MATLAB(temp)hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); (直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作圖xlabel(Hour),ylabel(Degrees Celsius)17xy機(jī)翼下輪廓線X035791 11 21 31 41 5Y01 . 21 . 72 . 02 . 12 .

10、 01 . 81 . 21 . 01 . 6例例 已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變每改變0.1時的時的y值。值。To MATLAB(plane)返回返回18二維插值的定義二維插值的定義 xyO O第一種（網(wǎng)格節(jié)點）：第一種（網(wǎng)格節(jié)點）：19 已知已知 m n個節(jié)點個節(jié)點 ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨設(shè)互不相同，不妨設(shè)bxxxam 21dyyycn 21 構(gòu)造一個二元函數(shù)構(gòu)造一個二元函數(shù)),(yxfz 通過全部已知節(jié)點通過全部已知節(jié)點,即即再用再用),(yxf計算插值，即計算插值，即).,(

11、*yxfz ),1 ,0;,1 ,0(),(njmizyxfijji 20第二種（散亂節(jié)點）：第二種（散亂節(jié)點）： yx0 021已知已知n個節(jié)點個節(jié)點),.,2 , 1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同， 構(gòu)造一個二元函數(shù)構(gòu)造一個二元函數(shù)),(yxfz 通過全部已知節(jié)點通過全部已知節(jié)點,即即),1 ,0(),(nizyxfiii 再用再用),(yxf計算插值，即計算插值，即).,(*yxfz 返回返回22 注意：注意：最鄰近插值一般不連續(xù)。具有連續(xù)性的最簡單的插值是分片線性插值。最鄰近插值最鄰近插值x y(x1, y1)(x1, y2)(x2, y1)(x2,

12、 y2)O O 二維或高維情形的最鄰近插值，與被插值點最鄰近的節(jié)點的函數(shù)值即為所求。返回返回23 將四個插值點（矩形的四個頂點）處的函數(shù)值依次簡記為： 分片線性插值分片線性插值xy (xi, yj)(xi, yj+1)(xi+1, yj)(xi+1, yj+1)O Of (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f424插值函數(shù)為：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x(fj23i121第二片(上三角形區(qū)域)：(x, y)滿足iii1ij1jy)xx(xxyyy插值函數(shù)為：

13、)xx)(ff ()yy)(ff (f)y, x(fi43j141注意注意：(x, y)當(dāng)然應(yīng)該是在插值節(jié)點所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)： (x, y)滿足返回返回25 雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：)dcy)(bax()y, x(f其中有四個待定系數(shù)，利用該函數(shù)在矩形的四個頂點（插值節(jié)點）的函數(shù)值，得到四個代數(shù)方程，正好確定四個系數(shù)。雙線性插值雙線性插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O返回返回26 要求要求x0,y0 x0,y0單調(diào)；單調(diào)；x

14、x，y y可取可取為矩陣，或為矩陣，或x x取取行向量，行向量，y y取為列向量，取為列向量，x,yx,y的值分別不能超出的值分別不能超出x0,y0 x0,y0的范圍。的范圍。z=interp2(x0,y0,z0,x,y,method)被插值點插值方法用用MATLAB作網(wǎng)格節(jié)點數(shù)據(jù)的插值作網(wǎng)格節(jié)點數(shù)據(jù)的插值插值節(jié)點被插值點的函數(shù)值nearestnearest 最鄰近插值最鄰近插值linearlinear 雙線性插值雙線性插值cubiccubic 雙三次插值雙三次插值缺省時缺省時, , 雙線性插值雙線性插值27例：測得平板表面例：測得平板表面3 3* *5 5網(wǎng)格點處的溫度分別為：網(wǎng)格點處的溫度

15、分別為： 82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 試作出平板表面的溫度分布曲面試作出平板表面的溫度分布曲面z=f(x,y)z=f(x,y)的圖形。的圖形。輸入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.2以平滑數(shù)據(jù),在x、y方向上每隔0.2個單位的地方進(jìn)行插值.28再輸入以下命令:x

16、i=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖. To MATLAB(wendu)29例例 山區(qū)地貌：山區(qū)地貌： 在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為在某山區(qū)測得一些地點的高程如下表。平面區(qū)域為 1200=x=4000,1200=y=3600)試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進(jìn)行比較。試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進(jìn)行比較。 X Y12001600200024002800320036004000120011301250128012301

17、040900500700160013201450142014001300700900850200013901500150014009001100106095024001500120011001350145012001150101028001500120011001550160015501380107032001500155016001550160016001600155036001480150015501510143013001200980 通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較。To MATLAB (moutain)返回返回30 插值函數(shù)插值函數(shù)gridda

18、ta格式為格式為: cz =griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散點數(shù)據(jù)的插值計算作散點數(shù)據(jù)的插值計算 要求要求cxcx取行向量，取行向量，cycy取為列向量取為列向量。被插值點插值方法插值節(jié)點被插值點的函數(shù)值nearestnearest 最鄰近插值最鄰近插值linearlinear 雙線性插值雙線性插值cubiccubic 雙三次插值雙三次插值v4- Matlab提供的插值方法提供的插值方法缺省時缺省時, , 雙線性插值雙線性插值31 例例 在某海域測得一些點在某海域測得一些點(x,y)(x,y)處的水深處的水深z z由下由下表給出，船的吃水

19、深度為表給出，船的吃水深度為5 5英尺，在矩形區(qū)域（英尺，在矩形區(qū)域（7575，200200）* *（-50-50，150150）里的哪些地方船要避免進(jìn)入。）里的哪些地方船要避免進(jìn)入。xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.59 9 8 8 9 4 932 ) 1( .150,50200,75. 2hd三次插值法作二維插值在矩形區(qū)域. 3 作海底曲面圖.1

20、 輸入插值基點數(shù)據(jù)To MATLAB hd1返回返回4.作出水深小于5的海域范圍,即z=5的等高線.33實驗作業(yè)實驗作業(yè) 山區(qū)地貌：山區(qū)地貌：在某山區(qū)測得一些地點的高程如下表：在某山區(qū)測得一些地點的高程如下表：( (平平面區(qū)域面區(qū)域1200=x=4000,1200=y=3600)1200=x=4000,1200=y=3600)，試作出該山區(qū)的試作出該山區(qū)的地貌圖和等高線圖，并對幾種插值方法進(jìn)行比較。地貌圖和等高線圖，并對幾種插值方法進(jìn)行比較。36003200280024002000160012001480 1500 1550 1510 1430 1300 1200 9801500 1550 1

21、600 1550 1600 1600 1600 15501500 1200 1100 1550 1600 1550 1380 10701500 1200 1100 1350 1450 1200 1150 10101390 1500 1500 1400 900 1100 1060 9501320 1450 1420 1400 1300 700 900 8501130 1250 1280 1230 1040 900 500 700Y/x1200 1600 2000 2400 2800 3200 3600 4000返回返回34數(shù)學(xué)建模與數(shù)學(xué)實驗數(shù)學(xué)建模與數(shù)學(xué)實驗后勤工程學(xué)院數(shù)學(xué)教研室 擬擬 合合 3

22、5實驗?zāi)康膶嶒災(zāi)康膶嶒瀮?nèi)容實驗內(nèi)容2 2、掌握用數(shù)學(xué)軟件求解擬合問題。、掌握用數(shù)學(xué)軟件求解擬合問題。1 1、直觀了解擬合基本內(nèi)容。、直觀了解擬合基本內(nèi)容。1 1、擬合問題引例及基本、擬合問題引例及基本理論。理論。4 4、實驗作業(yè)。、實驗作業(yè)。2 2、用數(shù)學(xué)軟件求解擬合問題。、用數(shù)學(xué)軟件求解擬合問題。3 3、應(yīng)用實例、應(yīng)用實例36擬擬 合合2.2.擬合的基本原理擬合的基本原理1. 擬合問題引例擬合問題引例37擬擬 合合 問問 題題 引引 例例 1 1溫度溫度t(t(0 0C) 20.5 32.7 C) 20.5 32.7 51.0 73.0 95.751.0 73.0 95.7電阻電阻R(R(

23、) 765 826 ) 765 826 873 942 1032873 942 1032已知熱敏電阻數(shù)已知熱敏電阻數(shù)據(jù)：據(jù)：求求60600 0C C時的電阻時的電阻R R。2040608010070080090010001100 設(shè)設(shè) R=at+bR=at+ba,ba,b為待定系為待定系數(shù)數(shù)38擬擬 合合 問問 題題 引引 例例 2 2 t (h) 0.25 0.5 1 t (h) 0.25 0.5 1 1.5 2 3 4 6 81.5 2 3 4 6 8c (c ( g/ml) g/ml) 19.21 18.15 15.36 14.10 12.89 19.21 18.15 15.36 14.

24、10 12.89 9.32 7.45 5.24 3.019.32 7.45 5.24 3.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0(t=0注注射射300300mg)mg)求血藥濃度隨時間的變化規(guī)律求血藥濃度隨時間的變化規(guī)律c(t).c(t).作半對數(shù)坐標(biāo)系作半對數(shù)坐標(biāo)系(semilogy)下的下的圖形圖形為待定系數(shù)kcectckt,)(002468100101102MATLAB(aa1)MATLAB(aa1)39曲曲 線線 擬擬 合合 問問 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n n個點個點（x

25、 xi i,y,yi i) ) i=1,n,i=1,n, 尋求一個函數(shù)（曲線）尋求一個函數(shù)（曲線）y=f(x),y=f(x), 使使 f(x)f(x) 在某種準(zhǔn)則下與所有數(shù)據(jù)點最為接近，即曲線擬合在某種準(zhǔn)則下與所有數(shù)據(jù)點最為接近，即曲線擬合得最好。得最好。 +xyy=f(x)y=f(x)(xi,yi)i i i 為點為點（x xi i,y,yi i) ) 與與曲線曲線 y=f(x) y=f(x) 的距離的距離40擬合與插值的關(guān)系擬合與插值的關(guān)系 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似，由于近似的要求不同，二者的數(shù)一個函數(shù)作為近似，由于近

26、似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。學(xué)方法上是完全不同的。 實例：實例：下面數(shù)據(jù)是某次實驗所得，希望得到X和 f之間的關(guān)系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)MATLAB(cn)問題：問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面解決方案：解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點，就是插值問題插值問題；41最臨近插值、線性插值、樣條插值與曲線

27、擬最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：合結(jié)果：0246810121416180510152025已知數(shù)據(jù)點spline三次多項式插值0246810121416180510152025已知數(shù)據(jù)點linest三次多項式插值0246810121416180510152025已知數(shù)據(jù)點nearest三次多項式插值42曲線擬合問題最常用的解法曲線擬合問題最常用的解法線性最小二乘法的線性最小二乘法的基本思路基本思路第一步: :先選定一組函數(shù)先選定一組函數(shù) r r1 1(x), r(x), r2 2(x), r(x), rm m(x), (x), mn,m0)k(0)模型假設(shè)模型假設(shè)1. 1. 機(jī)

28、體看作一個房室，室內(nèi)血藥濃度均勻機(jī)體看作一個房室，室內(nèi)血藥濃度均勻一室模型一室模型模型建立模型建立d/c(0) 3得：由假設(shè)-kcdtdc 2得：由假設(shè)ktevdtc)( 在此，在此，d=300mgd=300mg，t t及及c c（t t）在某些點處的）在某些點處的值見前表，需經(jīng)擬合求出參數(shù)值見前表，需經(jīng)擬合求出參數(shù)k k、v v用線性最小二乘擬合用線性最小二乘擬合c(t)c(t)ktevdtc)()/ln(,ln21vdakacyktvdc)/ln(ln2/,121aedvakatayMATLAB(lihe1)MATLAB(lihe1)計算結(jié)果：計算結(jié)果：)(02.15),/1 (2347.

29、 0lvhkd=300;t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2)程序：程序：用非線性最用非線性最小二乘擬合小二乘擬合c(t)c(t)給藥方案給藥方案 設(shè)設(shè)計計cc2c10t 設(shè)每次注射劑量D, 間隔時間 血藥濃度c(t) 應(yīng)c1 c(t) c2 初次劑量D0 應(yīng)加大,0DD給藥方案記為：給藥方案記為：kecc2112ln1cck2 2、)( ,1220ccDcD1 1、計算結(jié)果：計算結(jié)果：9

30、 . 3, 3 .225, 5 .3750DD)(4),(225),(3750hmgDmgD給藥方案：給藥方案：c c1 1=10,c=10,c2 2=25=25k=0.2347k=0.2347v=15.02v=15.0268故可制定給藥方案：故可制定給藥方案：)(4),(225),(3750hmgDmgD即即: : 首次注射首次注射375mg375mg， 其余每次注射其余每次注射225mg225mg， 注射的間隔時間為注射的間隔時間為4 4小時。小時。69估計水塔的流量估計水塔的流量2 2、解題思路、解題思路3 3、算法設(shè)計與編程、算法設(shè)計與編程1 1、問題、問題70 某居民區(qū)有一供居民用水

31、的園柱形水塔，一般可以通過測量其水位來估計水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時，水泵自動啟動向水塔供水，到設(shè)定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量通常水泵每天供水一兩次，每次約兩小時.水塔是一個高12.2米，直徑17.4米的正園柱按照設(shè)計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米時水泵停止工作表1 是某一天的水位測量記錄，試估計任何時刻（包括水泵正供水時）從水塔流出的水流量，及一天的總用水量71 表 1 水位測量記錄 （符號/表示水泵啟動）時刻(h)水位(cm)0 0.92 1.84 2.95 3.87 4.98 5.90 7

32、.01 7.93 8.97968 948 931 913 898 881 869 852 839 822時刻(h)水位(cm)9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93/ / 1082 1050 1021 994 965 941 918 892時刻(h)水位(cm)19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91866 843 822 / / 1059 1035 101872流量估計的解題思路流量估計的解題思路73 擬合水位擬合水位 時間函數(shù)時間函數(shù) 測量記錄看，一天有兩個供

33、水時段（以下稱第1供水時段和第2供水時段），和3個水泵不工作時段（以下稱第1時段t=0到t=8.97，第2次時段t=10.95到t=20.84和第3時段t=23以后）對第1、2時段的測量數(shù)據(jù)直接分別作多項式擬合，得到水位函數(shù)為使擬合曲線比較光滑，多項式次數(shù)不要太高，一般在36由于第3時段只有3個測量記錄，無法對這一時段的水位作出較好的擬合74 2、確定流量確定流量 時時間函數(shù)間函數(shù) 對于第1、2時段只需將水位函數(shù)求導(dǎo)數(shù)即可，對于兩個供水時段的流量，則用供水時段前后（水泵不工作時段）的流量擬合得到，并且將擬合得到的第2供水時段流量外推，將第3時段流量包含在第2供水時段內(nèi)753、一天總用水量的估計

34、一天總用水量的估計 總用水量等于兩個水泵不工作時段和兩個供水時段用水量之和，它們都可以由流量對時間的積分得到。76算法設(shè)計與編程算法設(shè)計與編程1 1、擬合第、擬合第1 1、2 2時段的水位，并導(dǎo)出流量時段的水位，并導(dǎo)出流量2 2、擬合供水時段的流量、擬合供水時段的流量3 3、估計一天總用水量、估計一天總用水量4 4、流量及總用水量的檢驗、流量及總用水量的檢驗77 1、擬合第擬合第1 1時段的水位，并導(dǎo)出流時段的水位，并導(dǎo)出流量量 設(shè)t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第第1 1時段時段各時刻的流量可如下得：1） c1=polyfitc1=polyfit（t t（1

35、1：1010），），h h（1 1：1010），），3 3）；）； %用3次多項式擬合第1時段水位，c1輸出3次多項式的系數(shù)2）a1=polydera1=polyder（c1c1）；）； % a1輸出多項式（系數(shù)為c1）導(dǎo)數(shù)的系數(shù) 3）tp1=0tp1=0：0.10.1：9 9； x1=-polyvalx1=-polyval（a1a1，tp1tp1）；）； % x1輸出多項式（系數(shù)為a1）在tp1點的函數(shù)值（取負(fù)后邊為正值），即tp1時刻的流量 MATLAB(llgj1)MATLAB(llgj1)4）流量函數(shù)為：流量函數(shù)為：1079.227173. 22356. 0)(2tttf78 2、擬合

36、第擬合第2 2時段的水位，并導(dǎo)出流時段的水位，并導(dǎo)出流量量 設(shè)t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第第2 2時段時段各時刻的流量可如下得：1） c2=polyfit(t(10.9:21),h(10.9:21),3)； %用3次多項式擬合第2時段水位，c2輸出3次多項式的系數(shù)2） a2=polyder(c2)； % a2輸出多項式（系數(shù)為c2）導(dǎo)數(shù)的系數(shù) 3）tp2=10.9:0.1:21; x2=-polyval(a2,tp2); % x2輸出多項式（系數(shù)為a2）在tp2點的函數(shù)值（取負(fù)后邊為正值），即tp2時刻的流量MATLAB(llgj2)MATLAB(llgj

37、2)4）流量函數(shù)為：流量函數(shù)為：8313. 17512. 87529. 00186. 0)(23ttttf79 3、擬合供水時段的流量擬合供水時段的流量 在第1供水時段（t=911）之前（即第1時段）和之后（即第2時段）各取幾點，其流量已經(jīng)得到，用它們擬合第1供水時段的流量為使流量函數(shù)在t=9和t=11連續(xù)，我們簡單地只取4個點，擬合3次多項式（即曲線必過這4個點），實現(xiàn)如下： xx1=-polyval（a1，8 9）； %取第1時段在t=8,9的流量 xx2=-polyval（a2，11 12）； %取第2時段在t=11,12的流量 xx12=xx1 xx2； c12=polyfit（8 9

38、 11 12，xx12，3）； %擬合3次多項式 tp12=9：0.1：11； x12=polyval（c12，tp12）； % x12輸出第1供水時段 各時刻的流量MATLAB(llgj3)MATLAB(llgj3)擬合的流量函數(shù)為：擬合的流量函數(shù)為：078.3555879.737207. 3)(2tttf80 在第2供水時段之前取t=20，20.8兩點的流水量，在該時刻之后（第3時段）僅有3個水位記錄，我們用差分得到流量，然后用這4個數(shù)值擬合第2供水時段的流量如下： dt3=diff（t(22：24)）； %最后3個時刻的兩兩之差 dh3=diff（h(22：24)）； %最后3個水位的兩

39、兩之差 dht3=-dh3./dt3； %t(22)和t(23)的流量 t3=20 20.8 t(22) t(23)； xx3=-polyval(a2，t3(1：2)，dht3)； %取t3各時刻的流量 c3=polyfit（t3，xx3，3）； %擬合3次多項式 t3=20.8：0.1：24； x3=polyval（c3，tp3）；% x3輸出第2供水時段 （外推至t=24）各時刻的流量MATLAB(llgj4)MATLAB(llgj4)擬合的流量函數(shù)為：擬合的流量函數(shù)為：8283.913077. 71405. 0)(2tttf81 3、一天總用水量的估計一天總用水量的估計 第1、2時段和第

40、1、2供水時段流量的積分之和，就是一天總用水量雖然諸時段的流量已表為多項式函數(shù)，積分可以解析地算出，這里仍用數(shù)值積分計算如下： y1=0.1*trapz(x1)； %第1時段用水量（仍按高 度計），0.1為積分步長 y2=0.1*trapz(x2)； %第2時段用水量 y12=0.1*trapz(x12)； %第1供水時段用水量 y3=0.1*trapz(x3)； %第2供水時段用水量 y=（y1+y2+y12+y3）*237.8*0.01； %一天總用水量（ ） 計算結(jié)果：計算結(jié)果：y1=146.2, y2=266.8, y12=47.4, y3=77.3，y=1250.4Lm3310MATLAB(llgjz)MATLAB(llgjz)82 4、流量及總用水量的檢驗流量及總用水量的檢驗 計算出的各時刻的流量各時刻的流量可用水位記錄的數(shù)值微分來檢驗用水量y1可用第1時段水位測量記錄中下降高度968-822=146來檢驗，類似地，y2用1082-822=260檢驗供水時段流量供水時段流量的一種檢驗方法檢驗方法如下：供水時段的用水量加上水位上升值260是該時段泵入的水量，除以時段長度得到水泵的功率（單位時間泵入的水量），而兩個供水時段水泵的功率