Matlab數(shù)據擬合教學文稿課件

Matlab數(shù)據擬合.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)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)grid例2用切削機床進行金屬品加工時,為了適當?shù)卣{整機床,需要測定刀具的磨損速度.在一定的時間測量刀具的厚度,得數(shù)據如表所示:切削時間t/h012345678刀具厚度y/cm30.029.128.428.128.027.727.527.227.0切削時間t/h910111213141516刀具厚度y/cm26.826.526.326.125.725.324.824.0解:描出散點圖,在命令窗口輸入: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,'*')解:描出散點圖,在命令窗口輸入: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=polyfit(t,y,1)a=-0.301229.3804holdony1=-0.3012*t+29.3804plot(t,y1),holdoff例2用切削機床進行金屬品加工時,為了適當?shù)卣{整機床,需要測定刀具的磨損速度.在一定的時間測量刀具的厚度,得數(shù)據如表所示:切削時間t/h012345678刀具厚度y/cm30.029.128.428.128.027.727.527.227.0切削時間t/h910111213141516刀具厚度y/cm26.826.526.326.125.725.324.824.0擬合曲線為:y=-0.3012t+29.3804例3一個15.4cm×30.48cm的混凝土柱在加壓實驗中的應力-應變關系測試點的數(shù)據如表所示2?/?N/m?1.55500?10?62.471000?10?62.931500?10?63.032.89?2?6?62000?102375?10??//?N/m?3.103?102.465?103331.219?10.517?101.953?10133已知應力-應變關系可以用一條指數(shù)曲線來描述,即假設??k?e1?k?2?表示應力,單位是N/m2;?表示應變.式中,已知應力-應變關系可以用一條指數(shù)曲線來描述,即假設??k?e1?k?2?表示應力,單位是N/m2;?表示應變.式中,解選取指數(shù)函數(shù)作擬合時,在擬合前需作變量代換,化為k1,k2的線性函數(shù).?于是,ln?lnk1?k2??令即?z?ln,a?k,alnk0?21?1?z?a??a01在命令窗口輸入: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-')已知應力-應變關系可以用一條指數(shù)曲線來描述,即假設??k?e1?k?2?表示應力,單位是N/m2;?表示應變.式中,令?z?ln,a?k,alnk,0?21?1?3則z?a??a01求得于是a??k?-494.5209,ak?ln?8.3009,0211k?4.0275?10,k?494.520912擬合曲線為:??4.0275?10?e3

-494.5209?在實際應用中常見的擬合曲線有:直線y?a0x?a1nn?1多項式y(tǒng)a一般n=2,3,不宜過高.?x?ax??La01n雙曲線(一支)a0y??a1x指數(shù)曲線y?aebx2.非線性曲線擬合:lsqcurvefit.x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)功能:根據給定的數(shù)據xdata,ydata(對應點的橫,縱坐標),按函數(shù)文件fun給定的函數(shù),以x0為初值作最小二乘擬合,返回函數(shù)fun中的系數(shù)向量x和殘差的平方和resnorm.例4已知觀測數(shù)據點如表所示x00.10.20.30.40.50.60.70.80.91y3.13.273.814.55.1867.058.569.6911.2513.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據點在最小二乘意義上充分接近.首先編寫存儲擬合函數(shù)的函數(shù)文件.functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.^2+x(3)*xdata.^3保存為文件nihehanshu.m例4已知觀測數(shù)據點如表所示x00.10.20.30.40.50.60.70.80.91y3.13.273.814.55.1867.058.569.6911.2513.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據點在最小二乘意義上充分接近.編寫下面的程序調用擬合函數(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)編寫下面的程序調用擬合函數(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)程序運行后顯示x=3.00224.03040.9404resnorm=0.0912例4已知觀測數(shù)據點如表所示x00.10.20.30.40.50.60.70.80.91y3.13.273.814.55.1867.058.569.6911.2513.17求三個參數(shù)a,b,c的值,使得曲線f(x)=aex+bx2+cx3與已知數(shù)據點在最小二乘意義上充分接近.說明:最小二乘意義上的最佳擬合函數(shù)為f(x)=3ex+4.03x2+0.94x3.此時的殘差是:0.0912.擬合函數(shù)為:f(x)=3ex+4.03x2+0.94x3.練習:1.已知觀測數(shù)據點如表所示x00.10.20.30.40.50.60.70.80.91y3.13.273.814.55.1867.058.569.6911.2513.17求用三次多項式進行擬合的曲線方程.2.已知觀測數(shù)據點如表所示x1.62.71.34.13.62.3y17.74913.1189.4110.834.50.644.9409.13652.436.9求a,b,c的值,使得曲線f(x)=aex+bsinx+clnx與已知數(shù)據點在最小二乘意義上充分接近.插值問題xx,y)(j?0,1,?n,其中已知n+1個節(jié)點(jj?????x),互不相同,不妨設ax0x1Lnb求任一插值點jx(?xj)??*處的插值y.*yy1y0*(xy,j)(j?0,1,Ln)j????節(jié)點可視為由y?g(x)產生,g表達式復雜,甚至無表達式x0x1x*xn實用插值方法1.分段線性插值??????xj-1xjxj+12.三次樣條插值x0xn機翼下輪廓線細木條:樣條用Matlab作插值計算1.分段線性插值:已有程序y=interp1(x0,y0,x)y=interp1(x0,y0,x,'linear')輸入:節(jié)點x0,y0,插值點x(均為數(shù)組,長度自定義);輸出:插值y(與x同長度數(shù)組).2.三次樣條插值:已有程序y=interp1(x0,y0,x,'spline')或y=spline(x0,y0,x)例5對1f?x??21?9x在[-1,1]上,用n=20的等距分點進行分段線性插值,繪制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,'*')1例6對y?在[-5,5]上,用n=11個等距分點作分段線21?x性插值和三次樣條插值,用m=21個插值點作圖,比較結果.解在命令窗口輸入: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)')1例6對y?在[-5,5]上,用n=11個等距分點作分段線21?x性插值和三次樣條插值,用m=21個插值點作圖,比較結果.xyy1y201.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例7在一天24h內,從零點開始每間隔2h測得的環(huán)境溫度為12,9,9,10,18,24,28,27,25,20,18,15,13(單位:)?C推測在每1s時的溫度.并描繪溫度曲線.解在命令窗口輸入: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-')例8在飛機的機翼加工時,由于機翼尺寸很大,通常在圖紙上只能標出部分關鍵點的數(shù)據.某型號飛機的機翼上緣輪廓線的部分數(shù)據如下:x04.749.051938577695114133y05.238.111.9716.1517.116.3414.6312.166.69x152171190y7.033.990例8在飛機的機翼加工時,由于機翼尺寸很大,通常在圖紙上只能標出部分關鍵點的數(shù)據.某型號飛機的機翼上緣輪廓線的部分數(shù)據如下: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)例9天文學家在1914年8月份的7次觀測中,測得地球與金星之間距離(單位:m),并取其常用對數(shù)值與日期的一組歷史數(shù)據如下所示,試推斷何時金星與地球的距離(單位:m)的對數(shù)值為9.9352.日期18202224262830距離9.96189.95449.94689.93919.93129.92329.9150對數(shù)解由于對數(shù)值9.9352位于24和26兩天所對應的對數(shù)值之間,所以對上述數(shù)據用三次樣條插值加細為步長為1的數(shù)據:解由于對數(shù)值9.9352位于24和26兩天所對應的對數(shù)值之間,所以對上述數(shù)據用三次樣條插值加細為步長為1的數(shù)據: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練習:1.設f?x??e,在區(qū)間[-2,2]上用10等分點作為節(jié)點,分別用三種插值方法:?x2(1)計算并輸出在該區(qū)間的20等分點的函數(shù)值.(2)輸出這個函數(shù)及兩個插值函數(shù)的圖形.(3)對輸出的數(shù)據和圖形進行分析.1.設f?x??e,在區(qū)間[-2,2]上用10等分點作為節(jié)點,分別用三種插值方法:?x2(1)計算并輸出在該區(qū)間的20等分點的函數(shù)值.zi=0.01830.03870.07730.14110.23690.36850.52730.69800.85210.95991.00000.95990.85210.69800.52730.36850.23690.14110.07730.03870.01831.設f?x??e,在區(qū)間[-2,2]上用10等分點作為節(jié)點,分別用兩種插值方法:?x2(2)輸出這個函數(shù)及兩個插值函數(shù)的圖形.練習:2.已知某型號飛機的機翼斷面下緣輪廓線上的部分數(shù)據如表所示:xy0031.251.772.092.1112.0121.8131.2141.0151.6假設需要得到x坐標每改變0.1時的y坐標,分別用兩種插值方法對機翼斷面下緣輪廓線上的部分數(shù)據加細,并作出插值函數(shù)的圖形.例5給藥方案一種新藥用于臨床之前,必須設計給藥方案.在快速靜脈注射的給藥方式下,所謂給藥方案是指,每次注射劑量多大,間隔時間多長.藥物進入機體后隨血液輸送到全身,在這個過程中不斷地被吸收,分布,代謝,最終排除體外.藥物在血液中的濃度,即單位體積血液中的藥物含量,稱血藥濃度.在最簡單的一室模型中,將整個機體看作一個房室,稱中心室,室內的血藥濃度是均勻的.快速靜脈注射后,濃度立即上升;然后逐漸下降.當濃度太低時,達不到預期的治療效果;血藥濃度太高,又可能導致藥物中毒或副作用太強.臨床上,每種藥物有一個最小有效濃度c1和一個最大治療濃度c2.設計給藥方案時,要使血藥濃度保持在c1-c2之間.設本題所研究藥物的最小有效濃度c1=10,最大治療濃度c2=25??g/ml?例5給藥方案顯然,要設計給藥方案,必須知道給藥后血藥濃度隨時間變化的規(guī)律.為此,從實驗和理論兩方面著手.在實驗方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血樣,測得血藥濃度c如表:血藥濃度c(t)的測試數(shù)據t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01例5給藥方案t:logc近似直線關系,即c(t)有按負指數(shù)規(guī)律減少的趨勢.例5給藥方案1.確定血藥濃度的變化規(guī)律假設:a)藥物向體外排除的速率與中心室的血藥濃度成正比,比例系數(shù)為k(>0),稱排除速率.b)中心室血液容積為常數(shù)V,t=0瞬時注入藥物的d劑量為d,血藥濃度立即為.dc??kc由假設a),中心室的血藥濃度c(t)應滿足微分方程dtd由假設b),方程的初始條件為:c?0??.Vd?kt求解得:c?t??e.即血藥濃度c(t)按指數(shù)規(guī)律下降.VV2.給藥方案設計簡單實用的給藥方案是:每隔一定時間?,重復注入固定劑量D,使血藥濃度c(t)呈周期性變化,并保持在c1-c2之間.yc2c10?x2.給藥方案設計簡單實用的給藥方案是:每隔一定時間?,重復注入固定劑量D,使血藥濃度c(t)呈周期性變化,并保持在c1-c2之間.為此,初次劑量需加大到D0.d?kt由式c?t??e得到:V1c2DV?c,D?Vc?c,??ln??0221kc1,??,顯然,當c1,c2給定后,要確定給藥方案?D0,D必須知道參數(shù)V和k.2.由實驗數(shù)據作曲線擬合以確定參數(shù)d?kt為了用線性最小二乘法擬合c?t??e的系數(shù)V和k,V先取對數(shù)得記dy?ln,ca?k,aln1?2?Vdlnc?ln?kt