数值分析第四章

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