數(shù)值分析上機(jī)作業(yè)

1、數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告選題：曲線擬合的最小二乘法指導(dǎo)老師：專業(yè)：學(xué)號(hào)：姓名：課題八曲線擬合的最小二乘法問(wèn)題提出從隨機(jī)的數(shù)據(jù)中找出其規(guī)律性，給出其近似表達(dá)式的問(wèn)題，在生產(chǎn)實(shí)踐和科學(xué)實(shí)驗(yàn)量存在，通常利用數(shù)據(jù)的最小二乘法求得擬合曲線。在某冶煉過(guò)程中，根據(jù)統(tǒng)計(jì)數(shù)據(jù)的含碳量與時(shí)間關(guān)系，試求含碳量y與時(shí)問(wèn)t的擬合曲線。t(分)0510152025303540455055_-4y(X10)01.272.162.863.443.874.154.374.514.584.024.64、要求1、用最小二乘法進(jìn)行曲線擬合；2、近似解析表達(dá)式為taita2t2a3t3；3、打印出擬合函數(shù)t,并打印出L與ytj的誤差，j1

2、.2,12；4、另外選取一個(gè)近似表達(dá)式，嘗試擬合效果的比較；5、*繪制出曲線擬合圖*。三、目的和意義1、掌握曲線擬合的最小二乘法；2、最小二乘法亦可用于解超定線代數(shù)方程組；3、探索擬合函數(shù)的選擇與擬合精度間的關(guān)系。四、計(jì)算公式對(duì)于給定的測(cè)量數(shù)據(jù)(xi,fi)(i=1,2,，n),設(shè)函數(shù)分布為my(x)ajj(x)j0特別的，取j(x)為多項(xiàng)式j(luò)(x)xj(j=o,1，m則根據(jù)最小二乘法原理，可以構(gòu)造泛函H(a0,ai,n(fii1majj(Xi)j0Hakk=0,1，m)則可以得到法方程0,求該解方程組,a。a(f,(f,m)m)m)(f,m)則可以得到解a0,a1，am,因此可得到數(shù)據(jù)的最小

3、二乘解f(x)majj(x)j0曲線擬合：實(shí)際工作中，變量間未必都有線性關(guān)系，如服藥后血藥濃度與時(shí)間的關(guān)系；疾病療效與療程長(zhǎng)短的關(guān)系；毒物劑量與致死率的關(guān)系等常呈曲線關(guān)系。曲線擬合是指選擇適當(dāng)?shù)那€類型來(lái)擬合觀測(cè)數(shù)據(jù)，并用擬合的曲線方程分析兩變量間的關(guān)系。五、結(jié)構(gòu)程序設(shè)計(jì)在程序結(jié)構(gòu)方面主要是按照順序結(jié)構(gòu)進(jìn)行設(shè)計(jì)，在進(jìn)行曲線的擬合時(shí)，為了進(jìn)行比較，在程序設(shè)計(jì)中，直接調(diào)用了最小二乘法的擬合函數(shù)polyfit,并且依次調(diào)用了plot>figure、holdon函數(shù)進(jìn)行圖象的繪制，最后調(diào)用了一個(gè)絕對(duì)信函數(shù)abs用于計(jì)算擬合函數(shù)與原有數(shù)據(jù)的誤差，進(jìn)行擬合效果的比較,，22,35.1 用一元三次多項(xiàng)

4、式t知a2ta3t進(jìn)行擬合計(jì)算解析表達(dá)式系數(shù)：a1,a2,a3t=0510152025303540455055;y=01.272.162.863.443.874.154.374.514.584.024.64;>>n=length(xi);f=0.34364.*10A(-4)*x.A3-5.2156.*10A(-3)*x.A2+0.26340.*x+0.017839;x=0:0.01:55;F=0.34364.*10.a(-4)*x.a3-5.2156.*10.a(-3)*x.A2+0.26340.*x+0.017839;fy=abs(f-y);fy2=fy.A2;Ew=max(fy

5、),E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'t*'),holdon,plot(t,F,'b-'),holdoff4332所得函數(shù)為0.3436410t5.215610t0.2634t0.013839運(yùn)行后屏幕顯示數(shù)據(jù)(xi,yi)與擬合函數(shù)f的最大誤差Ew平均誤差E1和均方根誤差E2及其數(shù)據(jù)點(diǎn)(xi，yi)和擬合曲線y=f(x)的圖形如圖5.1.Ew=0.4243E1=0.0911E2=0.14675.2用一兀四次多項(xiàng)式tata2t2a3t3a4進(jìn)行擬合:計(jì)算多項(xiàng)式系數(shù)：a1,a2,a3,a4xi=051015202

6、5303540455055;y=01.272.162.863.443.874.154.374.514.584.024.64;n=length(xi);x=0:0.01:55;f=0.6026.*10.A(-6)*x.A4-0.31918.*10.A(-4)*x.A3-0.0029323.*xA2+0.23807.*x+0.060449;x=0:0.01:55;F=0.6026.*10.a(-6)*x.a4-0.31918.*10.a(-4)*x.a3-0.0029323.*x.A2+0.23807.*x+0.060449;fy=abs(f-y);fy2=fy.A2;Ew=max(fy),E1=

7、sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'),holdon,plot(x,F,'b-'),holdoff所得函數(shù)為644332(t)0.602610t0.3191810t2.932310t0.23807t0.060449運(yùn)行后屏幕顯示數(shù)據(jù)(xi，yi)與擬合函數(shù)f的最大誤差Ew,平均誤差E1和均方根誤差E2及其數(shù)據(jù)點(diǎn)(xi，yi)和擬合曲線y=f(x)的圖形如圖5.2。Ew=0.3897E1=0.1034、E2=0.1429圖5.2一元四次多項(xiàng)式擬合曲線誤差圖，25.3用一元二次多項(xiàng)式ta1a2ta3t進(jìn)行擬合：計(jì)算多

8、項(xiàng)式系數(shù)：ai,a2,a3輸入程序：>>symsaia2a3x=0510152025303540455055;fi=a1.*x.A2+a2.*x+a3運(yùn)行后屏幕顯示關(guān)于ai,a2和a3的線性方程組：fi=a3,25*a1+5*a2+a3,100*a1+10*a2+a3,225*a1+15*a2+a3,400*a1+20*a2+a3,625*a1+25*a2+a3,900*a1+30*a2+a3,1225*a1+35*a2+a3,1600*a1+40*a2+a3,2025*a1+45*a2+a3,2500*a1+50*a2+a3,3025*a1+55*a2+a3編寫構(gòu)造誤差平方和的M

9、ATLAB?序：y=01.272.162.863.443.874.154.374.514.584.024.64;fi=a3,25*a1+5*a2+a3,100*a1+10*a2+a3,225*a1+15*a2+a3,400*a1+20*a2+a3,625*a1+25*a2+a3,900*a1+30*a2+a3,1225*a1+35*a2+a3,1600*a1+40*a2+a3,2025*a1+45*a2+a3,2500*a1+50*a2+a3,3025*a1+55*a2+a3;fy=fi-y;fy2=fy.A2;J=sum(fy.A2)運(yùn)行后屏幕顯示誤差平方和如下：J=(100*a1+10*a

10、2+a3-54/25卜2+(25*a1+5*a2+a3-127/100)A2+(225*a1+15*a2+a3-143/502+(400*a1+20*a2+a3-86/25)A2+(900*a1+30*a2+a3-83/202+(625*a1+25*a2+a3-387/100)A2+(1225*a1+35*a2+a3-437/1002+(1600*a1+40*a2+a3-451/100)A2+(2025*a1+45*a2+a3-229/502+(2500*a1+50*a2+a3-201/502+(3025*a1+55*a2+a3-116/252+a3A20為求a1,a2,a3,a4使J達(dá)到最小

11、，只需利用極值的必要條件為(kl2,3,4),得到關(guān)于a1，a2，a3，a4的線性方程組，這可以由下面的MATLAB?序完成，即輸入程序：>>symsa1a2a3J=(100*a1+10*a2+a3-54/252+(25*a1+5*a2+a3-127/100)A2+(225*a1+15*a2+a3-143/502+(400*a1+20*a2+a3-86/25)A2+(900*a1+30*a2+a3-83/202+(625*a1+25*a2+a3-387/100)A2+(1225*a1+35*a2+a3-437/1002+(1600*a1+40*a2+a3-451/100)A2+(2

12、025*a1+45*a2+a3-229/502+(2500*a1+50*a2+a3-201/502+(3025*a1+55*a2+a3-116/252+a3A2;Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3),運(yùn)行后屏幕顯示J分別對(duì)a1,a2,a3的偏導(dǎo)數(shù)如下：Ja11=49967500*a1+1089000*a2+25300*a3-217403/2Ja21=1089000*a1+25300*a2+660*a3-27131/10Ja31=25300*a

13、1+660*a2+24*a3-3987/50解線性方程組Ja11=0,Ja21=0,Ja31=0輸入下列程序：>>A=49967500,1089000,25300;1089000,25300,660;25300,660,24;B=217403/2,27131/10,3987/50;C=B/A,F=poly2sym(C)運(yùn)行后屏幕顯示擬合函數(shù)f及其系數(shù)C如下：C=-0.00240.20370.2305F=(64133*x)/-(52799*xA2)/52+32545/故所求的擬合曲線為：f(x2)0.0023805x20.20369x0.23047編寫下面的MATLAB?序估計(jì)其誤差

14、，并作出擬合曲線和數(shù)據(jù)的圖形。輸入程序：>>xi=0510152025303540455055;y=01.272.162.863.443.874.154.374.514.584.024.64;n=length(xi);f=-0.0023805.*x.A2+0.20369.*x+0.23047;x=0:0.01:55;F=-0.0023805.*x.A2+0.20369.*x+0.23047;fy=abs(f-y);fy2=fy.A2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt(sum(fy2)/n)plot(xi,y,'r*'),holdon,p

15、lot(x,F,'b-'),holdofflegend('數(shù)據(jù)點(diǎn)(xi,yi)',"擬合曲線y=f(x)'),xlabel(,x,),ylabel(,y,),title(,數(shù)據(jù)點(diǎn)(xi,yi)和擬合曲線y=f(x)的圖形,)運(yùn)行后屏幕顯示數(shù)據(jù)(xi，yi)與擬合函數(shù)f的最大誤差Ew,平均誤差E1和均方根誤差E2及其數(shù)據(jù)點(diǎn)(xi，yi)和擬合曲線y=f(x)的圖形如圖5.3所示：Ew=0.4437E1=0.1426E2=0.19830圖5.3一元二次多項(xiàng)式擬合曲線誤差圖六、結(jié)果討論和分析:由以上結(jié)果可知，擬合方程的選取至關(guān)重要，它決定了最大誤差、平均誤差以及均方根誤差的大小，即擬合曲線的接近程度。本次實(shí)驗(yàn)，最初所選取的,2,3擬合解析方程ta1ta2ta3t獲得較好的擬合，選用解析方程為ta1ta2t2a3t3a4t4的曲線擬合時(shí)，精確度有所下降。由此，擬合函數(shù)的選擇和擬合精度致密相關(guān)，最小二乘法如果想將曲線擬合的比較完美，必須應(yīng)用適當(dāng)?shù)哪M曲線，如果模擬曲線選擇不夠適當(dāng)，那么用最小二乘法計(jì)算完后，會(huì)發(fā)現(xiàn)擬合曲線誤差比較大，均方誤差也比較大，而如果擬合曲線選擇適當(dāng)，那么效