插值和擬合(講義)

1、插值與擬合在數(shù)學(xué)建模中的應(yīng)用插值與擬合在數(shù)學(xué)建模中的應(yīng)用數(shù)學(xué)學(xué)院數(shù)學(xué)學(xué)院 袁栩袁栩11一維插值一維插值22二維插值二維插值注：Hermite插值（略）Runge現(xiàn)象用用MATLAB作插值計(jì)算作插值計(jì)算一維插值函數(shù)：一維插值函數(shù)：yi=interp1(x，y，xi，method)插值方法插值方法被插值點(diǎn)被插值點(diǎn)插值節(jié)點(diǎn)插值節(jié)點(diǎn)xi處的處的插值結(jié)果插值結(jié)果nearest 最鄰近插值；最鄰近插值；linear 線性插值線性插值(缺?。?；缺?。?；spline 三次樣條插值；三次樣條插值；cubic 立方插值；立方插值； 缺省時(shí)缺省時(shí) 分段線性插值分段線性插值 注意：所有的插值注意：所有的插值方法都要求

2、方法都要求x是單調(diào)，是單調(diào)，并且并且xi不能夠超過不能夠超過x的的范圍范圍例例1.55,11)(2xxxg用三次樣條插值選取用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值個(gè)基點(diǎn)計(jì)算插值x0=linspace(-5,5,11);y0=1./(1+x0.2);x=linspace(-5,5,100);y=interp1(x0,y0,x,spline);x1=linspace(-5,5,100);y1=1./(1+x1.2);plot(x1,y1,k,x0,y0,+,x,y,r)此例，可以看出插值函數(shù)得到的函數(shù)圖像與原函數(shù)很接近。此例，可以看出插值函數(shù)得到的函數(shù)圖像與原函數(shù)很接近。 例例2 2：從：從1 1點(diǎn)

3、點(diǎn)1212點(diǎn)點(diǎn)的的1111小時(shí)內(nèi)，每隔小時(shí)內(nèi)，每隔1 1小時(shí)測量小時(shí)測量一次溫度，測得的溫度的數(shù)值依次為：一次溫度，測得的溫度的數(shù)值依次為：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424試估計(jì)每試估計(jì)每隔隔1/101/10小時(shí)的溫度值小時(shí)的溫度值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,b,hours,temps

4、,r:) %作圖xlabel(Hour),ylabel(Degrees Celsius)xy機(jī)翼下輪廓線例例3. 3. 已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變每改變0.1時(shí)的時(shí)的y值值x0=0 3 5 7 9 11 12 13 14 15 ;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ;x=0:0.1:15;y1=interp1(x0,y0,x);y2=interp1(x0,y0,x,spline);subplot(2,1,1)plot(x0,y0,k+,x,y1,r)gridtitle(piecewise linear)su

5、bplot(2,1,2)plot(x0,y0,k+,x,y2,r)gridtitle(spline) 要求要求x0, ,y0單調(diào)；單調(diào)；x，y可取可取為矩陣，或?yàn)榫仃?，或x取行向量，取行向量，y取為列向量，取為列向量，x,y的值分別不能超出的值分別不能超出x0, ,y0 0的范圍的范圍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ù)值nearest 最鄰近插值；最鄰近插值；linear 雙線性插值（缺省值）；雙線性插值（缺省值）；cubic 雙三次插值；雙三次插值；缺省時(shí)缺省時(shí) 雙線性插值

6、雙線性插值. .例例4.4.測得平板表面測得平板表面3 35 5網(wǎng)格點(diǎn)處的溫度分別為：網(wǎng)格點(diǎn)處的溫度分別為： 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) )的圖形的圖形輸入以下命令：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ù)，畫出粗糙的

7、溫度分布曲線圖.2以平滑數(shù)據(jù),在 x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖. 通過此例對(duì)最近鄰點(diǎn)插值、雙線性插值方法和雙三次插值方法的插值效果進(jìn)行比較figure(1);meshz(x,y,z)xlabel(X)ylabel(Y)zlabel(Z)xi=0:50:5600;yi=0:50:4800;figure(2)z1i=interp2(x,y,z,xi,yi,nearest);surfc(xi,yi,z1i

8、)xlabel(X),ylabel(Y),zlabel(Z)figure(3)z2i=interp2(x,y,z,xi,yi);surfc(xi,yi,z2i)xlabel(X),ylabel(Y),zlabel(Z) figure(4)z3i=interp2(x,y,z,xi,yi,cubic);surfc(xi,yi,z3i)xlabel(X),ylabel(Y),zlabel(Z)figure(5)subplot(1,3,1),contour(xi,yi,z1i,10,r);subplot(1,3,2),contour(xi,yi,z2i,10,r);subplot(1,3,3),con

9、tour(xi,yi,z3i,10,r); 插值函數(shù)插值函數(shù)griddata格式為格式為: cz =griddata（x，y，z，cx，cy，method）用用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算作散點(diǎn)數(shù)據(jù)的插值計(jì)算 要求要求cx取行向量，取行向量，cy取為列向量取為列向量被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值nearest最鄰近插值最鄰近插值linear 雙線性插值雙線性插值cubic 雙三次插值雙三次插值v4- MATLAB提供的插值方法提供的插值方法缺省時(shí)缺省時(shí), , 雙線性插值雙線性插值 例例6. 6. 在某海域測得一些點(diǎn)在某海域測得一些點(diǎn)( (x, ,y) )處的水深處的水深z由下由下

10、表給出，船的吃水深度為表給出，船的吃水深度為5 5英尺，在矩形區(qū)域（英尺，在矩形區(qū)域（7575，200200）（-50-50，150150）里的哪些地方船要避免進(jìn)）里的哪些地方船要避免進(jìn)入入 2.75,20050,150. 在矩形區(qū)域作二維插值三次插值法 .1 輸入插值基點(diǎn)數(shù)據(jù)x=129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 -8 -6 -8 -6 -8 -8 -9 -9

11、 -8 -8 -9 -4 -9;cx=75:0.5:200;cy=-50:0.5:150;cz=griddata(x,y,z,cx,cy,cubic);4.作出水深小于5的海域范圍,即z=5的等高線. 3 作海底曲面圖meshz(cx,cy,cz),rotate3dxlabel(X),ylabel(Y),zlabel(Z)figure(2),contour(cx,cy,cz,-5 -5);gridhold onplot(x,y,+)xlabel(X),ylabel(Y)用用MATLAB解擬合問題解擬合問題1.1.線性最小二乘擬合線性最小二乘擬合2.2.非線性最小二乘擬合非線性最小二乘擬合用用M

12、ATLAB作線性最小二乘擬合作線性最小二乘擬合1. 1. 作多項(xiàng)式作多項(xiàng)式f(x)=a1xm+ +amx+am+1擬合擬合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2. 2. 對(duì)超定方程組對(duì)超定方程組)(11nmyaRnmmn可得最小二乘意義下的解可得最小二乘意義下的解，用，用yRa3.3.多項(xiàng)式在多項(xiàng)式在x處的值處的值y可用以下命令計(jì)算：可用以下命令計(jì)算：y=polyval（a，x）輸出擬合多項(xiàng)式系數(shù)輸出擬合多項(xiàng)式系數(shù)a=a1, ,am , am+1 (數(shù)組數(shù)組) )）輸入同長度輸入同長度的數(shù)組的數(shù)組x，y擬合多項(xiàng)擬合多項(xiàng)式次數(shù)式次數(shù)即要求即要求 出二次多項(xiàng)式出二次

13、多項(xiàng)式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例7. 對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合1）輸入以下命令）輸入以下命令：x=0:0.1:1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;R=(x.2) x ones(11,1);A=Ry211211111 1xxRxx此時(shí)解法解法1 1用解超定方程的方法用解超定方程的方法2）計(jì)算結(jié)果）計(jì)算結(jié)果： = -9.8108 20.1293 -0.03170317.01293.208108.9)

14、(2xxxf1）輸入以下命令：）輸入以下命令： x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; A=polyfit(x,y,2) z=polyval(A,x); plot(x,y,k+,x,z,r) %作出數(shù)據(jù)點(diǎn)和擬合曲線的圖作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形形2）計(jì)算結(jié)果：）計(jì)算結(jié)果： = -9.8108 20.1293 -0.0317解法解法2用多項(xiàng)式擬合的命令用多項(xiàng)式擬合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1. 1. lsqcur

15、vefit已知數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)： xdata=（xdata1,xdata2,xdatan）, ydata=（ydata1,ydata2,ydatan） 用用MATLAB作非線性最小二乘擬合作非線性最小二乘擬合 MATLAB提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin兩個(gè)命令都要先建立兩個(gè)命令都要先建立M M文件文件fun.m，在其中定義函數(shù)，在其中定義函數(shù)f(x)，但兩者定義，但兩者定義f( (x) )的方式是不同的方式是不同的的, ,可參考例題可參考例題.21( ,) niiiF x xdataydata最小 用以求含參量

16、用以求含參量x（向量）的向量值函數(shù)（向量）的向量值函數(shù)F(x,xdata)=（F（x,xdata1）,F（x，xdatan）T中的參變量中的參變量x(x(向量向量),),使得使得 輸入格式為輸入格式為: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit(fun,x0,xdata,ydata,options); (3)x=lsqcurvefit(fun,x0,xdata,ydata,options,grad); (4) x,options=lsqcurvefit(fun,x0,xdata,ydata,); (5) x,o

17、ptions,funval=lsqcurvefit(fun,x0,xdata,ydata,);(6)x,options,funval,Jacob=lsqcurvefit(fun,x0,xdata, ydata,);fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)F(x,xdata) 的的M文件文件, 自變量為自變量為x和和xdata說明:x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知數(shù)據(jù)點(diǎn)已知數(shù)據(jù)點(diǎn)選項(xiàng)見無選項(xiàng)見無約束優(yōu)化約束優(yōu)化 lsqnonlin用以求含參量用以求含參量x（向量）的向量值函數(shù)（向量）的向量值函數(shù) f( (x)

18、)=(=(f1 1( (x),),f2 2( (x),), ,fn( (x)T T 中的參量中的參量x，使得，使得 最小最小 其中其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22221)()()()()(xfxfxfxfxfnT2. lsqnonlin已知數(shù)據(jù)點(diǎn)：已知數(shù)據(jù)點(diǎn)： xdata=（xdata1，xdata2，xdatan） ydata=（ydata1，ydata2，ydatan）輸入格式為：輸入格式為： 1) x=lsqnonlin（fun，x0）； 2)x=lsqnonlin（fun，x0，options）； 3)x= lsqno

19、nlin（fun，x0，optionsgrad）； 4)x，options=lsqnonlin （fun，x0，）； 5)x，options，funval=lsqnonlin（funx0，）；說明：x= lsqnonlin （fun，x0，options）；fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)f(x)的的M文件，文件，自變量為自變量為x迭代初值迭代初值選項(xiàng)見無選項(xiàng)見無約束優(yōu)化約束優(yōu)化100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc100.0221mi

20、n( , , )ejktjjF a b kabc 例例8. 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a，b，k0.0.2( )ektc tab該問題即解最優(yōu)化問題：該問題即解最優(yōu)化問題： 1 1）編寫）編寫M M文件文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中其中 x(1)=a; x(2)=b；x(3)=k;2）輸入命令）輸入命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.

21、39,6.50,6.59;x0=0.2,0.05,0.05;x=lsqcurvefit (curvefun1,x0,tdata,cdata)f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)1010.020.02T(e,e)ktktabab解法解法1 1. 用命令用命令lsqcurvefit3 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為：f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.25424）結(jié)論）結(jié)論：a=0.0063, b=-

22、0.0034, k=0.25421010.020.02T11(e,e)ktktabcabc 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）1）編寫編寫M M文件文件 curvefun2.m function f=curvefun2(x)tdata=100:100:1000;cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=cdata-x(1)-x(2)*exp(-0.02*x(3)*tdata) 2）輸入命令）輸入命令: x0=0.2,0.05,0.05

23、;x=lsqnonlin(curvefun2,x0)f= curvefun2(x)函數(shù)函數(shù)curvefun2的自變量是的自變量是x，cdata和和tdata是已知參數(shù)，故是已知參數(shù)，故應(yīng)將應(yīng)將cdata tdata的值寫在的值寫在curvefun2.m中中3 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為f =1.0e-003 =1.0e-003 * *(0.2322 -0.1243 -0.2495 -0.2413 -0.1668 (0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.0724 0.0241 0.

24、1159 0.2030 0.2792 x =0.0063 -0.0034 0.2542 x =0.0063 -0.0034 0.2542可以看出可以看出,兩個(gè)命令的計(jì)算結(jié)果是相同的兩個(gè)命令的計(jì)算結(jié)果是相同的.4）結(jié)論）結(jié)論：即擬合得即擬合得a=0.0063 0.0063 b=-0.0034 =-0.0034 k=0.2542=0.2542MATLAB解應(yīng)用問題實(shí)例解應(yīng)用問題實(shí)例1. 電阻問題電阻問題 2. 給藥方案問題給藥方案問題*3. 水塔流量估計(jì)問題水塔流量估計(jì)問題電阻問題電阻問題 溫度溫度t(C) 20.5 32.7 51.0 73.0 95.7電阻電阻R( ) 765 826 873

25、942 1032例例. 由數(shù)據(jù)由數(shù)據(jù)擬合擬合R=a1t +b方法方法1.1.用命令用命令 polyfit(x,y,m)t=20.5 32.5 51 73 95.7;r=765 826 873 942 1032;aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,k+,t,y,r)得到得到 a=3.3940, b=702.4918方法方法2.直接用直接用yRa結(jié)果相同結(jié)果相同t=20.5 32.5 51 73 95.7;r=765 826 873 942 1032; R=t ones(5,1);aa=Rr;a=aa(1)b=aa(2)

26、y=polyval(aa,t);plot(t,r,k+,t,y,r)一室模型一室模型：將整個(gè)機(jī)體看作一個(gè)房室，稱：將整個(gè)機(jī)體看作一個(gè)房室，稱中心室中心室，室內(nèi)血藥，室內(nèi)血藥濃度是均勻的快速靜脈注射后，濃度立即上升；然后迅速濃度是均勻的快速靜脈注射后，濃度立即上升；然后迅速下降當(dāng)濃度太低時(shí)，達(dá)不到預(yù)期的治療效果；當(dāng)濃度太高，下降當(dāng)濃度太低時(shí)，達(dá)不到預(yù)期的治療效果；當(dāng)濃度太高，又可能導(dǎo)致藥物中毒或副作用太強(qiáng)臨床上，每種藥物有一又可能導(dǎo)致藥物中毒或副作用太強(qiáng)臨床上，每種藥物有一個(gè)最小有效濃度個(gè)最小有效濃度c1和一個(gè)最大有效濃度和一個(gè)最大有效濃度c2設(shè)計(jì)給藥方案時(shí)，設(shè)計(jì)給藥方案時(shí)，要使血藥濃度要使血藥

27、濃度 保持在保持在c1c2之間本題設(shè)之間本題設(shè)c1=10ug/ml，c2=25ug/ml.擬擬 合合 問問 題題 實(shí)實(shí) 例例 2 2給藥方案給藥方案 一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案. 藥物進(jìn)入機(jī)體后通過血液輸送到全身，在這個(gè)過程中不藥物進(jìn)入機(jī)體后通過血液輸送到全身，在這個(gè)過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為濃度，即單位體積血液中的藥物含量，稱為血藥濃度血藥濃度 在實(shí)驗(yàn)方面在實(shí)驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注對(duì)某人用快速靜脈注射方式

28、一次注入該藥物入該藥物300mg后后,在一定時(shí)刻在一定時(shí)刻t(h)采集血藥采集血藥,測得血測得血藥濃度藥濃度c(ug/ml)如下表如下表: t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01 要設(shè)計(jì)給藥方案要設(shè)計(jì)給藥方案,必須知道給藥后血藥濃度隨必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律從實(shí)驗(yàn)和理論兩方面著手：時(shí)間變化的規(guī)律從實(shí)驗(yàn)和理論兩方面著手：給藥方案給藥方案 1. 在快速靜脈注射的給藥方式下，研究血藥濃度在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物

29、含量）的變化規(guī)律（單位體積血液中的藥物含量）的變化規(guī)律tc2cc1O問問題題2. 給定藥物的最小有效濃度和最大治療濃度，設(shè)計(jì)給定藥物的最小有效濃度和最大治療濃度，設(shè)計(jì)給藥方案：每次注射劑量多大；間隔時(shí)間多長給藥方案：每次注射劑量多大；間隔時(shí)間多長分分析析 理論：用一室模型研理論：用一室模型研究血藥濃度變化規(guī)律究血藥濃度變化規(guī)律 實(shí)驗(yàn)：對(duì)血藥濃實(shí)驗(yàn)：對(duì)血藥濃度數(shù)據(jù)作擬合，符度數(shù)據(jù)作擬合，符合負(fù)指數(shù)變化規(guī)律合負(fù)指數(shù)變化規(guī)律3.3.血液容積血液容積v, , t=0=0注射劑量注射劑量d, , 血藥濃度立即為血藥濃度立即為d/ /v. .2.2.藥物排除速率與血藥濃度成正比，比例系數(shù)藥物排除速率與血藥

30、濃度成正比，比例系數(shù) k(0)(0)模型假設(shè)模型假設(shè)1.1.機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻一室模型一室模型模型建立模型建立3 (0)/cd由假設(shè) 得：d2 -dckct由假設(shè) 得：( )ektdc tv 在此，在此，d=300mg，t及及c（t）在某些點(diǎn)處的值見前表，）在某些點(diǎn)處的值見前表，需經(jīng)擬合求出參數(shù)需經(jīng)擬合求出參數(shù)k、v.用線性最小二乘擬合用線性最小二乘擬合c(t)( )ektdc tv)/ln(,ln21vdakacyktvdc)/ln(ln2/,121aedvakatay計(jì)算結(jié)果：計(jì)算結(jié)果：)(02.15),/1 (2347. 0lvhkd=3

31、00;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)-用用lsqcurvefit2. 主程序如下主程序如下cleartdata=0.25 0.5 1 1.5 2 3 4 6 8;cdata=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01; x0=10,0.5;x=lsqcurvefit(

32、curvefun3,x0,tdata,cdata);f=curvefun3(x,tdata) x1. 1. 用用M M文件文件curvefun3.m定義函數(shù)定義函數(shù)function f=curvefun3(x,tdata)d=300f=(x(1)d)*exp(-x(2)*tdata) % x(1)=v; x(2)=kktevdtc)(給藥方案給藥方案 設(shè)計(jì)設(shè)計(jì)cc2c1Ot 設(shè)每次注射劑量D, 間隔時(shí)間 血藥濃度c(t) 應(yīng)c1 c(t) c2 初次劑量D0 應(yīng)加大,0DD給藥方案記為：給藥方案記為：kecc2112ln1cck2. )( ,1220ccDcD1. 計(jì)算結(jié)果：計(jì)算結(jié)果：9.3,

33、3.225,5.3750DD0375(mg),225(mg),4(h)DD給藥方案：給藥方案：c1=10,c2=25k=0.2347v=15.02故可制定給藥方案：故可制定給藥方案：0375mg,225mg,4hDD即即: 首次注射首次注射375mg， 其余每次注射其余每次注射225mg， 注射的間隔時(shí)間為注射的間隔時(shí)間為4h估計(jì)水塔的流量估計(jì)水塔的流量2. 解題思路解題思路3. 算法設(shè)計(jì)與編程算法設(shè)計(jì)與編程1. 問題問題 某居民區(qū)有一供居民用水的圓柱形水塔，一般可以通過測量其水位來估計(jì)水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時(shí)，水泵自動(dòng)啟動(dòng)向水塔供水，到設(shè)定的最高水位時(shí)停止供

34、水，這段時(shí)間無法測量水塔的水位和水泵的供水量通常水泵每天供水一兩次，每次約兩小時(shí).水塔是一個(gè)高12.2m，直徑17.4m的正圓柱按照設(shè)計(jì)，水塔水位降至約8.2m時(shí)，水泵自動(dòng)啟動(dòng)，水位升到約10.8m時(shí)水泵停止工作表1 是某一天的水位測量記錄，試估計(jì)任何時(shí)刻（包括水泵正供水時(shí)）從水塔流出的水流量，及一天的總用水量 表 1 水位測量記錄 （符號(hào)/表示水泵啟動(dòng)）時(shí)刻(h)水位(cm)0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97968 948 931 913 898 881 869 852 839 822時(shí)刻(h)水位(cm)9.98 10.92 10.

35、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 1018流量估計(jì)的解題思路流量估計(jì)的解題思路擬合水位擬合水位時(shí)間函數(shù)時(shí)間函數(shù)確定流量確定流量時(shí)間函數(shù)時(shí)間函數(shù)估計(jì)一天總用水量估計(jì)一天總用水量 擬合水位擬合水位時(shí)間函數(shù)時(shí)間函數(shù) 從測量記錄看，一天有兩個(gè)供水時(shí)段（以下稱第1供水時(shí)段和第2供水時(shí)段），和3個(gè)水泵不工作時(shí)段（

36、以下稱第1時(shí)段t=0到t=8.97，第2次時(shí)段t=10.95到t=20.84和第3時(shí)段t=23以后）對(duì)第1、2時(shí)段的測量數(shù)據(jù)直接分別作多項(xiàng)式擬合，得到水位函數(shù)為使擬合曲線比較光滑，多項(xiàng)式次數(shù)不要太高，一般在36由于第3時(shí)段只有3個(gè)測量記錄，無法對(duì)這一時(shí)段的水位作出較好的擬合 確定流量確定流量時(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) 一天總用水量的估計(jì)一天總用水量的估計(jì) 總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和，它們都

37、可以由流量對(duì)時(shí)間的積分得到算法設(shè)計(jì)與編程算法設(shè)計(jì)與編程1. 擬合第擬合第1、2時(shí)段的水位，并導(dǎo)出流量時(shí)段的水位，并導(dǎo)出流量2. 擬合供水時(shí)段的流量擬合供水時(shí)段的流量3. 估計(jì)一天總用水量估計(jì)一天總用水量4. 流量及總用水量的檢驗(yàn)流量及總用水量的檢驗(yàn) 1. 擬合第擬合第1時(shí)段的水位，并導(dǎo)出流量時(shí)段的水位，并導(dǎo)出流量 設(shè)t，h為已輸入的時(shí)刻和水位測量記錄（水泵啟動(dòng)的4個(gè)時(shí)刻不輸入），第第1時(shí)段時(shí)段各時(shí)刻的流量可如下得：1） c1=polyfit（t（1：10），），h（1：10），），3）；）； %用3次多項(xiàng)式擬合第1時(shí)段水位，c1輸出3次多項(xiàng)式的系數(shù)2）a1=polyder（c1）；）； % a

38、1輸出多項(xiàng)式（系數(shù)為c1）導(dǎo)數(shù)的系數(shù) 3）tp1=0：0.1：9； x1=-polyval（a1，tp1）；）；% x1輸出多項(xiàng)式（系數(shù)a1）在tp1點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp1時(shí)刻的流量 4）流量函數(shù)為：流量函數(shù)為：1079.227173. 22356. 0)(2tttf 擬合第擬合第2時(shí)段的水位，并導(dǎo)出流量時(shí)段的水位，并導(dǎo)出流量 設(shè)t，h為已輸入的時(shí)刻和水位測量記錄（水泵啟動(dòng)的4個(gè)時(shí)刻不輸入），第第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=po

39、lyder(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í)刻的流量4）流量函數(shù)為：流量函數(shù)為：8313. 17512. 87529. 00186. 0)(23ttttf 2. 擬合供水時(shí)段的流量擬合供水時(shí)段的流量 在第1供水時(shí)段（t=911）之前（即第1時(shí)段）和之后（即第2時(shí)段）各取幾點(diǎn)，其流量已經(jīng)得到，用它們擬合第1供水時(shí)段的流量為使流量函數(shù)在t=9和t=11連續(xù)，我們簡單地只取4個(gè)點(diǎn)，擬合3次多項(xiàng)式（即曲線必過這4個(gè)點(diǎn)），實(shí)

40、現(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 11 12,xx12，3);%擬合3次多項(xiàng)式 tp12=9：0.1：11; x12=polyval(c12,tp12); %x12輸出第1供水時(shí)段 各時(shí)刻的流量擬合的流量函數(shù)為：擬合的流量函數(shù)為：078.3555879.737207. 3)(2tttf 在第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è)水位的兩兩之差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í)刻的流量擬合的流量函數(shù)為：擬合的流量