南京郵電大學數(shù)值計算實踐報告

1、數(shù)值計算實踐I、方程求根一、實驗?zāi)康氖煜ず驼莆誑ewton法,割線法，拋物線法的方法思路，并能夠在 matlab上編 程實現(xiàn)二、問題描述(1) .給定一個三次方程，分別用Newton法,割線法，拋物線法求解.方程的構(gòu)造方法：(a)根:方程的根為學號的后三位乘以倒數(shù)第二位加1再除以1000.假設(shè)你的學號為 B06060141，則根為141*(4+1)/1000=0.564(b)方程:以你的學號的后三位數(shù)分別作為方程的三次項，二次項，一次項的系數(shù)，根 據(jù)所給的根以及三個系數(shù)確定常數(shù)項.例如：你的學號是B06060141，則你的方程是x3+4x2+x+a0=0的形式.方程的根為0.564,因止匕有0

2、.5643+4*0.5642+0.564+a0=0，于是 a0=-2.015790144你的方程為 x3+4x2+x-2.015790144=0.(2)假設(shè)方程是sinx+4x2+x+a0=0的形式(三個系數(shù)分別是學號中的數(shù)字)，重新 解決類似的問題(3)構(gòu)造一個五次方程完成上面的工作.四次方程的構(gòu)造：將三次多項式再乘以(x-p*)2得到對應(yīng)的五次多項式(p*為已經(jīng)確 定的方程的根，顯然，得到的五次方程有重根).(4)將(2)中的方程同樣乘以(x-p*)得到一個新的方程來求解注:(1)Newton法取0.5為初值，割線法以0,1為初值，拋物線法以0,0.5,1為初值, (2)計算精度盡量地取高

3、.終止準則：根據(jù)| pn - pnFW來終止可供研究的問題:（一）8的取值不同對收斂速度有多大的影響（二）將注（1）中的初值該為其它的初值，對收斂性以及收斂速度有無影響（二）能否求出方程的所有一的根（4）實驗報告的撰寫實驗報告包含的內(nèi)容：（一）實驗?zāi)康模ǘ﹩栴}描述（三）算法介紹（包 括基本原理）（四）程序（五）計算結(jié)果（六）結(jié)果分析（七）心得體會三、算法介紹在本問題中，我們用到了 newton法，割線法，拋物線法。1 .Newton法迭代格式為：f （xk）一“一國當初值與真解足夠靠近,newton迭代法收斂，對于單根,newton收斂速度 很快，對于重根，收斂較慢。2 .割線法：為了回避導

4、數(shù)值的計算，使用上的差商代替，得到割線法迭代 公式：f （Xk）（Xk -Xk j）Xk 1 = xk -f（Xk） -f （XkJ割線法的收斂階雖然低于newton法，但迭代以此只需計算一次函數(shù)值，不 需計算其導數(shù)，所以效率高，實際問題中經(jīng)常應(yīng)用。3 .拋物線法：可以通過三點做一條拋物線，產(chǎn)生迭代序列的方法稱為拋物 線法。其迭代公式為：p（X）= fXk fXk,Xk_i（X-Xk）fXk,Xk_i,Xk/（X-Xk）（X-Xk_i）其中fXk,Xk, fXk, Xk山Xkj是一階均差和二階均差。收斂速度比割線法更接近于newton法。對于本問題的解決就以上述理論為依據(jù)。終止準則為：| Xn

5、 -Xn|< 8本題中所有精度取1e-8。四、程序計算結(jié)果問題一根據(jù)所給的要求，可知待求的方程為：X3 +2x2 +x-0.3356 =0牛頓法建立newton_1.m的源程序，源程序代碼為：function y=newton_1（a,n,x0,nn,eps1）x（1）=x0;b=1;i=1;while(abs(b)>eps1*x(i)i=i+1;x(i)=x(i-1)-n_f(a,n,x(i-1)/n_df(a,n,x(i-1);b=x(i)-x(i-1);if(i>nn)errorreturn;endendy=x(i);建立n_f.m的源程序來求待求根的實數(shù)代數(shù)方程的函數(shù)

6、，源程序代碼為: function y=n_f(a,n,x)% 待求根的實數(shù)代數(shù)方程的函數(shù)y=0.0;for i=1:(n+1)y=y+a(i)*xA(n+1-i);end建立n_df.m的源程序來方程的一階導數(shù)的函數(shù)，源程序代碼為：function y=n_df(a,n,x)% 方程的一階導數(shù)的函數(shù)y=0.0;for i=1:ny=y+a(i)*(n+1-i)*xA(n-i);end在matlata件中執(zhí)行下列語句并得到的最終結(jié)果截圖：» a=l, 2, 1, -0.3366:» *3;» x0=0. 5;» nn=10C0：>> epsl

7、=le-8；>> y=nevt on_ 1 fa, nr xO, eps 1)0. 2240割線法建立gexian.m的源程序，源程序代碼為function x=gexian(f,x0,x1,e)if nargin<4,e=1e-4;endy=x0;x=x1;i=0;while abs(x-y)>ei=i+1;z=x-(feval(f,x)*(x-y)/(feval(f,x)-feval(f,y);y=x;x=z;endi在matlata件中執(zhí)行下列語句并得到的最終結(jié)果截圖:>> fun= inline ( x 3+2+x 2+x-0. 3356：>&

8、gt; gexian (fun, 0, 1, le-S)ans =0. 2240拋物線建立paowuxian.m的源程序，源程序代碼為：function x=paowuxian(f,x0,x1,x2,e)if nargin<4,e=1e-4;endx=x2;y=x1;z=x0;i=0;while abs(x-y)>ei=i+1;h1=y-z;h2=x-y;c1=(feval(f,y)-feval(f,z)/h1;c2=(feval(f,x)-feval(f,y)/h2;d=(c1-c2)/(h1+h2);w=c2+h2*d;xi=x-(2*feval(f,x)/(w+(w/abs(

9、w)*sqrt(wA2-4*feval(f,x)*d);z=y;y=x;x=xi;end在matlata件中執(zhí)行下列語句并得到的最終結(jié)果截圖:>> funFinlineCx'3+2*x 2+x-O,3356h);>> paowuxiandurij 必 0, 53 L le-8)CL 2240研究一：只改變初值由上述結(jié)果可知，方程的解在0.2附近，所以將牛頓法為0.2;割線法的初值設(shè)為0, 0.4;拋物線法的初值設(shè)為0, 0.2, 0.4;牛頓法根據(jù)問題1中牛頓法的程序，在matlab軟件中執(zhí)行下列語句并得到的最終結(jié)果截 圖：» H 2, 1,-0,33

10、56;» n=3;» x0=0.2;» epsl=le-8;» nn=lDOO;» y=nevrton_1 區(qū) 延 nn, epsl)0. 2240割線法根據(jù)問題1中割線程序，在matlab軟件中執(zhí)行下列語句并得到的最終結(jié)果截圖:>> finline ( x" 3+2*k * 2+x-0. 3356')>> gestianlfun, Cl, . % 1e-8)i -ans =0. 2240拋物線法根據(jù)問題1中拋物線法程序，在 matlab軟件中執(zhí)行下列語句并得到的最終結(jié)果 截圖：» fun=i

11、nline( x 2+x-O. 3356)：>> paowuxian(fun, 0* 2t 0* 4, 1 s-8)ans 二Q. 2240研究二只改變精度將精度由1e-8改為1e-50和1e-100觀察迭代次數(shù)有何變化牛頓法：根據(jù)問題1中的牛頓法的程序，在 matlab軟件中執(zhí)行下列語句并得到的最終結(jié)果截圖：精度為1e-50時» 卒1,2. L-0, 3356;» n=3;» xO=O. E:» nn=1000;>> epsl=1e_50;» y=nevrt on_ 1 (n, x0, nrij ep s 1)I 22

12、40精度為1e-100時» epsl=le-lC0;» y=nevrton_l (a, n3 kOj nn, epsl)0. 2240割線法根據(jù)問題1中的割線法的程序，在 matlab軟件中執(zhí)行下列語句并得到的最終結(jié)果截圖：精度為1e-50時>> fun= inline(工"3+2株"2+x-0. 3356 );>> gexian (fun, 0, 1, le-BO)g ans -精度為1e-100時>> fun= inline工-3+2"2+x-04 3356?);>> gexianffurij

13、 0? lf ls_100)ans =0.2240拋物線法根據(jù)問題1中的拋物線法的程序，在 matlab軟件中執(zhí)行下列語句并得到的最終 結(jié)果截圖：精度為1e-50時>> fun=inline ( x* 3+2*x"2+x-0.3356');» paowuxian(fun, 0, 0» 5j L le-505i =10ars =0. 2240精度為1e-100時>> paowuxian(funj 0 0* 53 1 j le-100)10ans =0. 2240研究結(jié)論在只改變初值時，當初值定得越靠近初值，迭代次數(shù)就越少 在只改變精度

14、時，當精度越來越大時，迭代次數(shù)并幾乎不變。綜上所述，初值對迭代次數(shù)的影響比較大，精度對迭代次數(shù)影響不大問題二問題描述根據(jù)所給的要求，可知待求的方程為：sin x + 2x2 +x-0.3356 = 0問題分析仍然利用（1）中方法求解這一問題，并利用圖解法找到初值，通過觀察圖像， 將newton法初值設(shè)為：0.1,割線法初值設(shè)為：0,0.2。拋物線法初值設(shè)為：0,0.1,02 圖像見下圖：43.532.521.510.50-0.500.10.20.30.40.50.60.70.80.91Newton 法問題一的牛頓法的求解只適用于線性方程，所以在問題二中用其他方法來求解方程。建立newton1.

15、m源程序，源程序代碼為：function x=newton1(fn,dfn,x0,e) if nargin<4,e=1e-4;endx=x0;x0=x+2*e;i=1 ;while abs(x0-x)>ei=i+1 ；x0=x;x=x0-feval(fn,x0)/feval(dfn,x0);end在matlab軟件中執(zhí)行下列語句并得到最終結(jié)果截圖» fun= inline (" sin (x) +2*x* 2-hc-O. 3356,)；» dfun=inline(J cos;» newt on! (fun, dfun (L I, le-8)i

16、 -5ans -0.1466割線法利用問題一中的割線法程序，在matlab軟件中執(zhí)行下列語句» gexiandurij Qj Q, 2,1 e-3)ans =0.1466拋物線法利用問題一中的拋物線法程序,在 matlab軟件中執(zhí)行下列語句>> paowuxian (fun, 0, 0. 1, 0. 2t 1 e-8)0.1466問題三問題描述按照題目要求對五次方程進行構(gòu)造為：(x3 +2x2 +x-0.3356)(x-0.224)2 =0問題分析仍然利用一中方法求解這一問題，并利用圖解法找到初值，通過觀察圖像， 將牛頓法的兩組初值為0以及0.5,割線法初值設(shè)為：0,0.

17、5以及0, 0.2。拋物線 法初值設(shè)為：兩組初值為 0,0.5,1以及0,0.25,05牛頓法利用問題二中的程序，在 matlab軟件中執(zhí)行下列語句：初值為0時> > fun=inline C (x 3+2*k 2+x-0. 3356)224 2');> > dfun= inlineC 2* (i-3+2*xR2-hc-0, 3356) * (x-Q. 224) + (3*x*2+4*x+l)* (x-0. 224) *Z );> > neirtonl (fuiij dfun, 0, le-8137ans =0. 2240初值為0.5時>>

18、; newt on 1 (fdfuiij 0.1 e-S)33ans =0. 2240割線法利用問題一中割線法的程序，在 matlab軟件中執(zhí)行下列語句: 初值為0, 0.5時» funFuilineC24*L-0, 3356) *(x-0, 224) ,2,» g exian (fun., Oj 0. 5j le-8)50ans =0. 2240初值為0, 0.2時>> fun=inline(, (x"3+2*x*2tK-0, 3356) * (xO. 224) *2");» gezianCfurij 0, 0, 2, le-8)

19、i -44ans -0. 2240拋物線法利用問題一中拋物線法的程序，在matlab軟件中執(zhí)行下列語句:初值為0, 0.5, 1時» fun=inlineC (x"3+2*x 2+x-0,3355)221);» paowuxian(funj 0, 0. 5；, lj le-B)i =63ans -0. 2240初值為0, 0.25, 0.5時» fun=inlineC (x'3+2*x 2-tx-O. 33S6) * (x-0. 224) 2J);» paowwcianCfurij 0j 0. 25, 0. 5,1 e-8)i -51a

20、naQ. 224Q問題四問題描述根據(jù)題目要求對方程進行構(gòu)造為:(sin x 2x2 x-0.3356)( x-0.224) =0問題分析仍然利用問題一中方法求解這一問題，并利用圖解法找到初值，通過觀察圖像， newton法初值選取了兩組初值為0以及0.5,割線法初值設(shè)為：0, 0.5和0, 0.3 拋物線法初值設(shè)為：兩組初值為0,0.2,0.4以及0, 0.5, 1。0.080.070.060.050.040.030.020.010-0.0100.050.10.150.20.250.3牛頓法利用問題二中的程序，在 matlab軟件中執(zhí)行下列語句:初值為0時» fun=inline C

21、 (sin(x)+2xx 2+-Hc-0. 3356) * (x-0. 224)' )：» dfun=inlinef (cos (k)+4*x+1)*(s-0. 224) + (sin(x)+2*s 2+so, 335E)');>> nevtonl (fun, dfun, 0： le-8)ans =0. 1466初值為0.5時» fun=inline(n (sin(x)+2*x"2i la 0. 3356) * (xO. 224)1);> > (ifun=inline C (cos (x)+44x+l) *(x-0, 224

22、) + (sin (x)+2*x''2+x-U4 3356)：> > nevtuni (fun, dfun, 0. 5, le-E)ans -0.2240割線法利用問題一中割線法的程序，在matlab軟件中執(zhí)行下列語句:初值為0, 0.5時» funpijtiline C (sin(x)+2*x*2-Hi-O4 3356) * tf-O* 224)1);> > gesianCfun, 0, 0, 5j le-8'113ans =初值為0, 0.3時» fun=inline C (sin(x)+2*x*2+x-0. 3356)*

23、(x-0.224)r)：>> gexian(fun, . 0. 3, le-S)i =10獨3 =0. 2240拋物線法利用問題一中拋物線法的程序，在matlab軟件中執(zhí)行下列語句:初值為0, 0.1, 0.2時>> fun=inline (? (sin(x)+2*x 2+-hc-04 3356) * tx-0. 224), )：>> paowuxianffun, 口，0. 1, 0. 2, 1 e-8)i =14ans -0,1465 + 0.00001初值為0, 0.5, 1時>> funFinline( (sin+2*x 2+x-0

24、7; 3355)*(x-04 224)：>> paowuxian(fun, 0, 0. 5, 1,1 e-8)14ans -五、計算結(jié)果及分析Newton 法割線法拋物線法問題一0.2240.2240.224問題二0.14660.14660.1466問題三0.2240.2240.224問題四0.224 和0.224 和0.224 和0.1465772585571010.14660.1466問題一迭代步數(shù)問題二迭代步數(shù)問題三迭代步數(shù)問題四迭代步數(shù)Newton 法653789割線法855013 10拋物線法956314 13結(jié)果分析將Newton法，割線法，拋物線法進行比較可以看到在本

25、文題中，三種方法計算得到的最終結(jié)果基本相同，但是迭代步數(shù)有較大差別，綜合看來Newton法迭代步數(shù)最少，割線法次之，拋物線法最次。在各個問題的研究中，我通常都會 采用不同的初值，發(fā)現(xiàn)不同初值會對應(yīng)不同的迭代次數(shù)， 另外針對問題一，我選 用了不同的精度，發(fā)現(xiàn)迭代次數(shù)并沒有很大的變化，因而一個初值的選定可以對 迭代次數(shù)產(chǎn)生很大的影響，而精度的改變對迭代次數(shù)的影響很小。對每個算法單獨來看，顯然選擇初值不同對于迭代步數(shù)影響較大，對于找到 根也會有影響。因此應(yīng)該先通過畫圖確定根的大致位置，給出在其附近的初值。六、心得體會在實現(xiàn)這三個算法的過程中，本身編程較易實現(xiàn)，最重要的是對算法本身的 理解，只有真正理

26、解算法的含義才能更快更好的實現(xiàn)程序。II、離散型最小二乘和連續(xù)型最小二乘問題一、實驗?zāi)康恼莆詹⒛軌蚶秒x散型最小二乘和連續(xù)性最小二乘求解問題、問題描述1:以函數(shù)錯誤！未找到引用源。生成的6個數(shù)據(jù)點(0.25,23.1) ,(1.0,1.68), (1.5,1.0), (2.0,0.84) , (2.4,0.826) , (5.0,1.2576)為例，運行程序得到不同階對應(yīng)的曲線擬合產(chǎn)生的多項式函數(shù)。2.例題：計算f(x)=exp(x)在卜1, 1上的二、三次最佳平方逼近多項式。并畫圖進行比較。三、問題分析問題1是離散最小二乘問題最小離散最小二乘就是根據(jù)一批有誤的數(shù)據(jù)(錯誤！未找到引用源。i=1

27、,2,，n確定參數(shù)，并通過均方誤差來比較曲線擬合的優(yōu)劣，在本題中通過 畫圖來比較不同階方程擬合效果的優(yōu)劣。選擇兩種方法實現(xiàn)離散最小二乘。方法一，建立normal equation(法方程組),求解k次多項式系數(shù)。法方程組構(gòu) 造方法：m0a。' xi =0mnan"xi=0m_ 、 0- yi xii =0mm-1-2a0xi .a xi -i =0i =0m_ n，1anxii旬m1二 xxi =0ma°v xini =0mm- a/' xin 1 ' . aj xi2ni =0i =0m=£i =0nYiX方法二：由于在matlab中存在

28、ployfit函數(shù)，可以即為方便的用k次多項式擬問題2是一個連續(xù)型最小二乘法的應(yīng)用實例。對于給定的函數(shù)f(x) w Ca, b如果存在S (x) Span( 0(x), 1(x), IH, n(x)使得2b2(P(x) f (x) -S (x) I dx = min f P(x) If (x) -s(x) Tdxa<<b一.一一4八一，、則稱$(刈是£(刈在集合Span(50(x), %(x), | * (x)中的最佳平萬逼近函數(shù)。n顯然，求最佳平萬逼近函數(shù)S (x) =£ aj ,p(x)的問題可歸結(jié)為求它的系數(shù)j 0a；,a；,，a；,使多元函數(shù)bnI (a

29、。，ai, " , a；)=dx(P(x) |f (x) -Z aj啊(x)a-1/一取得極小值，也即點(a0, a1*,，a：)是I (a。，an)的極點。由于I (a。，ai,an)是關(guān)于a。, ai,，an的二次函數(shù)，利用多元函數(shù)取得極值的必要條件,=0 (k = 0, 1,2, , n) a即Hbn =2 ：(x) f(x) -y aj j(x) L ：k(x)dx=0-aka |L «得方程組n b'、aj ?(x) (x) j(x)dx aj=0b=7(x)f (x) k(x)dx, (k=0,1,2, ,n) a如采用函數(shù)內(nèi)積記號b促k,%)= P(x

30、)%(x)Q(x)dx,q(f, k) =:?(x)f(x) k(x)dx,a那么，方程組可以簡寫為(Dn、 (),。問=(f, 1) (k 0,1,2,|l|,n)j f這是一個包含n + 1個未知元a0, ai,，an的n + 1階線性代數(shù)方程組，寫成矩陣形式為1'(%wjWH)iHi,%)WiWi)llllllllllll ©,%)(中 n,*)IIIIIIIll促呼n)、 促 1,Q)(中 nWn)Ja0ai、(f,九),(f,Q)(2)此方程組叫做求aj (j = 0, 1,2,，n)的法方程組。顯然，其系數(shù)行列式就是克萊姆行列式 Gn = Gn (<p0,物

31、，<Pn)o由于邛0,明，旺線性無關(guān)，故Gn0 0 ,于是上述方程組存在唯一解ak =ak(k =0,i，n)。從而月止了函數(shù) f (x)在 SpaK%(x),%(x),| Wn(x)中如果存在最佳平方逼近函數(shù)，則必是.(3)n_ *一* 2S(x)c aj j(x)j =0四、實驗程序問題一法一：離散最小二乘法function f=normalequation(n)syms tr=0;xx=0.25,i.0,i.5,2.0,2.4,5.0;yy=23.i,i.68,i.0,0.84,0.826,i.2576;x=zeros(n,n);y=zeros(n,i); for i=i:nfor

32、 k=i:6x(i,i)=xx(k(2*i-2)+x(i,i);y(i,i)=yy(k).*xx(k(i-i)+y(i,i);endendfor i=2:nfor j=i:i-ifor k=i:6x(i,j)=xx(k).A(i+j-2)+x(i,j);x(j,i)=x(i,j);endendendz=xy;for i=i:nend vpa(r,5)法二:用matlab中的ployfit函數(shù)對k次多項式進行擬合。建立number2.m文件，帶入如下：x=0.25,1.0,1.5,2.0,2.4,5.0;y=23.1,1.68,1.0,0.84,0.826,1.2576; plot(x,y,&#

33、39;o') p1=polyfit(x,y,1) %利用線性擬合 xi=-0.25:0.01:6.0;y1=polyval(p1,xi);plot(x,y,'o',xi,y1,'r:');hold on;p2=polyfit(x,y,2) % 二次擬合y2=polyval(p2,xi);plot(x,y,'o',xi,y2,'m');hold on;p3=polyfit(x,y,3) % 三次擬合y3=polyval(p3,xi);plot(x,y,'o',xi,y3,'b:');hold

34、on;p4=polyfit(x,y,4) % 四次擬合y4=polyval(p4,xi);plot(x,y,'o',xi,y4,'c');hold on;p5=polyfit(x,y,5) % 五次擬合y5=polyval(p5,xi);plot(x,y,'o',xi,y5,'g'); hold off;legend('初始點,'y=-2.3840*x+9.8499','初始點,'y=2.1793*xA2-15.8845*x+22.3999','初始點','y

35、=-2.0143*xA3+18.3499*xA2-45.2167*x+32.7386','初始點','y=1.4828*xA4-16.0435*xA3+56.3154*xA2-79.4994*x+39.6688',' 初始點 ','y=-0.9961*xA5+12.1743*xA4-55.2367*xA3+116.6262*xA2-116.7266*x+45.8090');問題二：最佳平方逼近算法(1)輸入被逼近函數(shù)f(x)和對應(yīng)的逼近區(qū)間a,b并選擇逼近函數(shù)系q(x)和權(quán)函數(shù)；(2)解方程組(1)或(2),其中方程組的系

36、數(shù)矩陣和右端的項由式(3)得至(3)由式(3)得到函數(shù)的最佳平方逼近。將上述算法編寫成MATLAB程序共需三個程序:第一個程序(函數(shù)名 Sapproach.m)計算最佳逼近函數(shù)的系數(shù)：function S=Sapproach(a,b,n)% 定義逼近函數(shù)global i;global j;if nargin<3 n=1;end% 判斷X=zeros(n+1,n+1);for i=0:nfor j=0:n;X(i+1,j+1)=quad(fan,a,b); % 求 fan 積分 end endY=zeros(n+1,1);for i=0:nY(i+1)=quad(yb,a,b);%求 yb

37、 積分end s=XY第二個程序(函數(shù)名：fan.m)function y=fan(x) global i;global j;y=(poly(x,i).*poly(x,j);；第三個程序(函數(shù)名：yb.m)function y=yb(x) global i;y=(poly(x,i).*exp(x);編寫多項式函數(shù)function y=poly(x,k)% 多項式函數(shù)if k=0y=ones(size(x);elsey=x.Ak;end五、實驗結(jié)果及圖像問題一：(1) 最小二乘法 normalequation.m運行結(jié)果為: 線性擬合：>> norjiial equat i on(

38、2)ans 二10. 68 - 2. SllS*t二次擬合：» normalequat ion(3) ans 二2. 5532M、2 - I6.962*t + 22.931三次擬合：» normalequationans =-2.2955*t13 + 19. 506*t-2 - 46. 507*+ + 33. 037四次擬合：» normalequal ion(b)ans =E6303*t'4 - 17. 152*t 3 + 5S. 393*t2 - 80, 932rt + 39. 917五次擬合：» no rm ale quation(6) a

39、ns =-1. 0853M"5 + 13. 027#t4 - 5；. 505+t"3 + 119.36粒 *2 - 113. 14*t + 46. 024用方法二得到的結(jié)果如下：number2.m運行結(jié)果為:pl =-2.9L1910. 6805p2 =2.5532-16.962022. 930919.5064-46. 507433, 03721. 6803 -17. 152253.3927 -80. 932439. 9168p5 =-1. 085313.0258 -57.5051 119. 3595 -118. 139545.0236C因此得到線性擬合多項式為：y - -

40、2.3840X 9.8499二次擬合多項式為：y =2.1793x2 -15.8845x 22.3999三次擬合多項式為：y = -2.0143x3 18.3499x2 -45.2167x 32.7386四次擬合多項式為：y =1.4828x4 -16.0435x3 +56.3154x2 -79.4994x+39.6688五 次 擬 合 多 項 式 為y = -0.9661x5 12.1743x4 -55.2367x3 116.6262x2 -116.7266x 45.8090繪制的各多項式擬合圖像為：50403020100-10-2000,511.522.533.544.55線性擬合 二次擬

41、合 三次擬合 四次擬合 五次擬合 原數(shù)據(jù)點問題2:當求二次逼近時，有以下結(jié)果:>> a approach (-1, lj 2)0. 9963I, 10360. 6367繪制兩者的圖形：> > fun='exp(x)'> > fplot(fun,-1,1)> > hold on> > xi=-1:0.1:1;> > yi=polyval(S,xi);>> plot(xi,yi,'r:')得到以下結(jié)果：當三次逼近時，有以下結(jié)果» sapproach(1, lj 3) s =

42、0.99630,99800,53670.1761繪制圖形如下：> > fun='exp(x)'> > fplot(fun,-1,1)> > hold on> > xi=-1:0.1:1;> > yi=polyval(S,xi);>> plot(xi,yi,'r:')得到如下所得圖形：0 S -C 6-0.40 2 Q Q 20.40.60 8六、結(jié)果分析顯然離散最小二乘中次數(shù)較高擬合的效果較好，但由于本問題中初始點較 少，并不能完全反應(yīng)函數(shù)本身的特點，同時在本問題中，選擇了兩種方式，結(jié)果 接

43、近，也可以互相驗證正確性。在連續(xù)最小二乘法的實驗中較順利的達到了預期的結(jié)果。從試驗結(jié)果看出三次逼近沒有二次逼近效果理想，驗證了最佳平方逼近理論。七、心得體會在離散最小二乘與連續(xù)最小二乘程序設(shè)計的過程中，要么通過均差來表示函 數(shù)擬合的優(yōu)劣，要么通過圖像來評價，本問題中選擇了圖像，更加清晰直觀。III、數(shù)值積分一、實驗?zāi)康氖煜げ⒄莆諗?shù)值積分的方法，重要訓練復化梯形公式，復化 simpson公式以及 romberg 積分。二、問題描述22問題三數(shù)值積分橢圓周長的計算。考慮橢圓與+%=1，為計算其周長，只要計算a b其第一象限的長度即可.一 .、一一. x = a cost用參數(shù)方程可以表示為(0 &

44、lt; t M二/ 2),y = bsint計算公式為.a2 sin21 b2 cos2 tdt0為計算方便，我們可以令a = 1，即計算下面的積分V2 :222sin t b cos tdt=f0 /l +(b2 -1)cos2 tdt/ /2 -212-22 ,(I. ;a sin t b cos tdt0=a(sin21 +(b)2cos tdt可以歸結(jié)為上面的形式 采用復化梯形公式，復化 Simpson公式以及Romberg積分的方法計算積分72：-I (b) = p , 1 (b - 1)cos tdt給出通用程序，該通用程序可以計算任何一個函數(shù)在任意一個區(qū)間在給定的精度下的數(shù)值積分

45、。程序輸出為計算出的數(shù)值積分值以及計算函數(shù)值的次數(shù)。利用你的程序計算5個橢圓的周長。三、算法介紹首先利用給出的各迭代公式，設(shè)計程序。在 matlab對話框中輸入要計算的函數(shù)， 給出區(qū)間和精度。bhn 1旻化梯形的迭代公式為：a f(X)dX=21f(a)2 f (Xj) . f(b)2j 1bhr ,、f)-1 ,、,、,、復化 simpson迭代公式為：f f (x)dx =- f (a)+2£ f(x2j)+4£ f (x2j)+ f (b)a3j 1j 1Romberg迭代公式為： RJ=R+ Rk,jRk"四、算法程序復化梯形公式和復化simpso心式我們

46、放在jifen.m中。(標記后的程序可用來把b看為變量時的算法實現(xiàn))%復化梯形公式function y=jifen(f,n,a,b)(說明：f表示任一函數(shù)，n精度，a, b為區(qū)間)fi=f(a)+f(b);h=(b-a)/n;d=1;%function f=jifen(n,a,b,c)%syms t%y=sqrt(1+(cA2-1)*cos(t)A2);%ya=subs(y,t,a);%yb=subs(y,t,b);%fi=ya+yb;for i=1:n-1x=a+i*h;fi=fi+2*f(x);d=d+1;%yx=subs(y,t,x);%fi=fi+2*yx;endf4=h/2*fi,d

47、%復化simposo心式f1=0;f2=0;dd=1;for i=1:n-1dd=dd+1;if rem(i,2)=0; x1=a+i*h; f1=f1+f(x1);else rem(i,2)=0; x2=a+i*h; f2=f2+f(x2); endendf3=(h/3)*(f(a)+4*f1+2*f2+f(b),ddromberg 積分建立romberg.m文件function y=romberg(f,n,a,b)(說明：f表示任一函數(shù)，n精度，a, b為區(qū)問)z=zeros(n,n);h=b-a;z(1,1)=(h/2)*(f(a)+f(b);f1=0;for i=2:nfor k=1:

48、2A(i-2)f1=f1+f(a+(k-0.5)*h);endz(i,1)=0.5*z(i-1,1)+0.5*h*f1;h=h/2;f1=0;for j=2:iz(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1)/(4A(j-1)-1);endendz,n五、實驗結(jié)果(1)對于復化梯形公式和復化simpson式，我們運行下列語句并得到結(jié)果:>> fun=inline(? sqrt 11+' C. 5 2-1 j. *coe (t). 2 1 );>> j if ent fun, 8j Ojpi/2)f4 =1.2111d -f3 =L2U1d

49、d =8題中將橢圓中的未知量a取為1, b取為0.5。f4為復化梯形公式得到的1/4個橢 圓周長，f3為復化simpson公式得到的1/4個橢圓周長。從而得到橢圓的周長為 4*1.2111=4.8444。(2)對于romberg,運行下列語句并最終得到結(jié)果為：» fun=inline(' sqrt(1+10. S"21). *cos(t), -2)J):» romberg fun, 31 0t pi/0. 5)2 -3. 141600000003. 14163.1416000004. 71245.23606.3756000004.83984. 88234.

50、85874. 8505000084J24,84574.84324. 84304,84290004.84424, 84424.S4414. 84414.84424.6442004.84424. 84424. 94424. 84424. 84424.E4424. 844204.84424, 34424.84424. 84424,84424.64424.84424,8442n =S題中將橢圓中的未知量a取為1, b取為0.5。用romberg算法對橢圓的周長進行 計算的到橢圓的周長為4.8442。六、結(jié)果分析我們計算了當橢圓長軸為1,短軸為0.5時的周長。通過上述三種方 法的計算可以看到，結(jié)果相差不

51、大。根據(jù)橢圓周長的一個計算公式我 們可以得到L=4.283。因此三種方法都較好的接近真值。七、心得體會我們應(yīng)該熟練掌握這三種方法，才能在編程時正確快速的寫出迭代公 式。同時在一種思想的前提下可以尋找多種方法實現(xiàn)算法，互相驗證,從而得到更準確的結(jié)果。IV、newton差值和hermite差值方法一、實驗?zāi)康恼莆詹⒛軌蚶胣ewton差值和herm讓e差值方法解決問題。二、問題描述二/222I(b) = f0 也 +(b2 -1)cos2 tdt22 bcos t上述函數(shù)的導數(shù)為 I'(b)： dt0221 (b - 1)cos t采用三種方法中最好的方法計算這一積分(1)利用數(shù)值積分白方

52、法給出I(b)在b = 0,0.1,0.2,|, 1 (可以直接計算精確值的，用精確值)，用Newton插值方法得到5個橢圓的周長(2)利用數(shù)值積分的方法給出I(b),I '(b)在b = 0,0.1, 0.2川|,1(可以直接計算精確值的，用精確值)，用Hermle插值方法得到5個橢圓的周長(3)選做題：利用I(b)以及導數(shù)更多的值來進行插值，插值誤差會有什么變化？(4)選做題：采用其它的插值方法改進插值的效果。三、算法介紹血定，對于給定的b值都對應(yīng)著一個橢圓，在本問題中用newton®值法和hermae 得到的多項式代替橢圓周長公式中的，1+(b21)cos2t進行積分，

53、首先畫出圖像，選擇初始點。圖彳O勺實現(xiàn)代碼見picture1.m。newtonj雨值法迭代公式：nPn (X)= f Xo 二 f Xo, XiXk ( X - Xo )( X - Xi )(X - Xk)k =1Hermite法迭代公式：2n 1H 2n 1 (X) = f Zo% f Zi Zk( X - Zo) (X - Zk)k 1作圖時建立picture.m文件畫出和其導數(shù)圖像。(注：此圖像為b=0.5時)0.80.60.40.20-0.2-0.40.20.40.60.81.21.41.61.8sqrt(1+(0.5 2-1).*cos(x). 2)0.750./(1-0.75.*c

54、os(x).2).(1/2).*cos(x).*sin(x)根據(jù)圖像，我們選取我們選取0,0.3,0.6,0.9,1.2,1.5為初始點四、算法程序(1)建立 newtondedail.m文件。function z=newtondedai1(f,n)syms xia=zeros(n,n);x=0 0.3 0.6 0.9 1.2 1.5;y=feval(f,x);a(:,1)=y;for i=2:nfor j=2:ia(i,j)=(a(i,j-1)-a(i-1,j-1)/(x(1,i)-x(1,i-j+1);endendt=xi-x(1,1);p=a(1,1);for i=2:np=p+a(i,

55、i)*t;t=t*(xi-x(1,i);endp=collect(vpa(p)(2)建立 hermite3.m文件。function z=hermite3(fn,dfn,n)syms xia=zeros(2*n,2*n);x=0.3 0.6 0.9;xx=0.3 0.3 0.6 0.6 0.9 0.9;y=feval(fn,x);yy=feval(dfn,x);for i=1:nz(1,2*i)=x(1,i);z(1,2*i-1)=x(1,i);a(2*i,1)=y(1,i);a(2*i-1,1)=y(1,i);a(2*i,2)=yy(1,i);endfor i=1:nif 2*i<6a(2*i+1,2)=(a(2*i+1,1)-a(2*i,1)/(z(1,2*i+1)-z(1,2*i);endendfor i=3:(2*n)for j=3:ia(iJ)=(a(iJ-1)-a(i-1J-1)/(z(1,i)-z(1,i-j+1);endendp=a(1.1);t=xi-xx(1,1);for i=2:2*np=a(i,i)*t+p;t=t*(xi-xx(1,i);endp=vpa(collect(p) %p=vpa