第5講 matlab數(shù)據(jù)擬合

1、1擬擬 合合2.2.擬合的基本原理擬合的基本原理1. 擬合問題引例擬合問題引例2擬擬 合合 問問 題題 引引 例例 1 1溫度溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻電阻R( ) 765 826 873 942 1032已知熱敏電阻數(shù)據(jù)：已知熱敏電阻數(shù)據(jù)：求求60600C時(shí)的電阻時(shí)的電阻R。2040608010070080090010001100 設(shè)設(shè) R=at+ba,b為待定系數(shù)為待定系數(shù)3曲曲 線線 擬擬 合合 問問 題題 的的 提提 法法已知一組（二維）數(shù)據(jù)，即平面上已知一組（二維）數(shù)據(jù)，即平面上 n個(gè)點(diǎn)個(gè)點(diǎn)（xi,yi) i=1,n, 尋求一個(gè)函數(shù)（曲線）尋

2、求一個(gè)函數(shù)（曲線）y=f(x), 使使 f(x) 在某種準(zhǔn)則下與所在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。 +xyy=f(x)(xi,yi)i i 為點(diǎn)為點(diǎn)（xi,yi) 與與曲線曲線 y=f(x) 的距離的距離4擬合與插值的關(guān)系擬合與插值的關(guān)系 函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。的。 實(shí)例：實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到X和 f之間的關(guān)系？x124791 21

3、31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1問題：?jiǎn)栴}：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，這就是數(shù)據(jù)擬合數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點(diǎn)，就是插值問題插值問題；5最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：0246810121416180510152025已知數(shù)據(jù)點(diǎn)spline三次多項(xiàng)式插值0246810121416180510152025已知

4、數(shù)據(jù)點(diǎn)linest三次多項(xiàng)式插值0246810121416180510152025已知數(shù)據(jù)點(diǎn)nearest三次多項(xiàng)式插值6用用MATLAB解擬合問題解擬合問題1 1、線性最小二乘擬合、線性最小二乘擬合2 2、非線性最小二乘擬合、非線性最小二乘擬合7用用MATLAB作線性最小二乘擬合作線性最小二乘擬合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 x處

5、的值處的值y y可用以下命令計(jì)算：可用以下命令計(jì)算： y=polyvaly=polyval（a a，x x）輸出擬合多項(xiàng)式系數(shù)輸出擬合多項(xiàng)式系數(shù)a=a1, am , am+1 (數(shù)組數(shù)組) )）輸入同長(zhǎng)度輸入同長(zhǎng)度的數(shù)組的數(shù)組X，Y擬合多項(xiàng)擬合多項(xiàng)式次數(shù)式次數(shù)8即要求即要求 出二次多項(xiàng)式出二次多項(xiàng)式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例 對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合xi0.10.20.40.50.60.70.80.91yi1.9783.286.167.347.669.589.489.3011.29

6、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=Ry11 11211121xxxxR此時(shí)解法解法1 1用解超定方程的方法用解超定方程的方法2）計(jì)算結(jié)果）計(jì)算結(jié)果： = -9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf101）輸入以下命令：）輸入以下命令： x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.4

7、8 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)(2xxxf111. lsqcurvefit1. lsqcurvefit已知數(shù)據(jù)點(diǎn)數(shù)據(jù)點(diǎn)： xdata=xdata=（xdata1，xdata2，xdataxdatan n），）， ydata=ydat

8、a=（ydataydata1 1，ydataydata2 2，ydataydatan n） 用用MATLAB作非線性最小二乘擬合作非線性最小二乘擬合 Matlab Matlab的提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：的提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。兩個(gè)命令都要先建立。兩個(gè)命令都要先建立M-M-文件文件fun.mfun.m，在其中定義函數(shù)在其中定義函數(shù)f(x)f(x)，但兩者定義，但兩者定義f(x)f(x)的方式是不同的的方式是不同的, ,可參可參考例題考例題.最小 ),(21niiiydataxdata

9、xF lsqcurvefitlsqcurvefit用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù)F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F(xiàn) F（x x，xdataxdatan n）T T中的參變量中的參變量x(x(向量向量),),使得使得 12 輸入格式為輸入格式為: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefit (fun,x0,xdat

10、a,ydata,options,grad); (4) x, options = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, options,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,ydat

11、a,options);迭代初值迭代初值已知數(shù)據(jù)點(diǎn)已知數(shù)據(jù)點(diǎn)選項(xiàng)見無(wú)選項(xiàng)見無(wú)約束優(yōu)化約束優(yōu)化13 lsqnonlin用以求含參量用以求含參量x x（向量）的向量值函數(shù)（向量）的向量值函數(shù) f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的參量中的參量x x，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22221)()()()()(xfxfxfxfxfnT2. lsqnonlin已知數(shù)據(jù)點(diǎn)：已知數(shù)據(jù)點(diǎn)： xdata=xdata=（xdata1，xdata2

12、，xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）14輸入格式為：輸入格式為： 1） x=lsqnonlin（fun，x0）； 2） x= lsqnonlin （fun，x0，options）； 3） x= lsqnonlin （fun，x0，options，grad）； 4） x，options= lsqnonlin （fun，x0，）； 5） x，options，funval= lsqnonlin （fun，x0，）；說明：x= lsqnonlinlsqnonlin （fun，x0，options）

13、；）；fun是一個(gè)事先建立的是一個(gè)事先建立的定義函數(shù)定義函數(shù)f(x)的的M-文件，文件，自變量為自變量為x迭代初值迭代初值選項(xiàng)見無(wú)選項(xiàng)見無(wú)約束優(yōu)化約束優(yōu)化15100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc210102. 0),(minjjktcbeakbaFj 例例2 用下面一組數(shù)據(jù)擬合用下面一組數(shù)據(jù)擬合 中的參數(shù)中的參數(shù)a，b，kktbeatc2 . 0 . 0)(該問題即解最優(yōu)化問題：該問題即解最優(yōu)化問題：16 1 1）編寫）編寫M-M-文件文件 curvefu

14、n1.mcurvefun1.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:1000tdata=100:100:1000cdata=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0.2

15、,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tktktbeabea),(10102. 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit173 3）運(yùn)算結(jié)果為）運(yùn)算結(jié)果為：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.

16、0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542 x = 0.0063 -0.0034 0.25424）結(jié)論）結(jié)論：a=0.0063, b=-0.0034, k=0.254218Tktktcbeacbea),(102. 0102. 0101 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）1）編寫編寫M-M-文件文件 curvefun2.mcurvefun2.m function f=curvef

17、un2(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=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata2）輸入命令）輸入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x)函數(shù)函數(shù)curvefun2的自變量是的自變量是x，cdata和和tdata是已知參數(shù)，故應(yīng)是已知參數(shù)，故應(yīng)將將cdata tdata的值寫在的值寫在curvefun2.m中中193 3）運(yùn)算結(jié)果為）

18、運(yùn)算結(jié)果為 f =1.0e-003 f =1.0e-003 * *(0.2322 -0.1243 -0.2495 -0.2413 (0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.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 b=-0.0

19、034 k=0.25420.0063 b=-0.0034 k=0.254220電阻問題電阻問題 溫度溫度t(0C) 20.5 32.7 51.0 73.0 95.7電阻電阻R( ) 765 826 873 942 1032例例. 由數(shù)據(jù)由數(shù)據(jù)擬合擬合R=a1t+a2方法方法1.1.用命令用命令 polyfit(x,y,m)得到得到 a1=3.3940, a2=702.4918方法方法2.直接用直接用yRa結(jié)果相同。結(jié)果相同。21 在實(shí)驗(yàn)方面在實(shí)驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注對(duì)某人用快速靜脈注射方式一次注入該藥物入該藥物300mg后后,在一定時(shí)刻在一定時(shí)刻t(小時(shí)小時(shí))采集血藥采集血藥,

20、測(cè)測(cè)得血藥濃度得血藥濃度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)和理論兩方面著手：3.3.血液容積血液容積v, t=0v, t=0注射劑量注射劑量d, d, 血藥濃度立即為血藥濃度立即為d/v.d/v.2.2.藥物排除速率與血藥濃度成正比，比例系數(shù)藥物排除速率與血藥濃度成正比，比例系數(shù)

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

22、 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)給藥方案給藥方案 設(shè)計(jì)設(shè)計(jì)cc2c10t 設(shè)每次注射劑量D, 間隔時(shí)間 血藥濃度c(t) 應(yīng)c1 c(t) c2 初次劑量D0 應(yīng)加大,0DD給藥方案記為：給藥方案記為：kecc2112ln1cck2、)( ,1220ccDcD1、計(jì)算結(jié)果：計(jì)算結(jié)果：9 . 3, 3 .225, 5 .3750DD)(4),(225),(3750hmgDmgD給藥方案：給藥方案：c1=10,c2=25k=0.2347v=15.0225故可制定給藥方案：故可制定給藥方案：)(4),(225),(3750hm