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

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

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

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

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

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

6、平均選取11個(gè)點(diǎn)作插值個(gè)點(diǎn)作插值(xch12)3.在在-6,6中平均選取中平均選取21個(gè)點(diǎn)作插值個(gè)點(diǎn)作插值(xch13)12比分段線性插值更光滑。比分段線性插值更光滑。xyxi-1 xiab 在數(shù)學(xué)上，光滑程度的定量描述是：函數(shù)(曲線)的k階導(dǎo)數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。 光滑性的階次越高，則越光滑。是否存在較低次的分段多項(xiàng)式達(dá)到較高階光滑性的方法？三次樣條插值就是一個(gè)很好的例子。三次樣條插值三次樣條插值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個(gè)基點(diǎn)計(jì)算插值個(gè)基點(diǎn)計(jì)算插值(ych)返回返回To MATLABych(larg1)15用用MATLABMATLAB作插值計(jì)算作插值計(jì)算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xixi

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

9、，2727，2424。試估計(jì)每隔。試估計(jì)每隔1/101/10小時(shí)的小時(shí)的溫度值。溫度值。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時(shí)的時(shí)的y值。值。To MATLAB(plane)返回返回18二維插值的定義二維插值的定義 xyO O第一種（網(wǎng)格節(jié)點(diǎn)）：第一種（網(wǎng)格節(jié)點(diǎn)）：19 已知已知 m n個(gè)節(jié)點(diǎn)個(gè)節(jié)點(diǎn) ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨設(shè)互不相同，不妨設(shè)bxxxam 21dyyycn 21 構(gòu)造一個(gè)二元函數(shù)構(gòu)造一個(gè)二元函數(shù)),(yxfz 通過(guò)全部已知節(jié)點(diǎn)通過(guò)全部已知節(jié)點(diǎn),即即再用再用),(yxf計(jì)算插值，即計(jì)算插值，即).,(

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

12、 y2)O O 二維或高維情形的最鄰近插值，與被插值點(diǎn)最鄰近的節(jié)點(diǎn)的函數(shù)值即為所求。返回返回23 將四個(gè)插值點(diǎn)（矩形的四個(gè)頂點(diǎn)）處的函數(shù)值依次簡(jiǎn)記為： 分片線性插值分片線性插值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é)點(diǎn)所形成的矩形區(qū)域內(nèi)。顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達(dá)式如下：第一片(下三角形區(qū)域)： (x, y)滿足返回返回25 雙線性插值是一片一片的空間二次曲面構(gòu)成。雙線性插值函數(shù)的形式如下：)dcy)(bax()y, x(f其中有四個(gè)待定系數(shù)，利用該函數(shù)在矩形的四個(gè)頂點(diǎn)（插值節(jié)點(diǎn)）的函數(shù)值，得到四個(gè)代數(shù)方程，正好確定四個(gè)系數(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可取可取為矩陣，或?yàn)榫仃?，或x x取取行向量，行向量，y y取為列向量，取為列向量，x,yx,y的值分別不能超出的值分別不能超出x0,y0 x0,y0的范圍。的范圍。z=interp2(x0,y0,z0,x,y,method)被插值點(diǎn)插值方法用用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值nearestnearest 最鄰近插值最鄰近插值linearlinear 雙線性插值雙線性插值cubiccubic 雙三次插值雙三次插值缺省時(shí)缺省時(shí), , 雙線性插值雙線性插值27例：測(cè)得平板表面例：測(cè)得平板表面3 3* *5 5網(wǎng)格點(diǎn)處的溫度分別為：網(wǎng)格點(diǎn)處的溫度

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個(gè)單位的地方進(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ū)測(cè)得一些地點(diǎn)的高程如下表。平面區(qū)域?yàn)樵谀成絽^(qū)測(cè)得一些地點(diǎn)的高程如下表。平面區(qū)域?yàn)?1200=x=4000,1200=y=3600)試作出該山區(qū)的地貌圖和等高線圖，并對(duì)幾種插值方法進(jìn)行比較。試作出該山區(qū)的地貌圖和等高線圖，并對(duì)幾種插值方法進(jìn)行比較。 X Y12001600200024002800320036004000120011301250128012301

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

18、ta格式為格式為: cz =griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算作散點(diǎn)數(shù)據(jù)的插值計(jì)算 要求要求cxcx取行向量，取行向量，cycy取為列向量取為列向量。被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值nearestnearest 最鄰近插值最鄰近插值linearlinear 雙線性插值雙線性插值cubiccubic 雙三次插值雙三次插值v4- Matlab提供的插值方法提供的插值方法缺省時(shí)缺省時(shí), , 雙線性插值雙線性插值31 例例 在某海域測(cè)得一些點(diǎn)在某海域測(cè)得一些點(diǎn)(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、 輸入插值基點(diǎn)數(shù)據(jù)To MATLAB hd1返回返回4.作出水深小于5的海域范圍,即z=5的等高線.33實(shí)驗(yàn)作業(yè)實(shí)驗(yàn)作業(yè) 山區(qū)地貌：山區(qū)地貌：在某山區(qū)測(cè)得一些地點(diǎn)的高程如下表：在某山區(qū)測(cè)得一些地點(diǎn)的高程如下表：( (平平面區(qū)域面區(qū)域1200=x=4000,1200=y=3600)1200=x=4000,1200=y=3600)，試作出該山區(qū)的試作出該山區(qū)的地貌圖和等高線圖，并對(duì)幾種插值方法進(jìn)行比較。地貌圖和等高線圖，并對(duì)幾種插值方法進(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í)驗(yàn)數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)后勤工程學(xué)院數(shù)學(xué)教研室 擬擬 合合 3

22、5實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)?zāi)康膶?shí)驗(yàn)內(nèi)容實(shí)驗(yàn)內(nèi)容2 2、掌握用數(shù)學(xué)軟件求解擬合問(wèn)題。、掌握用數(shù)學(xué)軟件求解擬合問(wèn)題。1 1、直觀了解擬合基本內(nèi)容。、直觀了解擬合基本內(nèi)容。1 1、擬合問(wèn)題引例及基本、擬合問(wèn)題引例及基本理論。理論。4 4、實(shí)驗(yàn)作業(yè)。、實(shí)驗(yàn)作業(yè)。2 2、用數(shù)學(xué)軟件求解擬合問(wèn)題。、用數(shù)學(xué)軟件求解擬合問(wèn)題。3 3、應(yīng)用實(shí)例、應(yīng)用實(shí)例36擬擬 合合2.2.擬合的基本原理擬合的基本原理1. 擬合問(wèn)題引例擬合問(wèn)題引例37擬擬 合合 問(wèn)問(wèn) 題題 引引 例例 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時(shí)的電阻時(shí)的電阻R R。2040608010070080090010001100 設(shè)設(shè) R=at+bR=at+ba,ba,b為待定系為待定系數(shù)數(shù)38擬擬 合合 問(wèn)問(wèn) 題題 引引 例例 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)求血藥濃度隨時(shí)間的變化規(guī)律求血藥濃度隨時(shí)間的變化規(guī)律c(t).c(t).作半對(duì)數(shù)坐標(biāo)系作半對(duì)數(shù)坐標(biāo)系(semilogy)下的下的圖形圖形為待定系數(shù)kcectckt,)(002468100101102MATLAB(aa1)MATLAB(aa1)39曲曲 線線 擬擬 合合 問(wèn)問(wèn) 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n n個(gè)點(diǎn)個(gè)點(diǎn)（x

25、 xi i,y,yi i) ) i=1,n,i=1,n, 尋求一個(gè)函數(shù)（曲線）尋求一個(gè)函數(shù)（曲線）y=f(x),y=f(x), 使使 f(x)f(x) 在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。得最好。 +xyy=f(x)y=f(x)(xi,yi)i i i 為點(diǎn)為點(diǎn)（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)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)一個(gè)函數(shù)作為近似，由于近

26、似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。學(xué)方法上是完全不同的。 實(shí)例：實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到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)問(wèn)題：?jiǎn)栴}：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：解決方案：若不要求曲線（面）通過(guò)所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，這就是數(shù)據(jù)擬合數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過(guò)所給所有數(shù)據(jù)點(diǎn)，就是插值問(wèn)題插值問(wèn)題；41最臨近插值、線性插值、樣條插值與曲線

27、擬最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：合結(jié)果：0246810121416180510152025已知數(shù)據(jù)點(diǎn)spline三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)linest三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)nearest三次多項(xiàng)式插值42曲線擬合問(wèn)題最常用的解法曲線擬合問(wèn)題最常用的解法線性最小二乘法的線性最小二乘法的基本思路基本思路第一步: :先選定一組函數(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、體看作一個(gè)房室，室內(nèi)血藥濃度均勻機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻一室模型一室模型模型建立模型建立d/c(0) 3得：由假設(shè)-kcdtdc 2得：由假設(shè)ktevdtc)( 在此，在此，d=300mgd=300mg，t t及及c c（t t）在某些點(diǎn)處的）在某些點(diǎn)處的值見(jiàn)前表，需經(jīng)擬合求出參數(shù)值見(jiàn)前表，需經(jīng)擬合求出參數(shù)k k、v v用線性最小二乘擬合用線性最小二乘擬合c(t)c(t)ktevdtc)()/ln(,ln21vdakacyktvdc)/ln(ln2/,121aedvakatayMATLAB(lihe1)MATLAB(lihe1)計(jì)算結(jié)果：計(jì)算結(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è)計(jì)計(jì)cc2c10t 設(shè)每次注射劑量D, 間隔時(shí)間 血藥濃度c(t) 應(yīng)c1 c(t) c2 初次劑量D0 應(yīng)加大,0DD給藥方案記為：給藥方案記為：kecc2112ln1cck2 2、)( ,1220ccDcD1 1、計(jì)算結(jié)果：計(jì)算結(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， 注射的間隔時(shí)間為注射的間隔時(shí)間為4 4小時(shí)。小時(shí)。69估計(jì)水塔的流量估計(jì)水塔的流量2 2、解題思路、解題思路3 3、算法設(shè)計(jì)與編程、算法設(shè)計(jì)與編程1 1、問(wèn)題、問(wèn)題70 某居民區(qū)有一供居民用水

31、的園柱形水塔，一般可以通過(guò)測(cè)量其水位來(lái)估計(jì)水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時(shí)，水泵自動(dòng)啟動(dòng)向水塔供水，到設(shè)定的最高水位時(shí)停止供水，這段時(shí)間無(wú)法測(cè)量水塔的水位和水泵的供水量通常水泵每天供水一兩次，每次約兩小時(shí).水塔是一個(gè)高12.2米，直徑17.4米的正園柱按照設(shè)計(jì)，水塔水位降至約8.2米時(shí)，水泵自動(dòng)啟動(dòng)，水位升到約10.8米時(shí)水泵停止工作表1 是某一天的水位測(cè)量記錄，試估計(jì)任何時(shí)刻（包括水泵正供水時(shí)）從水塔流出的水流量，及一天的總用水量71 表 1 水位測(cè)量記錄 （符號(hào)/表示水泵啟動(dòng)）時(shí)刻(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時(shí)刻(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時(shí)刻(h)水位(cm)19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91866 843 822 / / 1059 1035 101872流量估計(jì)的解題思路流量估計(jì)的解題思路73 擬合水位擬合水位 時(shí)間函數(shù)時(shí)間函數(shù) 測(cè)量記錄看，一天有兩個(gè)供

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

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

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

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

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

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

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

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