昆明理工數(shù)值分析大作業(yè)線性方程組的迭代法資料

1、 工科研究生數(shù)值分析上機實驗實驗報告題目：線性方程組的迭代法學院：專業(yè)：年級：201 級學生姓名：學號：任課教師： II1、 實驗?zāi)康?1、 通過上機計算體會迭代法求解線性方程組的特點，并能和消去法比較； 2、運用所學的迭代法算法，解決各類線性方程組，編出算法程序； 3、體會上機計算時，終止步驟 或 k大于所給的迭代次數(shù)，對迭代法斂散性的意義； 4、體會初始解x0 ,松弛因子的選取，對計算結(jié)果的影響。 二、問題簡述： 分別選用Jacobi 迭代法，Gauss-Seidol迭代法和SOR方法計算下列線性方程組的解。 方程組一:方程組二(設(shè)對稱正定陣系數(shù)陣線方程組)：方程組三(三對角形線性方程組

2、)：3、 本題目用到的三種基本解線性方程組的迭代法簡述：設(shè)有線性方程組：Ax=b x為n維向量 A為行列標分別為i,j的元素aij素構(gòu)成的n階系數(shù)矩陣Jacobi 迭代法：將A分解為三部分：A=-=D-L-U.Jacobi 迭代法之迭代公式如下：x(k+1)=D-1(L+U)x(k)+D-1bGauss-Seidol迭代法迭代公式:x(k+1)=(D-L)-1Ux(k)+(D-L)-1bSOR方法迭代公式:()x(k+1)=(D-L)-1(1-)D+U)x(k)+(D-L)-1b其中為SOR方法權(quán)因子。根據(jù)實驗要求設(shè)計計算流程圖：Jacobi 迭代法：B= D-1（L+U）l 將系數(shù)矩陣A三角

3、分解為L，U和D，即A=D-L-Ul 設(shè)置迭代值x（k-1），x（k）l 設(shè)置初始迭代向量x0量l 選取精確度。l 定義迭代次數(shù)k=1l 讀取b迭代系數(shù)矩陣譜半徑小于1 否 是帶入初次迭代向量x0 .得到下一次迭代的值：x（k）| x（k）- x（k-1）|<求范數(shù)| x（k）- x（k-1）|;記錄迭代次數(shù)k=k+1;x（k） = D-1（L+U）x（k-1）+D-1b帶入上一次次迭代向量x（k-1） . 否 是輸出最終計算值x*=x（k）;以及終迭代次數(shù)k; 終止；Gauss-Seidol迭代法l 將系數(shù)矩陣A三角分解為L，U和D，即A=D-L-Ul 設(shè)置迭代值x（k-1），x（k）

4、l 設(shè)置初始迭代向量x0量l 選取精確度。l 定義迭代次數(shù)k=1l 讀取b：B= (D-L)-1U迭代系數(shù)矩陣譜半徑小于1 否 是帶入初次迭代向量x0 .得到下一次迭代的值：x（k）| x（k）- x（k-1）|<求范數(shù)| x（k）- x（k-1）|;記錄迭代次數(shù)k=k+1;帶入上一次次迭代向量x（k-1） .x（k） = (D-L)-1Ux(k-1)+(D-L)-1b 否 是輸出最終計算值x*=x（k）;以及終迭代次數(shù)k; 終止；l 將系數(shù)矩陣A三角分解為L，U和D，即A=D-L-Ul 設(shè)置迭代值x（k-1），x（k）l 設(shè)置初始迭代向量x0量l 選取精確度及權(quán)因子、。l 定義迭代次數(shù)

5、k=1l 讀取b SOR方法：B= (D-L)-1(1-)D+U)迭代系數(shù)矩陣譜半徑小于1 否 是帶入初次迭代向量x0 .得到下一次迭代的值：x（k）| x（k）- x（k-1）|<求范數(shù)| x（k）- x（k-1）|;記錄迭代次數(shù)k=k+1;帶入上一次次迭代向量x（k-1） .x（k） = (D-L)-1(1-)D+U)x (k-1)+(D-L)-1b 否 是輸出最終計算值x*=x（k）;以及終迭代次數(shù)k; 終止；五、實驗內(nèi)容：本次實驗共使用三種迭代方法采用不同精確度和相關(guān)因子對三種不同類型的線性方程組進行了上機編程計算。除了在不同精度下采用三種迭代方法進行計算不同方法同一方程的迭代次

6、數(shù)不同，不同精度同一方程的迭代次數(shù)也不同。其中SOR方法因權(quán)因子選取不同而迭代次數(shù)不同，本實驗采用同一精確值分別對三個方程對應(yīng)選取不同的三組權(quán)因子進行了計算。并且采用二分法逼近的方法尋找了相對最佳權(quán)因子。六、計算結(jié)果：通過編程和迭代得到以下計算結(jié)果Jacobi 迭代法：方程一：方程一用Jacobi 迭代法計算時，迭代系數(shù)矩陣譜半徑大于1，故迭代法發(fā)散。方程二：方程三：分析：Gauss-Seidol迭代法：方程一：方程一用Gauss-Seidol迭代法：算時，迭代系數(shù)矩陣譜半徑大于1，故迭代法發(fā)散。方程二：方程三：分析：SOR 方法：方程一：方程一用SOR 迭代法計算時,以0.001，從(0,2

7、)選取權(quán)因子 ,迭代系數(shù)矩陣譜半徑均大于1，選權(quán)因子為0，=1，且隨著增大程增大趨勢，故判斷迭代法發(fā)散。從而方程一用此三種方法均無法計算！方程二：方程三：選取迭代向量初值：w=0.8w=0.9w=1.0e=e-1e=e-2e=e-3e=e-1e=e-2e=e-3e=e-1e=e-2x12.0082.01842.00182.00152.0122.00122.00091.9995x20.88321.02871.00181.00131.01251.00071.00021.0063x3-3.087-2.9782-2.9996-3-2.9936-3.0008-3.001-3.0052x40.04990.

8、0047-0.0016-0.0016-0.006-0.0021-0.0019-0.0076x51.06340.98020.99710.99750.98320.99760.99820.9909x6-1.9766-2.0371-2.003-2.0024-2.0175-2.0017-2.0011-2.0027x73.00682.96722.99792.99842.99032.99922.99953.0002x8-0.0182-0.0151-0.0011-0.0008-0.0023-0.0003-0.00010.0005x90.95240.99830.99960.99971.00030.999911.

9、0003x10-1.0335-0.9989-1.0001-1.0001-0.9997-1-1-0.9999迭代次數(shù)3 5825726w=1.0w=1.1w=1.2e=e-3e=e-1e=e-2e=e-3e=e-1e=e-2e=e-3x11.99981.99962.00431.99791.99741.99251.9989x20.99940.99940.99411.00061.0011.00651.0002x3-3.0008-3.0007-2.998-3.0004-3.0005-3.0036-3.0004x4-0.0006-0.0003-0.00140.00060.00080.00210.0002

10、x50.999811.00230.99980.99960.99760.9999x6-2-2-2.0006-2-2-1.9989-2x7332.9999333.00013x800000-0.00010x91111111x10-1-1-1-1-1-1-1迭代次數(shù)7267279分析：七、結(jié)果分析八、實驗結(jié)論附件： 三種迭代法的迭代程序：4Jacobi 迭代法：A=A;b=ba;x0=x0a; %帶已知數(shù)據(jù)(a:aready)D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣B=D(L+U);f=Db;x=B*x0

11、+f;n=; %迭代次數(shù)if vrho(B)<1 k=10e-1 10e-2 10e-3;for eps=k; i=1; n(i)=1; while norm(x-x0)>=eps x0=x; x=B*x0+f; n(i)=n+1; end endelse n=0;endGauss-Seidol迭代法：A=A;b=ba;x0=x0a; %帶已知數(shù)據(jù)(a:aready)D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣B=(D-L)U;f=(D-L)b;x=B*x0+f;n=; %迭代次數(shù)if v

12、rho(B)<1 k=10e-1 10e-2 10e-3;for eps=k; i=1; n(i)=1; while norm(x-x0)>=eps x0=x; x=B*x0+f; n(i)=n+1; end endelse n1=0;endSOR 方法：A=Amn;b=ba;x0=x0a; %帶已知數(shù)據(jù)(a:aleady)D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣w=0;wm=(0.8:0.1:1.2)；x=zeros(m 15);i=1;j=1;n=zeros(5 3);c=zero

13、s(2 15);for w=wm i=i+1; B=(D-w*L)(1-w)*D+w*U); f=(D-w*L)(w*b); x(:,j)=B*x0+f; if vrho(B)<1 k=10e-1 10e-2 10e-3; for eps=k; while norm(x(:,j)-x0)>=eps x0=x(:,j); x(:,j)=B*x0+f; n(i,j)=n1+1; end c(2,j)=w; j=j+1; end else if vrho(B)>=1 x(:,j)=0; c(1,j)=vrho(B); j=j+3; endendEnd clearx=-1*1,1,1

14、,1,1,1,1,1,1,;c=4*ones(1,10);A=diag(c)+diag(x,1)+diag(x,-1); b=7 5 -13 2 6 -12 14 -4 5 -5 'x0=2 1 -3 0 1 -2 3 0 1 -1 'D=diag(diag(A); %求A的對角矩陣L=-tril(A,-1); %求A的下三角陣U=-triu(A,1); %求A的上三角陣w=0;wm=(0.8:0.1:1.2);x=zeros(10 15);i=1;j=1;n=zeros(5 3);c=zeros(2 15);x0=2 1 -3 0 1 -2 3 0 1 -1 '>> >> for w=wm B=(D-w*L)(1-w)*D+w*U); f=(D-w*L)(w*b); x(:,j)=B