數(shù)學(xué)建模中Matlab數(shù)據(jù)擬合應(yīng)用

數(shù)學(xué)建模中Matlab數(shù)據(jù)擬合應(yīng)用第一頁(yè)，共44頁(yè)。例1已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分別用3次和6次多項(xiàng)式曲線擬合這些數(shù)據(jù)點(diǎn).x=0:0.1:1y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)編寫Matlab程序如下:第二頁(yè)，共44頁(yè)。t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)holdonplot(t,s,'r-','linewidth',2)plot(t,s,'b--','linewidth',2)gridx=0:0.1:1y=[-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)第三頁(yè)，共44頁(yè)。例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(shù)據(jù)如表所示:切削時(shí)間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時(shí)間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm第四頁(yè)，共44頁(yè)。解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')第五頁(yè)，共44頁(yè)。解:描出散點(diǎn)圖,在命令窗口輸入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')a=-0.301229.3804holdonplot(t,y1),holdoffa=polyfit(t,y,1)y1=-0.3012*t+29.3804第六頁(yè)，共44頁(yè)。例2用切削機(jī)床進(jìn)行金屬品加工時(shí),為了適當(dāng)?shù)卣{(diào)整機(jī)床,需要測(cè)定刀具的磨損速度.在一定的時(shí)間測(cè)量刀具的厚度,得數(shù)據(jù)如表所示:切削時(shí)間t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度y/cm切削時(shí)間t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度y/cm擬合曲線為:y=-0.3012t+29.3804第七頁(yè)，共44頁(yè)。例3一個(gè)15.4cm×30.48cm的混凝土柱在加壓實(shí)驗(yàn)中的應(yīng)力-應(yīng)變關(guān)系測(cè)試點(diǎn)的數(shù)據(jù)如表所示1.552.472.933.03已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.2.89第八頁(yè)，共44頁(yè)。已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.解選取指數(shù)函數(shù)作擬合時(shí),在擬合前需作變量代換,化為k1,k2的線性函數(shù).于是,令即第九頁(yè)，共44頁(yè)。在命令窗口輸入:x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6]y=[3.103*1.0e+32.465*1.0e+31.953*1.0e+31.517*1.0e+31.219*1.0e+3]z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=[1.552.472.933.032.89]plot(x,w,'*')y1=exp(8.3009)*x.*exp(-494.5209*x)plot(x,w,'*',x,y1,'r-')第十頁(yè)，共44頁(yè)。已知應(yīng)力-應(yīng)變關(guān)系可以用一條指數(shù)曲線來描述,即假設(shè)式中,表示應(yīng)力,單位是N/m2;表示應(yīng)變.擬合曲線為:令則求得于是第十一頁(yè)，共44頁(yè)。在實(shí)際應(yīng)用中常見的擬合曲線有:直線多項(xiàng)式一般n=2,3,不宜過高.雙曲線(一支)指數(shù)曲線第十二頁(yè)，共44頁(yè)。2.非線性曲線擬合:lsqcurvefit.功能:x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)根據(jù)給定的數(shù)據(jù)xdata,ydata(對(duì)應(yīng)點(diǎn)的橫,縱坐標(biāo)),按函數(shù)文件fun給定的函數(shù),以x0為初值作最小二乘擬合,返回函數(shù)fun中的系數(shù)向量x和殘差的平方和resnorm.第十三頁(yè)，共44頁(yè)。例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.首先編寫存儲(chǔ)擬合函數(shù)的函數(shù)文件.functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.^2+x(3)*xdata.^3保存為文件nihehanshu.m第十四頁(yè)，共44頁(yè)。例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.編寫下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;ydata=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];x0=[0,0,0];[x,resnorm]=lsqcurvefit(@nihehanshu,x0,xdata,ydata)第十五頁(yè)，共44頁(yè)。編寫下面的程序調(diào)用擬合函數(shù).xdata=0:0.1:1;ydata=[3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17];x0=[0,0,0];[x,resnorm]=lsqcurvefit(@nihehanshu,x0,xdata,ydata)程序運(yùn)行后顯示x=3.00224.03040.9404resnorm=0.0912第十六頁(yè)，共44頁(yè)。例4已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三個(gè)參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.說明:最小二乘意義上的最佳擬合函數(shù)為f(x)=3ex+4.03x2+0.94x3.此時(shí)的殘差是:0.0912.第十七頁(yè)，共44頁(yè)。f(x)=3ex+4.03x2+0.94x3.擬合函數(shù)為:第十八頁(yè)，共44頁(yè)。練習(xí):1.已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求用三次多項(xiàng)式進(jìn)行擬合的曲線方程.2.已知觀測(cè)數(shù)據(jù)點(diǎn)如表所示xy1.617.72.7491.313.14.1189.43.6110.82.334.50.644.9409.13652.436.9求a,b,c的值,使得曲線f(x)=aex+bsinx+clnx與已知數(shù)據(jù)點(diǎn)在最小二乘意義上充分接近.第十九頁(yè)，共44頁(yè)。插值問題已知n+1個(gè)節(jié)點(diǎn)其中互不相同,不妨設(shè)求任一插值點(diǎn)處的插值節(jié)點(diǎn)可視為由產(chǎn)生,g

表達(dá)式復(fù)雜,甚至無表達(dá)式第二十頁(yè)，共44頁(yè)。1.分段線性插值xjxj-1xj+1x0xn實(shí)用插值方法機(jī)翼下輪廓線2.三次樣條插值細(xì)木條:樣條第二十一頁(yè)，共44頁(yè)。輸入:節(jié)點(diǎn)x0,y0,插值點(diǎn)x(均為數(shù)組,長(zhǎng)度自定義);輸出:插值y(與x同長(zhǎng)度數(shù)組).1.分段線性插值:已有程序

y=interp1(x0,y0,x)y=interp1(x0,y0,x,’linear’)2.三次樣條插值:已有程序

y=interp1(x0,y0,x,’spline’)

或

y=spline(x0,y0,x)用Matlab作插值計(jì)算第二十二頁(yè)，共44頁(yè)。例5對(duì)在[-1,1]上,用n=20的等距分點(diǎn)進(jìn)行分段線性插值,繪制f(x)及插值函數(shù)的圖形.解在命令窗口輸入:x=-1:0.1:1y=1./(1+9*x.^2)xi=-1:0.1:1yi=interp1(x,y,xi)plot(x,y,'r-',xi,yi,'*')第二十三頁(yè)，共44頁(yè)。例6對(duì)在[-5,5]上,用n=11個(gè)等距分點(diǎn)作分段線性插值和三次樣條插值,用m=21個(gè)插值點(diǎn)作圖,比較結(jié)果.解在命令窗口輸入:n=11,m=21x=-5:10/(m-1):5y=1./(1+x.^2)z=0*xx0=-5:10/(n-1):5y0=1./(1+x0.^2)y1=interp1(x0,y0,x)y2=interp1(x0,y0,x,'spline')[x'y'y1'y2']plot(x,z,'r',x,y,'k:',x,y1,'b',x,y2,'g')gtext('Piece.-linear.'),gtext('Spline'),gtext('y=1/(1+x^2)')第二十四頁(yè)，共44頁(yè)。01.00001.00001.00000.50000.80000.75000.82051.00000.50000.50000.50001.50000.30770.35000.29732.00000.20000.20000.20002.50000.13790.15000.14013.00000.10000.10000.10003.50000.07550.07940.07454.00000.05880.05880.05884.50000.04710.04860.04845.00000.03850.03850.0385例6對(duì)在[-5,5]上,用n=11個(gè)等距分點(diǎn)作分段線性插值和三次樣條插值,用m=21個(gè)插值點(diǎn)作圖,比較結(jié)果.xyy1y2第二十五頁(yè)，共44頁(yè)。解在命令窗口輸入:例7在一天24h內(nèi),從零點(diǎn)開始每間隔2h測(cè)得的環(huán)境溫度為12,9,9,10,18,24,28,27,25,20,18,15,13(單位:)推測(cè)在每1s時(shí)的溫度.并描繪溫度曲線.t=0:2:24T=[129910182428272520181513]plot(t,T,'*')ti=0:1/3600:24T1i=interp1(t,T,ti)plot(t,T,'*',ti,T1i,'r-')T2i=interp1(t,T,ti,'spline')plot(t,T,'*',ti,T1i,'r-',ti,T2i,'g-')第二十六頁(yè)，共44頁(yè)。例8在飛機(jī)的機(jī)翼加工時(shí),由于機(jī)翼尺寸很大,通常在圖紙上只能標(biāo)出部分關(guān)鍵點(diǎn)的數(shù)據(jù).某型號(hào)飛機(jī)的機(jī)翼上緣輪廓線的部分?jǐn)?shù)據(jù)如下:x04.749.051938577695114133y05.238.111.9716.1517.116.3414.6312.166.69x152171190y7.033.990第二十七頁(yè)，共44頁(yè)。例8在飛機(jī)的機(jī)翼加工時(shí),由于機(jī)翼尺寸很大,通常在圖紙上只能標(biāo)出部分關(guān)鍵點(diǎn)的數(shù)據(jù).某型號(hào)飛機(jī)的機(jī)翼上緣輪廓線的部分?jǐn)?shù)據(jù)如下:x=[04.749.051938577695114133152171190]y=[05.238.111.9716.1517.116.3414.6312.169.697.033.990]xi=[0:0.001:190]yi=interp1(x,y,xi,'spline')plot(xi,yi)第二十八頁(yè)，共44頁(yè)。例9天文學(xué)家在1914年8月份的7次觀測(cè)中,測(cè)得地球與金星之間距離(單位:m),并取其常用對(duì)數(shù)值與日期的一組歷史數(shù)據(jù)如下所示,試推斷何時(shí)金星與地球的距離(單位:m)的對(duì)數(shù)值為9.9352.日期18202224262830距離對(duì)數(shù)9.96189.95449.94689.93919.93129.92329.9150解由于對(duì)數(shù)值9.9352位于24和26兩天所對(duì)應(yīng)的對(duì)數(shù)值之間,所以對(duì)上述數(shù)據(jù)用三次樣條插值加細(xì)為步長(zhǎng)為1的數(shù)據(jù):第二十九頁(yè)，共44頁(yè)。解由于對(duì)數(shù)值9.9352位于24和26兩天所對(duì)應(yīng)的對(duì)數(shù)值之間,所以對(duì)上述數(shù)據(jù)用三次樣條插值加細(xì)為步長(zhǎng)為1的數(shù)據(jù):x=[18:2:30]y=[9.96189.95449.94689.93919.93129.92329.9150]xi=[18:1:30]yi=interp1(x,y,xi,'spline')A=[xi;yi]A=18.000019.000020.000021.000022.000023.000024.000025.000026.000027.000028.000029.000030.00009.96189.95819.95449.95069.94689.94309.93919.93529.93129.92729.92329.91919.9150第三十頁(yè)，共44頁(yè)。練習(xí):1.設(shè)在區(qū)間[-2,2]上用10等分點(diǎn)作為節(jié)點(diǎn),分別用三種插值方法:(1)計(jì)算并輸出在該區(qū)間的20等分點(diǎn)的函數(shù)值.(2)輸出這個(gè)函數(shù)及兩個(gè)插值函數(shù)的圖形.(3)對(duì)輸出的數(shù)據(jù)和圖形進(jìn)行分析.第三十一頁(yè)，共44頁(yè)。1.設(shè)在區(qū)間[-2,2]上用10等分點(diǎn)作為節(jié)點(diǎn),分別用三種插值方法:(1)計(jì)算并輸出在該區(qū)間的20等分點(diǎn)的函數(shù)值.zi=0.01830.03870.07730.14110.23690.36850.52730.69800.85210.95991.00000.95990.85210.69800.52730.36850.23690.14110.07730.03870.0183第三十二頁(yè)，共44頁(yè)。1.設(shè)在區(qū)間[-2,2]上用10等分點(diǎn)作為節(jié)點(diǎn),分別用兩種插值方法:(2)輸出這個(gè)函數(shù)及兩個(gè)插值函數(shù)的圖形.第三十三頁(yè)，共44頁(yè)。練習(xí):2.已知某型號(hào)飛機(jī)的機(jī)翼斷面下緣輪廓線上的部分?jǐn)?shù)據(jù)如表所示:假設(shè)需要得到x坐標(biāo)每改變0.1時(shí)的y坐標(biāo),分別用兩種插值方法對(duì)機(jī)翼斷面下緣輪廓線上的部分?jǐn)?shù)據(jù)加細(xì),并作出插值函數(shù)的圖形.xy0031.251.772.092.1112.0121.8131.2141.0151.6第三十四頁(yè)，共44頁(yè)。例5給藥方案。

一種新藥用于臨床之前,必須設(shè)計(jì)給藥方案.在快速靜脈注射的給藥方式下,所謂給藥方案是指,每次注射劑量多大,間隔時(shí)間多長(zhǎng).藥物進(jìn)入機(jī)體后隨血液輸送到全身,在這個(gè)過程中不斷地被吸收,分布,代謝,最終排除體外.藥物在血液中的濃度,即單位體積血液中的藥物含量,稱血藥濃度.在最簡(jiǎn)單的一室模型中,將整個(gè)機(jī)體看作一個(gè)房室,稱中心室,室內(nèi)的血藥濃度是均勻的.快速靜脈注射后,濃度立即上升;然后逐漸下降.當(dāng)濃度太低時(shí),達(dá)不到預(yù)期的治療效果;血藥濃度太高,又可能導(dǎo)致藥物中毒或副作用太強(qiáng).臨床上,每種藥物有一個(gè)最小有效濃度c1和一個(gè)最大治療濃度c2.設(shè)計(jì)給藥方案時(shí),要使血藥濃度保持在c1-c2之間.設(shè)本題所研究藥物的最小有效濃度c1=10,最大治療濃度c2=25第三十五頁(yè)，共44頁(yè)。例5給藥方案。

顯然,要設(shè)計(jì)給藥方案,必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律.為此,從實(shí)驗(yàn)和理論兩方面著手.在實(shí)驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時(shí)刻t(小時(shí))采集血樣,測(cè)得血藥濃度c如表:

血藥濃度c(t)的測(cè)試數(shù)據(jù)t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01第三十六頁(yè)，共44頁(yè)。