基于MATLAB實(shí)現(xiàn)的最小二乘擬合-數(shù)值計(jì)算方法

實(shí)驗(yàn)二最小二乘擬合【實(shí)驗(yàn)內(nèi)容】1、理解最小二乘擬合的原理；2、理解適定、欠定、超定三類線性方程組的區(qū)別，會(huì)求其解；3、理解和掌握“左除運(yùn)算”；4、完成最小二乘擬合的程序（過程或函數(shù)）5、會(huì)使用畫圖命令6、與系統(tǒng)函數(shù)進(jìn)行對(duì)比分析和學(xué)習(xí)?！緦?shí)驗(yàn)方法與步驟】一、回答下面的問題1、什么是超定方程，其對(duì)應(yīng)的珍貴方程組是什么，系數(shù)特點(diǎn)是什么？2、最小二乘擬合方法和插值方法的相同點(diǎn)和不同點(diǎn)，各自用在什么情況下？請(qǐng)舉例說明你在其它實(shí)驗(yàn)課中遇到的使用最小二乘擬合方法的場(chǎng)景；3、生成超定方程時(shí)，大家發(fā)現(xiàn)了Matlab的矩陣存儲(chǔ)順序是行優(yōu)先（C語言規(guī)則），還是列優(yōu)先（Fortran語言規(guī)則）？為什么呢？4、時(shí)間測(cè)試下，請(qǐng)對(duì)自己編寫的最小二乘擬合代碼和系統(tǒng)內(nèi)置函數(shù)的執(zhí)行效率比較和分析。在自己的代碼中，將循環(huán)體改寫成等價(jià)的矩陣運(yùn)算，執(zhí)行時(shí)間是否有所縮短？二、完成下列計(jì)算，寫出代碼1、用矩陣逆、左除和rref分別求p74第9題，測(cè)試執(zhí)行時(shí)間2、畫圖命令，畫出習(xí)題5第16、17題的擬合點(diǎn)圖形，進(jìn)行擬合函數(shù)選型3、編寫最小二乘擬合程序，求出習(xí)題5第16、17題的擬合函數(shù)，并繪圖4、執(zhí)行系統(tǒng)函數(shù)polyfit，與自己程序的計(jì)算結(jié)果進(jìn)行比對(duì)和分析。第1頁(yè)共1頁(yè)2015年5月4日星期一制作：王定第2第2頁(yè)共#頁(yè)2015年5月4日星期一制作：王定【實(shí)驗(yàn)結(jié)果】解：一、1:超定方程定義方程個(gè)數(shù)大于未知量個(gè)數(shù)的方程組對(duì)應(yīng)正規(guī)方程組對(duì)于方程組Ra=y,R為nXm矩陣，如果R列滿秩，且n>m。則方程組沒有精確解，此時(shí)稱方程組為超定方程組。系數(shù)特點(diǎn)R為nXm矩陣，R列滿秩，且n>m解：一、2最小二乘法與插值法方法相同點(diǎn)不同點(diǎn)適用情況最小都是從在知道一組實(shí)驗(yàn)數(shù)據(jù)(不準(zhǔn)確值)中數(shù)據(jù)通常是實(shí)驗(yàn)二乘一組數(shù)尋找自變量x和因變量y之間的函數(shù)測(cè)定出來的，出法據(jù)去求關(guān)系y=F(x),用擬合曲線S(x)去逼近于對(duì)精度的要函數(shù)實(shí)驗(yàn)數(shù)據(jù)，來描述自變量x和因變量求，不一定保留y二f(x)y之間的函數(shù)關(guān)系全部數(shù)據(jù)的近似在用最小二乘法求擬合曲線時(shí)，不要表達(dá)式求曲線過給定的一組點(diǎn)。插值在不知道函數(shù)y=f(x)解析式，知道函數(shù)在［a,b］區(qū)間上一系列點(diǎn)的函數(shù)值(準(zhǔn)確值)的前提下，構(gòu)造插值函數(shù)P(x)來代替f(x)，來求非插值節(jié)點(diǎn)的函數(shù)值。構(gòu)造插值函數(shù)時(shí)，要求插值函數(shù)過插值節(jié)丿點(diǎn)a<x0<x1<...<b已知函數(shù)y二f(x)在區(qū)間［a,b］上(n+1)個(gè)互異點(diǎn)x(i=0,1,...,n)上的函數(shù)值y.解：一、3：優(yōu)先問題matlab和Fortran一樣，都是列優(yōu)先的.與之相反的是C++；驗(yàn)證如下：clear結(jié)果：clca=a=ones(3,5)11111sum(a)%母列的兀素相加1111111111ans=33333由于在數(shù)學(xué)發(fā)展過程史中向量默認(rèn)均指列向量，這樣設(shè)計(jì)可以直接服務(wù)已有知識(shí)體系，比如aTa是向量的內(nèi)積，若用MATLAB寫就是a'*a，不存在新的學(xué)習(xí)成本。解：一、4：執(zhí)行效率與時(shí)間分析A=[32;1-1];y=[-l1]';結(jié)果tic,x=inv(A)*y;tocElapsedtimeis0.000000seconds.tic,x=A\y;tocElapsedtimeis0.000000seconds.tic,x=y'/A';tocElapsedtimeis0.000000seconds.說明：x=A\Y比inv(A)*y的計(jì)算效率平均高大約50倍，測(cè)試耗費(fèi)時(shí)間長(zhǎng)短與機(jī)器性能及當(dāng)前并發(fā)程序數(shù)量荷重都有關(guān)系。

解：二、1：矩陣逆、左除和rref分別求p74第9題%juzhenni%zuochuA=[4-2-4;-21710;-4109];A=[4-2-4;-21710;-4109];b=[103-7]';b=[103-7]';x=inv(A)*bx=A\b>>juzhennix=x=2211-1-1>>A=[00.83240.76750.9831;0.8324>>A=[00.83240.76750.9831;0.83240.69300.64000.8190;0.76750.64000.69300.64000.8190;0.76750.64000.59110.59110.7580;0.98310.81900.75800.7580;0.98310.81900.75800.0055];0.0055];b=[8.8997;7.4144;6.8428;4.9171];b=[8.8997;7.4144;6.8428;4.9171];x=A\bx=inv(A)*bx=x=1.0e+002*1.0e+002*0.002735359618640.002735359618642.071386399032832.07138639903280-2.17702240185198-2.177022401852000.036252319058550.03625231905855(3)rrefA=[4-2-4;-21710;-4109];b=[103-7]';>>C=[A,b];U0=rref(C)%用行階梯形解線性方程組，U0的取后列就疋解U0=10020101001-1%rrefA=[00.83240.76750.9831;0.83240.69300.64000.8190;0.76750.64000.59110.7580;0.98310.81900.75800.0055];B=[8.8997;7.4144;6.8428;4.9171];C=[A,B];U0=rref(C);%用行階梯形解線性方程組，U0的最后一列解釋方程的解U0=1.00000000.273501.000000207.1386001.00000-217.70220001.00003.6252解：二、2：畫出習(xí)題5第16、17題的擬合點(diǎn)圖形習(xí)題￡第仃題最小二乘法實(shí)現(xiàn)解：二、3：最小二乘擬合程序，求出習(xí)題5第16、17題的擬合函數(shù)%習(xí)題5第16題，最小二乘法實(shí)現(xiàn)x=[-1.00-0.75-0.50-0.2500.250.500.751.00]';y=[-0.22090.32950.88261.43922.00032.56453.13343.70614.2836]';V=[x.A0,x,x.A2];A=V\y%超定萬程組A=polyfit(x,y,2);plot(x,y,'*')holdonx1=-1.00:0.02:1.00;y1=polyval(A,x1);plot(x1,y1);gridxlabel('x軸'),ylabel('y軸'),title('習(xí)題5第16題最小二乘法實(shí)現(xiàn)')；A=2.00012.25160.0313%習(xí)題5第17題，最小二乘法實(shí)現(xiàn)x=[-3:1:3]';y=[1000012]';V=[x.a0,x,x.a2];A=V\y%超定萬程組A=polyfit(x,y,2);plot(x,y,'*')holdonx1=-3:0.02:3;y1=polyval(A,x1);plot(x1,y1);gridxlabel('x軸'),ylabel('y軸'),title('習(xí)題5第17題最小二乘法實(shí)現(xiàn)')；A=-0.14290.17860.1786解：二、4:執(zhí)行系統(tǒng)函數(shù)polyfit，與自己程序的計(jì)算結(jié)果進(jìn)行比對(duì)和分析。>>x=[0.10.40.50.70.70.9]';>>y=[0.610.920.991.521.472.0