迭代法試驗(yàn)報(bào)告_米瑞琪

1、數(shù)值實(shí)驗(yàn)報(bào)告n實(shí)驗(yàn)名稱迭代法求解Poisson方程實(shí)驗(yàn)時(shí)間2016年4月10日姓名米瑞琪班級信息1303學(xué)心1309010304成績一、實(shí)驗(yàn)?zāi)康模瑑?nèi)容1、理解并掌握三類迭代方法(Jacobi, Gauss-Siedel,SOR)迭代方法的構(gòu)造原理：2、了解矩陣在matlab中不同的存儲方式會導(dǎo)致不同的計(jì)算效率，學(xué)會采用最高效的方式存儲矩陣；3、學(xué)會在計(jì)算機(jī)上實(shí)現(xiàn)迭代法，并比較不同方法的效率與誤差;二、算法描述(一)Jacobi 迭代對于線性方程組：Au = b (1)為了構(gòu)造迭代格式，可將上式改寫為：u = Tu + d (2)假定A的對角元均不為0,將A分裂為：A = D-B (3)其中D為

2、：D = diag(Aii，4zz，433,，4”)(4)則(1)式可寫為：Du = Bu + b (5)形成Jacobi迭代：uk+1 =+ D-lb (6)為了在迭代格式中不出現(xiàn)B,將B=D-A代入(6)式有：uk+1 = (I-。)1? + D-lb (7)為了保證迭代格式的收斂，需要使迭代矩陣的譜半徑p(T) 1.(二)SOR迭代對于(1)式中的線性方程組，將A分裂為：A = D-L-U (8)其中L, U分別為A對應(yīng)的上三角與下三角矩陣的負(fù)值。由Gauss-Siedel迭代格式的中間步驟：D1 聲 1 = Luk+1 + Ruk + f (9)引入非零參數(shù)3作為松弛因子，則有：Uk+

3、1 = / + 3(小+1 - #)(10)將(9)式代入(10)式，整理之后可以得到SOR迭代格式：uk+1 = (D - a)L)-1(l- 3)。+ a)Ruk + (D -儂工廠1“/ (11) 計(jì)算可得最佳松弛因子為：opt2i +Vi-p2(12)(13)(14)其中_p = cosnh(三)Gauss-Siedel 迭代在SOR方法中令松弛因子3 = 1,則有114+1 = #+1則有迭代格式：Duk+1 = Luk+1 + Ruk +/(15)可見Jacobi迭代與Gauss-Siedel迭代的區(qū)別在于Jacobi迭代將A分解為D,L+R, Gauss-Siedel 迭代將A分

4、解為D-L, Ro三.程序代碼按照上述三種格式分別在matlab中實(shí)現(xiàn)，代碼如下：(一)Jacobi 迭代說明：u為最終的運(yùn)算結(jié)果，A為系數(shù)矩陣，tol為允許的最大迭代步數(shù)，steps記錄 當(dāng)前已經(jīng)經(jīng)過的迭代步數(shù),eps為規(guī)定的計(jì)算精度，uO為初值 huiction u,steps = Jacobian iteration( A,uO,b,cps,tol )%u is tlie final result, A is the iteration matrix， tol is tlie pcmiittcd biggest step,%steps shows how many steps of th

5、e iteration,eps is the precision of tlie%resiilt%u0為初值ifnargin5tol=100000;endifnargin=l)fjprintf(invalid iteration:tlie spectral radius is larger tliaii 1);u=uO:k=0;rctiini;end%d=Db;ciT=A*uO-b;%誤差的計(jì)克采用前后兩步迭代之間向量的差，提高計(jì)算效率wliile max(abs(eiT)eps&stcps=tol)fjMintitcration exceeds limitjacobian);stepsu=u

6、O;return;endu=iiO;end(二)SOR迭代fiuiction u,steps = SOR_Itcration( A,uO,b,eps,tol )if nargin5 tol=10000;endif nargincps&stcps=tolfpnntfCiteration exceeds limitation,SOR);endend以上代碼給出的是一般的SOR迭代方法，具有一般性。實(shí)際上針對Poisson方程，mu是可以預(yù)先計(jì)第 出來的。(三)Gauss-Siedel格式fiuiction u,steps = Gauss_Sicdel Jtcration( A,uO,b,eps,t

7、ol )if nargin5 tol=10000;endif nargincps&stepstol fprintf(fsteps exceed limit,G-S*) return;end end (四)主程序 clc clear %網(wǎng)格剖分 xa=0;xb=l;N1 =64 ;hl =(xb -x a)/Nl;x=xa+l:(Nl-l)*hl;ya=0;yb=l;N2=64;li2=(yb-ya)/N2; y=ya+l:(N2-l)*li2; %gencrate matrix A e2=ones(Nl-l,l);Kl=spdiags(e2,-2*e2,e2,-l,04JSn-lJSri-l)

8、; c3=ones(N2-l,l);Il =spdiags(e3,02-1 JST2-1); A=kion(Kl ,11);A=-A/hlA2;%generate Matrix B c4=oncs(Nl-l,l);I2=spdiags(c4,0Tl-lJSl-l); e5=-2*ones(N2-l ,1);e6=oncs(N2-l,l);K2=spdiags(c6,e5,c6,-l,0,lJSr2-l,N2-l); B=kron(I2JC2);B=B/h2 八 2; %gcnerate g %generate tlie vector of function f f=(x,y)-(2*piA2)

9、*exp(pi*(x+y)*(sin(pi*x)*cos(pi*y)+cos(pi*x)*sin(pi*y); fvec=ones(Nl-l)*(N2-l),l);fori=l:N14for j=l:N2-lf_vec(i-l)*(N2.1)4-j)=f(x(i),y(j);endendX,Y=incshgnd(x,y);uO=zcros0cngtli(fvcc),l);%分別采用三種送代方法求解緊湊格式的近似解ul,stepl=Jacobiaii_Iteration(A+B,uO,fvec);u2,step2=Gauss_Siedel_Itcration(A+BjiO,fvec);u3,st

10、ep3=SOR_Itcration(A+B,uO,fvec);%calciilatc tlic en orlire al= (x ,y) cxp(pi*(x+y)* sin(pi*x).* sin(pi *y);fori=l:Nl-lu_m(i-l)*(N2-l)+l:i*(N2-l)=ii_real(x(i),y);endii_v=ii_inr;ciT_J=inax (ab s(ul -u_v);ciT_G=inax (abs(u2-u_v);en _S=inax (ab s (ii3 -u_v);%月計(jì)算結(jié)果顯示在一張表中res_m=l stcpl err_J;2 stcp2 cn_G;3

11、 stcp3 cit_S;用 rintffMcthod: 1 Jacobiai】2GS3 SOR);fprintethod step error1); resinsol=rcshapc(ul J42-1 mesh(X,Y,sol)四.數(shù)值結(jié)果針對課本上93頁的問題3,程序運(yùn)行的結(jié)果為：表1三種迭代方法性能比較stepserrorJacobian68160. 035788G-S37280. 032068SOR1810. 028627從數(shù)值運(yùn)行結(jié)果可以看出，在誤差基本相同的情況下，三種迭代方法的收斂速度存在較大差異，其 中Jacobian迭代法需要迭代6846次，約為G-S的二倍，而SOR方法只迭

12、代了 184步就得到了具有 更小誤差的解，可見SOR方法在構(gòu)造上比較復(fù)雜，但是實(shí)際運(yùn)行效果卻是最佳的。最終迭代結(jié)果可 以繪制出圖像：圖1步長為1/64時(shí)利用迭代法的計(jì)算結(jié)果針對在黑板上布置的問題，其運(yùn)行結(jié)果為:迭代步數(shù)誤差Jacobian50910. 010989G-S29520. 005419SOR1730. 00047從運(yùn)算結(jié)果上來看，Jacobi迭代方法與GS的步數(shù)之間仍然具有:倍關(guān)系，SOR迭代的收斂速度更 快并且具有更高的精度。五.計(jì)算中出現(xiàn)的問題，解決方法及體會在計(jì)算過程中，困擾我的最大問題就是程序運(yùn)行的效率很低，在最開始一個(gè)程序需要運(yùn)行1到2 分鐘的時(shí)間，后來在同學(xué)的幫助下我們找出了以下問題：矩陣的存儲方式：采用稀疏存儲spdiags與diag相比會大大提高計(jì)算效率，這是因?yàn)閐iag中的0元素參與運(yùn)算導(dǎo)致程序運(yùn)行效率低下：誤差的計(jì)算方式：最初我采用的誤差計(jì)算方式是用en-Aiib進(jìn)行計(jì)算，這種誤差的計(jì)克方式導(dǎo)致遲遲無法收斂