MATLAB實現(xiàn)非線性曲線擬合最小二乘法

1、非線性曲線擬合最小二乘法、問題提出設(shè)數(shù)據(jù)(x,y),(i=0丄2,3,4).由表3-1給出，表中第四行為lny=丁，可以iiii看出數(shù)學(xué)模型為yaehx，用最小二乘法確定a及b。i01234xi1.001.251.501.752.00yi5.105.796.537.458.46yi1.6291.7561.8762.0082.135、理論基礎(chǔ)根據(jù)最小二乘擬合的定義:在函數(shù)的最佳平方逼近中f(x)eCa,b，如果f(x)只在一組離散點集x,i=0,l,m,上給定，這就是科學(xué)實驗中經(jīng)常見到的實驗i數(shù)據(jù)(x,y),i=0,1,m的曲線擬合，這里y=f(x)，i=0,1,m,要求一個TOC o 1-5

2、h ziiii函數(shù)y=S*(x)與所給數(shù)據(jù)(x,y),i=0丄,m擬合，若記誤差ii8=S*(x)一y,i=0,1,m,=(6,6,，6)T，設(shè)(x),(x),(x)是Ca,biii01m01n上線性無關(guān)函數(shù)族，在=span(x),(x),(x)中找一函數(shù)S*(x)，使誤差01n平方和2iii=0i=0S(x)ei=08|2=區(qū)62二S*(x)一y2=minS(x)一y2,這里S(x)=a(x)+a(x)HFa(x)(n0是a,b上的權(quán)函數(shù)，它表示不同點(x,f(x)處的數(shù)據(jù)比重不ii同。用最小二乘法求擬合曲線的問題，就是求形如S(x)的一個函數(shù)y=S*(x)，使62=藝(x)S(x)-f(x

3、)2取得最小。它轉(zhuǎn)化為求多元函數(shù)TOC o 1-5 h z12iiii=0I(a,a,a)=藝(x)工a(x)一f(x)201nijjiii=0j=0的極小點(a*,a*,a*)問題。再由求多元函數(shù)極值的必要條件，有01n=2Xw(x)工a(x)一f(x)(x)=0(k=0,l,n)daijjiikiki=0j=0此題中假設(shè)(x)=1，由已知所給數(shù)據(jù)點(x,y)畫出圖形，根據(jù)離散點的ii位置觀察出它們所擬合的曲線圖形應(yīng)類似于指數(shù)函數(shù)的曲線圖形，故設(shè)擬合曲線的函數(shù)為y=aebx。本題編程過程中，令f=y,z1=a,z2=b,令擬合曲線中對應(yīng)x的函數(shù)值與y的差的平方和為J,即J=sum(fy42)

4、;分別求iiJ關(guān)于z1,z2的偏導(dǎo)，簡化后并令其分別為0得一關(guān)于z1,z2的二元非線性方程組,最后利用fsolve命令求得z1,z2的值分別為z1=3.0751z2=0.5052故得到擬合曲線為y=3.0751e0.5052x為證明曲線擬合的正確性，我們將離散點(x,y)與所得的擬合曲線iiy=3.0751e0.5052x畫于同一圖形中，圖形如下：8.5軸y例8的數(shù)據(jù)點(x(i),y(i)和擬合曲線y=f(x)的圖形87.576.565.511.11.21.31.41.51.61.71.81.92x軸5四、結(jié)果分析根據(jù)實驗內(nèi)容求得擬合曲線yaebx中未知數(shù)a,b分別為a=3.0751b=0.5

5、052即擬合曲線為y3.0751e0.5052x。由圖形知擬合成功！參考文獻數(shù)值分析，李慶揚，王能超，易大義，2001,清華大學(xué)出版社(第四版)。數(shù)值方法，關(guān)治，陸金甫，2006，清華大學(xué)出版社。數(shù)值分析與實驗學(xué)習(xí)指導(dǎo)，蔡大用，2001，清華大學(xué)出版社。數(shù)值分析與實驗，薛毅，2005，北京工業(yè)大學(xué)出版社.附錄程序1：symsz1z2x=1.00:0.25:2.00;y=5.10,5.79,6.53,7.45,8.46;f=z1*exp(z2.*x)fy=f-y;J=sum(fy42);Ja=diff(J,z1);Jb=diff(J,z2);Ja1=simple(Ja),Jb1=simple(J

6、b),程序2：functiony1,y2=fun(z)y1=2*z(1)*exp(2*z(2)-51/5*exp(z(2)+2*z(1)*exp(5/2*z(2)-579/50*exp(5/4*z(2)+2*z(1)*exp(3*z(2)-653/50*exp(3/2*z(2)+2*z(1)*exp(7/2*z(2)-149/10*exp(7/4*z(2)+2*exp(4*z(2)*z(1)-423/25*exp(2*z(2);y2=-1/200*z(1)*(-400*z(1)*exp(2*z(2)+2040*exp(z(2)-500*z(1)*exp(5/2*z(2)+2895*exp(5/4*z(2)-600*z(1)*exp(3*z(2)+3918*exp(3/2*z(2)-700*z(1)*exp(7/2*z(2)+5215*exp(7/4*z(2)-800*z(1)*exp(4*z(2)+6768*exp(2*z(2);作圖程序：x=1.00:0.25:2.00;y=5.10,5.79,6.53,7.45,8.46;f=3.0751*exp(0.5052*x);plot(x