數(shù)值實踐 線性方程組迭代

1、 數(shù)值計算實踐課程設(shè)計報告課題名稱線性方程組的迭代算法探究與比較完成時間姓名班級學(xué)號成績一、實驗?zāi)康模瑑?nèi)容 二、相關(guān)背景知識介紹 三、程序代碼四、數(shù)值結(jié)果 五、計算結(jié)果的分析 六、計算中出現(xiàn)的問題，解決方法及體會一. 實驗?zāi)康募皟?nèi)容1.1實驗?zāi)康膶τ诰€性方程組的求解，主要有直接法與迭代法。對于階數(shù)不太高的線性方程組，用直接法比較方便。但是，隨著科學(xué)技術(shù)的飛速發(fā)展，需要求解的問題的規(guī)模越來越大，求解效率就顯得尤為重要。因此，有必要深入研究等古典迭代法，與此同時，迭代法也為我們在尋求其他快速的迭代法時提供了思路，本次課設(shè)意在探究以上四種方法對于不同線性方程組求解過程中在運行時間，計算誤差等方面的異

2、同，比較各方法間的優(yōu)缺點，從而得出效率較高的算法推廣應(yīng)用。1.2實驗內(nèi)容（1）確定某一線性方程組，判斷該線性方程組是否滿足的運行條件，若同時滿足則進行下一步。（2）輸入線性方程組系數(shù)矩陣，常數(shù)項，迭代初值以及要求的計算精度。（3）將上述值代入迭代算法中進行計算，得出滿足計算精度的最終解，迭代次數(shù)，以及程序運行時間。（4）通過對比分析輸出值中迭代次數(shù)與運行時間的大小以及最終解的誤差，綜合比較后得到分析結(jié)論。二相關(guān)背景知識介紹迭代法是一種逐步逼近的方法，尤其是隨著高速計算機的出現(xiàn)和廣泛應(yīng)用，使得迭代法更加容易實現(xiàn)。求解線性方程組的迭代法具有存儲空間小，程序設(shè)計簡單等優(yōu)點。對于大型線性方程組（尤其是

3、大型工程問題所產(chǎn)生的線性方程組），迭代法更加優(yōu)越。2.1 Jacobi迭代將線性方程組中的系數(shù)矩陣分成三部分設(shè)，選取為的對角線元素部分，即選取（對角矩陣），從而得到解的雅可比迭代法 (1)其中。稱為解的雅可比迭代法的迭代矩陣。2.2 Gauss-Seidel迭代選取分裂矩陣為的下三角部分，即選取（下三角矩陣），于是得到解的高斯-塞德爾迭代法 (2)其中稱為解得高斯-塞德爾迭代法的迭代矩陣。2.3 SOR迭代選取分裂矩陣為帶參數(shù)的下三角矩陣，其中為可選擇的松弛因子。于是，由可構(gòu)造一個迭代法，其迭代矩陣為 (3)從而得到解的逐次超松弛迭代法，簡稱方法。解的方法為 (4)其中2.4 CG迭代設(shè)對稱正

4、定，若中的向量組滿足，則稱它是中的一個，或稱。共軛梯度法算法如下：設(shè)從出發(fā)，已迭代計算至，選取搜索方向 (5)由與關(guān)于，即得，從而搜索方向 (6)再沿方向求的極小點，即，經(jīng)計算得 (7)由此得方程組的近似解序列。三程序代碼 3.1 嚴(yán)格對角占優(yōu)function D=SDDM(A)%A為線性方程組的系數(shù)矩陣m,n=size(A);a=ones(m,1);s=zeros(m,1);D=0;for i=1:m a(i)=abs(A(i,i);%將A中對角線元素取絕對值endfor i=1:m A(i,i)=0;%將A中對角線元素賦值為0endfor i=1:m for j=1:n s(i)=abs(

5、A(i,j)+s(i);%計算新的A中每行元素絕對值的和 endendfor i=1:m if a(i)>s(i) D=1; else D=0; break; endend 3.2 判斷矩陣為零矩陣function c=ZERO(A)m,n=size(A);for i=1:m for j=1:n if(A(i,j)=0)%判斷矩陣A是否為零矩陣 c=0; return end endendc=1;3.3 Jacobi迭代function x,N,R,tol=Jacobi(A,b,eps)%A為方程組的系數(shù)矩陣，b為方程組右端項，x0為初值，eps為精度tic%程序運行起始時間m,n=si

6、ze(A);x0=zeros(m,1);%默認(rèn)初值D=diag(diag(A);%產(chǎn)生A的對角矩陣L=-tril(A,-1);%產(chǎn)生A的下三角矩陣U=-triu(A,1);%產(chǎn)生A的上三角矩陣B=inv(D)*(L+U);f=inv(D)*b;x=B*x0+f;N=1;%迭代次數(shù)R,S=max(abs(eig(B);%計算特征值的最大值，R為最大值即譜半徑D=SDDM(A);%判斷系數(shù)矩陣A是否為嚴(yán)格對角占優(yōu)矩陣if D=1 disp('雅可比迭代收斂'); while(norm(x-x0)>=eps) x0=x; x=B*x0+f; N=N+1; end tol=nor

7、m(x-x0);%誤差估計else if R>1%譜半徑大于1，迭代不收斂 disp('雅克比迭代不收斂'); x=NaN; N=NaN; tol=NaN; return else disp('雅可比迭代收斂'); while(norm(x-x0)>=eps) x0=x; x=B*x0+f; N=N+1; end tol=norm(x-x0); endendtoc%程序運行結(jié)束時間3.4 Gauss-Seidel迭代function x,N,R,tol=GaussSeidel(A,b,eps)%A為方程組的系數(shù)矩陣，b為方程組右端項，x0為初值，ep

8、s為精度tic%程序運行起始時間m,n=size(A);x0=zeros(m,1);D=diag(diag(A);%產(chǎn)生A的對角矩陣L=-tril(A,-1);%產(chǎn)生A的下三角矩陣U=-triu(A,1);%產(chǎn)生A的上三角矩陣B=inv(D-L)*U;f=inv(D-L)*b;x=B*x0+f;N=1;%迭代次數(shù)R,S=max(abs(eig(B);%計算特征值的最大值，R為最大值即譜半徑D=SDDM(A);%判斷系數(shù)矩陣A是否為嚴(yán)格對角占優(yōu)矩陣if D=1 disp('高斯-塞德爾迭代收斂'); while(norm(x-x0)>=eps) x0=x; x=B*x0+f

9、; N=N+1; end tol=norm(x-x0);%誤差估計else if R>1%譜半徑大于1，迭代不收斂 disp('高斯-塞德爾迭代不收斂'); x=NaN; N=NaN; tol=NaN; return else disp('高斯-塞德爾迭代收斂'); while(norm(x-x0)>=eps) x0=x; x=B*x0+f; N=N+1; end tol=norm(x-x0); %輸出誤差 endendtoc%程序運行結(jié)束時間3.5 SOR迭代function x,N,wopt,tol=SOR(A,b,eps,w)%A為方程組的系數(shù)

10、矩陣，b為方程組右端項，x0為初值，eps為精度,w為松弛因子tic%程序運行起始時間m,n=size(A);x0=zeros(m,1);%默認(rèn)初值D=diag(diag(A);%產(chǎn)生A的對角矩陣L=-tril(A,-1);%產(chǎn)生A的下三角矩陣U=-triu(A,1);%產(chǎn)生A的上三角矩陣m,n=size(A);if m=n k=eig(A);%k為矩陣A的特征值else returnendif A'=A %判斷A是否為對稱矩陣 y=1;else y=0;endif k>0 %判斷是否為正定矩陣 z=1;else z=0;end B=inv(D)*(L+U);v,g=eig(B)

11、;%計算雅可比迭代矩陣的特征值，特征向量；g為特征值R,S=max(max(abs(g);%計算特征值的最大值，R為最大值即譜半徑wopt=2/(1+sqrt(1-R2);%計算最優(yōu)松弛因子Lw=inv(D-wopt*L)*(1-wopt)*D+wopt*U);f=wopt*inv(D-wopt*L)*b;x=Lw*x0+f;N=1; if w>0&w<2 if w<=1 %判斷是否滿足嚴(yán)格對角占優(yōu)，同時是否滿足0<w<=1 D=SDDM(A); if D=1 c=1; disp('對角占優(yōu)矩陣'); end else if y&z

12、 %判斷是否為對稱正定矩陣 c=1; disp('對稱正定矩陣'); else if R>1 %判斷譜半徑是否大于1 disp('SOR迭代不收斂'); else c=1; end end endelse disp('SOR迭代不收斂');endif c=1 disp('SOR迭代收斂'); while(norm(x-x0)>=eps) x0=x; x=Lw*x0+f; N=N+1; end tol=norm(x-x0); %誤差估計else x=NaN; N=NaN; wopt=NaN; tol=NaN;endtoc

13、%程序運行結(jié)束時間3.6 CG迭代function x,N,tol=CG(A,b,eps)%A為方程組的系數(shù)矩陣，b為方程組右端項，x0為初值，eps為精度tic %程序運行起始時間m,n=size(A);x0=zeros(m,1);%默認(rèn)初值if m=n k=eig(A);%k為矩陣A的特征值else returnendif A'=A %判斷A是否為對稱矩陣 y=1;else y=0;endif k>0 %判斷A是否為正定矩陣 z=1;else z=0;endN=0;r0=b-A*x0;p=r0;x=x0;R=A*p;if y&z for k=2:100 c=ZERO(

14、r0); if c=0&&dot(p,R)=0 alpha=dot(r0,r0)/dot(p,R); x0=x; x=x+alpha*p; r=r0-alpha*R; belta=dot(r,r)/dot(r0,r0); p=r+belta*p; R=A*p; r0=r; N=N+1; if norm(x-x0)<=eps break; end else break; end end tol=norm(x-x0); %誤差估計else x=NaN; N=NaN; tol=NaN; returnendtoc %程序運行結(jié)束時間四數(shù)值結(jié)果例1：用迭代法和迭代法解線性方程組解：

15、>> A=2 -1 1;1 1 1;1 1 -2;b=1;1;1;eps=1e-6;>> x,N,R,tol=Jacobi(A,b,eps)雅克比迭代不收斂x = NaNN = NaNR = 1.1180tol = NaN>> A=2 -1 1;1 1 1;1 1 -2;b=1;1;1;eps=1e-6;>> x,N,R,tol=GaussSeidel(A,b,eps)高斯-塞德爾迭代收斂Elapsed time is 0.004516 seconds.x = 0.6667 0.3333 0N = 21R = 0.5000tol = 6.7435

16、e-007例2：利用不同的迭代法求解線性方程組解：>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-6;>> x,N,R,tol=Jacobi(A,b,eps)雅可比迭代收斂Elapsed time is 0.001369 seconds.x = -4.0000 3.0000 2.0000N = 24R = 0.5061tol = 8.9360e-007>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-8;>> x,N,R,tol=Jacobi(A,b,eps)雅可比迭代

17、收斂Elapsed time is 0.001581 seconds.x = -4.0000 3.0000 2.0000N = 31R = 0.5061tol = 7.4992e-009>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-6;>> x,N,R,tol=GaussSeidel(A,b,eps)高斯-塞德爾迭代收斂Elapsed time is 0.000797 seconds.x = -4.0000 3.0000 2.0000N = 11R = 0.2000tol = 6.6794e-007>> A=5 2

18、1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-8;>> x,N,R,tol=GaussSeidel(A,b,eps)高斯-塞德爾迭代收斂Elapsed time is 0.000827 seconds.x = -4.0000 3.0000 2.0000N = 14R = 0.2000tol = 4.8364e-009>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-6;w=0.9;>> x,N,wopt,tol=SOR(A,b,eps,w)對角占優(yōu)矩陣SOR迭代收斂Elapsed time is

19、0.001053 seconds.x = -4.0000 3.0000 2.0000N = 15wopt = 1.0738tol = 4.2418e-007>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-8;w=0.9;>> x,N,wopt,tol=SOR(A,b,eps,w)對角占優(yōu)矩陣SOR迭代收斂Elapsed time is 0.001094 seconds.x = -4.0000 3.0000 2.0000N = 18wopt = 1.0738tol = 5.2489e-009>> A=5 2 1;-1 4

20、 2;2 -3 10;b=-12;20;3;eps=1e-6;>> x,N,tol=CG(A,b,eps)不是對稱正定矩陣x = NaNN = NaNtol = NaN>> A=5 2 1;-1 4 2;2 -3 10;b=-12;20;3;eps=1e-8;>> x,N,tol=CG(A,b,eps)不是對稱正定矩陣x = NaNN = NaNtol = NaN例3：利用不同的迭代法求解線性方程組，其中解：>> A=4 3 0;3 4 -1;0 -1 4;b=3;5;-5;eps=1e-6;>> x,N,R,tol=Jacobi(A

21、,b,eps)雅可比迭代收斂Elapsed time is 0.001527 seconds.x = -0.0000 1.0000 -1.0000N = 62R = 0.7906tol = 9.9300e-007>> A=4 3 0;3 4 -1;0 -1 4;b=3;5;-5;eps=1e-6;>> x,N,R,tol=GaussSeidel(A,b,eps)高斯-塞德爾迭代收斂Elapsed time is 0.001013 seconds.x = 0.0000 1.0000 -1.0000N = 28R = 0.6250tol = 9.1401e-007>&

22、gt; A=4 3 0;3 4 -1;0 -1 4;b=3;5;-5;eps=1e-6;w=1.3;>> x,N,wopt,tol=SOR(A,b,eps,w)對稱正定矩陣SOR迭代收斂Elapsed time is 0.000765 seconds.x = 0.0000 1.0000 -1.0000N = 13wopt = 1.2404tol = 2.1651e-007>> A=4 3 0;3 4 -1;0 -1 4;b=3;5;-5;eps=1e-6;>> x,N,tol=CG(A,b,eps)Elapsed time is 0.001865 seconds.x = -0.0000 1.0