1、Jacobi迭代法與迭代法與Seidel迭代法迭代法迭代法的收斂性分析迭代法的收斂性分析超松馳迭代算法超松馳迭代算法平面溫度場實驗平面溫度場實驗數值分析 Ch4 131581079321321321xxxxxxxxx例例4.1特點特點: :系數矩陣主系數矩陣主對角元均不為零對角元均不為零 15/1310/89/7015/115/110/1010/19/19/10)0(3)0(2)0(1)1(3)1(2)1(1xxxxxx計算格式計算格式 X(1)=B X(0) + f 15/ )13(10/ )8(9/ )7(213312321xxxxxxxxx取取 X(0) = 000 X* 1.0000

2、1.0000 1.0000 X(0) 0 0 0 X(1) 0.7778 0.8000 0.8667 X(2)0.96300.96440.9778 X(3)0.99290.99350.9952計算格式計算格式: X(k+1)=BX(k)+f準確解準確解 X(4) 0.99870.99880.9991 nnnnnnnnnnbxaxaxabxaxaxabxaxaxa22112222212111212111雅可比雅可比迭代法迭代法(i = 1,2,n; k=1,2,)取初始向量取初始向量X(0)=x1(0) x2(0) xn(0)T, 迭代計算迭代計算 nijkjijijkjijiiikixaxab

3、ax1)(11)()1(1injjijbxa1(i = 1,2,n)迭代法適用于解迭代法適用于解大型稀疏方程組大型稀疏方程組(萬階以上的方程組萬階以上的方程組, ,系數矩陣中零元素占很系數矩陣中零元素占很大比例大比例, ,而非零元按某種模式分布而非零元按某種模式分布)背景背景: : 電路分析、邊值問題的數值解和數學電路分析、邊值問題的數值解和數學物理方程物理方程問題問題: (1)如何構造迭代格式如何構造迭代格式？ (2)迭代格式是否收斂迭代格式是否收斂？ (3)收斂速度如何收斂速度如何？ (4)如何進行誤差估計如何進行誤差估計？injjijbxa 1(i = 1,2,n)高斯高斯- -賽德爾賽

4、德爾迭代法迭代法 nijkjijijkjijiiikixaxabax1)(11)1()1(1(i = 1,2,n; k =1,2,)取初始向量取初始向量x(0)=x1(0) x2(0) xn(0)T, 迭代計算迭代計算例例 131581079321321321xxxxxxxxx15/ )13()1(2)1(1)1(3 kkkxxx 15/1310/89/700010/1009/19/10115/115/10110/1001)(3)(2)(1)1(3)1(2)1(1kkkkkkxxxxxx9/ )7()(3)(2)1(1kkkxxx 10/ )8()(3)1(1)1(2kkkxxx 15/ )1

5、3(10/ )8(9/ )7(213312321xxxxxxxxx 000)0(3)0(2)0(1xxx雅可比迭代算法雅可比迭代算法A=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;er=1;k=0;while er0.00005 er=0;k=k+1; for i=1:3 s=0;t=x(i);x(i)=0; for j=1:3 s=s+A(i,j)*x(j); end x(i)=t; y(i)=(b(i)-s)/A(i,i); er=max(abs(x(i)-y(i),er); end x=y;xend 131581079321321321xxxxx

6、xxxx0.7778 0.8000 0.86670.9630 0.9644 0.97190.9929 0.9935 0.99520.9987 0.9988 0.99910.9998 0.9998 0.99981.0000 1.0000 1.00001.0000 1.0000 1.0000高斯高斯- -賽德爾迭代算法賽德爾迭代算法 131581079321321321xxxxxxxxxA=9 -1 -1;-1 10 -1;-1 -1 15;b=7;8;13;x=0;0;0;er=1;k=0;while er0.00005 er=0;k=k+1; for i=1:3 s=0;t=x(i);x(i)

7、=0; for j=1:3 s=s+A(i,j)*x(j); end x(i)=(b(i)-s)/A(i,i); er=max(abs(x(i)-t),er); end xend 0.7778 0.8778 0.9770 0.9839 0.9961 0.9987 0.9994 0.9998 0.9999 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000雅可比雅可比 迭代法的矩陣表示迭代法的矩陣表示將方程組將方程組AX = b 的系數矩陣的系數矩陣 A 分解分解 A = D U L nnaaaD2211 0001,121nnnaaaL 000, 1112nnna

8、aaUAX = b = DX(k+1) = (U+L)X(k) + bX(k+1)=D-1(U+L)X(k)+D-1b記記BJ = D-1(U+L) X(k+1)=BJX(k)+fJ雅可比迭代矩陣雅可比迭代矩陣 0002122111212211nnnnnnJaaaaaaaaaB 0/0/02122222211111112nnnnnnnnJaaaaaaaaaaaaB nnnJabababf/222111高斯高斯- -賽德爾迭代法的矩陣表示賽德爾迭代法的矩陣表示 nijkjijijkjijikiiixaxabxa1)(11)1()1( nijkjijiijkjijxabxa1)(1)1(i = 1

9、,2,n) )()(2)(1, 111221)1()2(2)1(121222111000knkknnnnknkknnnnxxxaaabbbxxxaaaaaa (D L)X(k+1) = b + UX(k) X(k+1) = (D L)-1b + (D L)-1UX(k) 記記 BG-S=(D L)-1U, fG-S=(D L)-1b 高斯高斯- -賽德爾迭代格式賽德爾迭代格式: X(k+1)=BG-SX(k)+fG-S 000, 1112121222111nnnnnnnSGaaaaaaaaaB nnnnnSGbbbaaaaaaf21121222111矩陣分裂導出矩陣分裂導出 的迭代法的迭代法A

10、 = M N (要求要求M為可逆矩陣為可逆矩陣)AX =b (M N )X = b MX = NX + b X(k+1) = (M-1N) X(k) + M-1b取取 M = D 雅可比迭代法雅可比迭代法 A =D (D A) X(k+1) = D-1(D A) X(k) + b X(k+1) = X(k)+D-1b AX(k) 思考思考: 取取 將導致什么樣的迭代法？將導致什么樣的迭代法？IM 1 平面點列平面點列: *)(2)(1limyxxxkkk0)()(lim2*2)(22*1)(1 xxxxkkkX (k)Rn : X(1), X (2), , X (k) , *)(limXXkk

11、 0|lim2*)( XXkk*)(limXXkk 0|lim*)( XXkk利用向量范數等價性利用向量范數等價性, 對任意范數對任意范數 | | )1(2)1(1xx )2(2)2(1xx )(2)(1kkxxA X = b (MN )X = b M X = N X + b記記 (k) = X(k) X* ( k = 0, 1, 2, 3, )則有則有 (k+1) = B (k) (k) = B (k-1) ( k = 1, 2, 3, )計算格式計算格式: X(k+1) = B X(k) + f ( B = M-1N ) X(k+1) X*= B(X(k) X*) 設方程組的精確解設方程組

12、的精確解為為 X*,則有則有X* = B X* + f 0lim0lim)( kkkBk (1) (k) = B (k-1)=B2 (k-2)=Bk (0)1|0lim)( Bkk (2)0lim*)( XXkk*)(limXXkk 迭代格式迭代格式 X(k+1) = B X(k) + f 收斂收斂證證: 由由 (k) = B (k-1),得得 | (k)| | B| | (k-1)| ( k = 1, 2, 3, )0lim)( kk 所以所以命題命題 若若|B|1,則迭代法則迭代法 X(k+1) =B X(k) +f 收斂收斂| (k)| | B|k | (0)| 0|lim|lim)0(

13、)( kkkkB| B| 1注注1: AX = b X = BX + f ( I B )X = f X = ( I B )-1 f 注注2: 若若 則則0lim kkB( I - B)-1 = I + B + B2 + + Bk + 事實上事實上 ( I - B)( I + B + B2 + + Bk ) =I Bk+1注注3: X(k) =B X(k-1) + f = B(B X(k-2) + f) + f = = Bk X(0) + ( I + B + + Bk-1)f ( I B )-1 f 定義定義4.1 A=(aij)nn, 如果如果則稱則稱A為嚴格對角占優陣為嚴格對角占優陣. ni

14、jjijiiaa1| . 0)1(, 0)0()1 , 0(yyxxyy例例1 1 常微分方程邊值問題常微分方程邊值問題 求在求在 x1=0.1, x2=0.2, , x9=0.9 處的數值解處的數值解 yj-1 + (2 + h2) yj yj+1 = xj h2 ( j= 1,2,9) 2922219212222112112hxhxhxyyyhhh高斯高斯-賽德爾迭代格式賽德爾迭代格式:212)(1)1(12)1(hxyyhyjkjkjkj 0 0.2 0.4 0.6 0.8 1 0 0.02 0.04 0.06 誤差限設置誤差限設置：10-5。迭代次數迭代次數k=60，error0 =

15、1.2742e-004 1sinhsinh)(xxxy 定理定理4.3 若若Ax=b的系數矩陣的系數矩陣A是嚴格對角占優是嚴格對角占優矩陣矩陣,則則Jacobi迭代和迭代和Seidel迭代均收斂迭代均收斂證證: 由于矩陣由于矩陣A嚴格對角占優嚴格對角占優由由A矩陣構造矩陣構造Jacobi迭代矩陣迭代矩陣BJ = D-1(D A) 第第 i 行絕對值求和行絕對值求和 nijjijiiaa1|1所以所以1 |1max|11 nijjijiiniJaaB1|11 nijjijiiaa nijjijiiaa1|矩陣矩陣B 的譜的譜設設n階方陣階方陣B 的的n個特征值為個特征值為: n ,21則稱集合則

16、稱集合,21n 為為B 的譜的譜. 記為記為 ch B矩陣矩陣B的譜半徑的譜半徑|max)(1knkB 注注1: 當當B是對稱矩陣時是對稱矩陣時, |B|2 = (B) 注注2: 對對 Rnn 中的范數中的范數| |,有有 (B) | B |特征值取模最大特征值取模最大定理定理4.1 迭代法迭代法 X(k+1) = B X(k) + f 收斂收斂 譜半徑譜半徑(B) 1例例2 2 線性方程組線性方程組 A X = b, 分別取系數矩陣為分別取系數矩陣為 1221112211A 2111111122A分析分析Jacobi 迭代法和迭代法和 Seidel 迭代法的斂散性迭代法的斂散性結論結論: A

17、1的的Jacobi迭代矩陣譜半徑小于迭代矩陣譜半徑小于1 收斂收斂 A2的的Jacobi迭代矩陣譜半徑大于迭代矩陣譜半徑大于1 發散發散 A1的的Seidel迭代矩陣譜半徑大于迭代矩陣譜半徑大于1 發散發散 A2的的Seidel迭代矩陣譜半徑小于迭代矩陣譜半徑小于1 收斂收斂Ans= 1.2604e-0051)( JB D=diag(diag(A1);B1=D(D-A1);max(abs(eig(B1)A1=1,2,-2;1,1,1;2,2,1收斂收斂A2=2,-1,1;1,1,1;1,1,-2D=diag(diag(A2)B2=D(D-A2)max(abs(eig(B2)Ans= 1.118

18、01180. 1)( JB 發散發散DL=tril(A1)B1=DL(DL-A1)max(abs(eig(B1)Ans= 22)( SB 發散發散DL=tril(A2)B2=DL(DL-A2)max(abs(eig(B2)Ans= 1/22/1)( SB 收斂收斂定理定理4.2 :設設X*為方程組為方程組 AX=b 的解的解若若|B|1,則對迭代格式則對迭代格式 X(k+1) = B X(k) + f 有有|1|*|)1()()( kkkXXBBXX(1)|1|*|)0()1()(XXBBXXkk (2)證證 由由|B|0, 記記 xTLTx = a , 則有則有xTAx=xT(D L LT)

19、x=p a a =p 2a 02/16xLDxxLxTTT)( xxLLDT 1)(xLDxLT)( 設設 為為BG-S的任一特征值的任一特征值, x 為其特征向量為其特征向量,則則 1)2(2222222 aappaapapa apaxLDxxLxTTT )( 3/16稱稱 R= ln 為迭代法的漸近收斂速度為迭代法的漸近收斂速度)(B 所以所以, 迭代矩陣迭代矩陣 BG-S 的譜半徑的譜半徑 (BG-S) 1,從而當從而當 方程組方程組 Ax=b的系數矩陣的系數矩陣A 是實對稱正定矩陣時是實對稱正定矩陣時, Gauss-Seidel迭代法收斂迭代法收斂)()1()()1(1)()1(bUX

20、LXDXXkkkk )1(1)(11)1()()1( nijkjijijkjijiiikikixaxabaxx (i=1,2, n; k = 1,2,3, )超松馳超松馳(SOR)迭代法迭代法Gauss-Seidel迭代格式迭代格式 nijkjijijkjijiiikixaxabax1)(11)1()1(14/16最佳松馳因子選取最佳松馳因子選取5/162)(112JB 為為Jacobi迭代譜半徑迭代譜半徑)(JB 定理定理4.8 若若 A 是對稱正定矩陣是對稱正定矩陣, 則當則當0 2 時時 SOR 迭代法解方程組迭代法解方程組 Ax = b 是收斂的是收斂的successive overr

21、elaxationProf. David M. Young1954 美國數學科學學報美國數學科學學報Iterative methods for solving partial differential equations of elliptic . 0)1(, 0)0()1 , 0(yyxxyy常微分方程邊值問題常微分方程邊值問題求在求在 x1=0.1, x2=0.2, , x9=0.9 處的數值解處的數值解 yj-1 + (2 + h2) yj yj+1 = xj h2 ( j= 1,2,9) SOR迭代迭代 格式格式2)1()(1)1(122)()1(kjkjjkjkjyyxhhyy hhBJ cos22)(2 2)(112JB Gauss 消元法消元法n誤差誤差51.1902e-00494.4146e-005191.1047e-005294.9106e-006JACOBIn