数值分析（计算方法）讲义4

数数数值值值分分分析析析讲讲义义 适适适用用用范范范围围围:理工科研究生，数学专业本科生 王 震 西京学院 应用数学研究所 2012年6月 王震,男,1981年生,硕士,讲师,陕西省工业与应用数学学会理事,主要研究方向为常微分方程与动力系 统,非线性控制系统理论及应用. Email:williamwangz@yeah.net. Xiji ng Uni v. 王震:数值分析 -1- Lecture 4 Iterative Methods for Solving LEs Ax = b; A 2 Rn�n且非奇异同时，aii 6= 0，令A = D � L� U，其中 D = diag(a11 � � � ann)， L = 26666666664 0 �a21 0 �a31 �a32 . . . ... ... . . . 0 �an1 �an2 � � � �an;n�1 0 37777777775 , U = 26666666664 0 �a12 �a13 � � � �a1n 0 �a23 � � � �a2n . . . . . . ... 0 �an�1;n 0 37777777775 1 Jocobi Iterative Method Ax = b) Dx = (L+ U)x+ b,又aii 6= 0,D�1存在,所以 X = D�1(L+ U)x+D�1b = Tx+ C T = D�1(L+ U)C = D�1b Ax = b,给定一个初始近似值X(0),从而可以得到一个迭代序列 �xfkg 1 k=1 ， 其中 x(k) = Tx(k�1) + C; k = 1; 2; : : : x (k) i = bi aii + nX j=1 j 6=i (�aijx (k�1) j aii ); i = 1; 2; : : : ; n Xiji ng Uni v. 西京学院应用数学研究所 -2- 2 Gauss-Seidel 思思思想想想:对Jacobi迭代的改进，在计算第i个变量x(k)i 时,x(k)1 � � � x(k)i�1已经 计算，故可将Jacobi迭代中的x(k�1)1 � � � x(k�1)i�1 用x(k)1 � � � x(k)i�1代替,即 Dx(k) = Lx(k) + Ux(k�1) + b 故 x(k) = (D � L)�1Ux(k�1) + (D � L)�1b 3 SOR 思思思想想想:在GS迭代中 x (k) 1 = 1 a11 (b1 � a12x(k�1)2 � a13x(k�1)3 � � � � a1nx(k�1)n ) ... x (k) i = 1 aii (bi � i�1P j=1 aijx (k) j � nP j=i+1 aijx (k�1) j ) ... x (k) n = 1ann (bn � an1x (k) 1 � an2x(k)2 � � � � an;n�1x(k)n�1) 即 x (k) i = 1 aii (bi � i�1X j=1 aijx (k) j � nX j=i aijx (k�1) j + aiix (k�1) i ) = x (k�1) i + 1 aii (bi � i�1P j=1 aijx (k) j � nP j=i aijx (k�1) j ) 即 x(k)i = x(k�1)i +4x(k)i 添加松弛因子!，有 x(k)i = x(k�1)i + !4 x(k)i 即 x (k) i = (1� !)x(k�1)i + ![ 1 aii (bi � i�1X j=1 aijx (k) j � nX j=i aijx (k�1) j )] Xiji ng Uni v. 王震:数值分析 -3- 故 x(k) = (1� !)x(k�1) + !D�1[b+ Lx(k) + Ux(k�1)] 即 Dx(k) = (1� !)Dx(k�1) + ![b+ Lx(k) + Ux(k�1)] x(k) = (D � !L)�1[(1� !)D + !U ]x(k�1) + !(D � !L)�1b 4 SSOR (Symmetric Successive Over Relaxation) 思思思想想想:在SOR迭代中，计算x的各分量从第1个到第n个，将次序倒过 来，有 x (k) n = x (k�1) n + !ann (bn � nP j=1 anjx (k�1) j ) ... x (k) i = x (k�1) i + ! aii (bi � iP j=1 aijx (k�1) j � nP j=i+1 aijx (k) j ) ... x (k) 1 = x (k�1) 1 + ! a11 (b1 � a11x(k�1)1 � nP j=2 a1jx (k) j ) 这样将SOR与逆序SOR结合，首先，x(k�1)通过SOR计算x(k)，将x(k)赋 给 x(k�1)，然后运用逆序SOR进行迭代，即可得SSOR. SOR：bx(k) i = x (k�1) i + ! aii (bi � i�1P j=1 aijbx(k)j � nP j=i aijx (k�1) j ) 逆序SOR：x(k)i = bx(k)i + !aii (bi � iP j=1 aijbx(k)j � nP j=i+1 aijx (k) j ) 8<: x^(k) = x(k�1) + !D�1(b+ Lx^(k) + (U �D)x(k�1))x(k) = x^(k) + !D�1(b+ Lx^(k) �Dx^(k) + Ux(k)) EX1：：： 26664 4 3 0 3 4 �1 0 �1 4 37775 26664 x1 x2 x3 37775 = 26664 24 30 �24 37775 写出J、GS、SOR、SSOR迭代格式。 Xiji ng Uni v. 西京学院应用数学研究所 -4- 5 迭迭迭代代代法法法的的的收收收敛敛敛性性性 定定定义义义 1 A 2 Rn�n满足jaiij > nP j=1 j 6=i jaijj; i = 1; 2; � � � ; n(或jajjj > nP i=1 i6=j jaijj)，称A为按行(列)严格对角占优. 定定定义义义 2 A 2 Rn�n满足jaiij > P j 6=i jaijj (或jajjj > P i6=j jaijj)，称A为按 行(列)对角占优. 定定定义义义 3 A 2 Rn�n，当n = 1，A的唯一元素为0，则A可约，否则 不可约.当 n � 2，令N = f1; 2; � � � ; ng，如果存在集合K，K � N;K 6= N，当i =2 K; j 2 K，有aij = 0，则A可约，否则不可约. A 2 Rn�n，如果A不能通过置换矩阵P化为 0@ A11 A12 0 A22 1A，其 中A11为r阶方阵，A22为n� r阶方阵，可称A为不可约矩阵. 定定定理理理 1 由A为按行(列)严格对角占优矩阵可得aii 6= 0，且det(a) 6= 0. Pr: 显然aii 6= 0.用反证法证明det(a) 6= 0. 如果det(A) = 0， 则Ax=0有非零解，令非零解x = (x1; x2; � � � ; xn)T，即存在xk，有 jxkj = max 16i6n jxij > 0 又 ak1x1 + � � �+ aknxn = 0，故 jakkxkj = �������� nX j=1 j 6=k akjxj �������� � nX j=1 j 6=k jakjj jxjj � jxkj nX j=1 j 6=k jakjj 即 jakkj � nP j=1 j 6=k jakjj 矛盾. 定定定理理理 2 A 2 Rn�n，则 Xiji ng Uni v. 王震:数值分析 -5- 1. 8x 2 Rn; lim n!1 Anx = 0 2. lim n!1 An = 0 3. �(A) < 1 9>>>=>>>;条件等价. Pr: 1 =) 2:8x 2 Rn; lim n!1 Anx = 0，不妨取x为Rn的n个单位向 量"1; "2; � � � ; "n，有 lim n!1 An"i = 0,故 lim n!1 AnI = 0，即 lim n!1 An = 0. 2 =) 1：kAnxk 6 kAnk kxk，又 lim n!1 An = 0) 8x 2 Rn; lim n!1 Anx = 0 2() 3： 8A 2 Rn�n,则9P有，P�1AP = J = 0BBB@ J1 . . . Jr 1CCCA 故 A = PJP�1 Ak = PJkP�1,An = P 8>>><>>>: Jn1 . . . Jnr 9>>>=>>>;P �1,又 Ji = 8>>>>>><>>>>>>: �i 1 . . . . . . . . . 1 �i 9>>>>>>=>>>>>>; ni�ni ni为几何重数. Jni = 8>>>>>><>>>>>>: �ni C 1 n� n�1 i � � � Cni�1n �n+1�nii . . . . . . ... . . . C1n� n�1 i �ni 9>>>>>>=>>>>>>; ni�ni 故 lim n!1 An , lim n!1 Jni = 0, j�ij < 1, �(A) < 1 定定定理理理 3 A 2 Rn�n，�(A) < 1) det(I � A) 6= 0 Pr: 如果det(I � A) = 0，则A有特征值为1，显然与�(A) < 1矛盾. J-迭迭迭代代代，，，GS-迭迭迭代代代，，，SOR-迭迭迭代代代，，，SSOR-迭迭迭代代代. 将AX = b! X = BX + f，选取x(0)，有x(k) = Bx(k�1) + f . Xiji ng Uni v. 西京学院应用数学研究所 -6- 定定定理理理 4 x(k) = Bx(k�1) + f，对任意x(0)，有1 迭代() �(B) < 1;2 kBk < 1)收敛. Pr: 1、收敛=) �(B) < 1 如果fx(k)g收敛，设x(k) ! x�，即x� = Bx� + f x� � x(k) = B(x� � x(k�1)) = � � � = Bk(x��x(0))， lim k!1 Bk(x��x(0)) = lim k!1 (x��x(k)) = 0， 由定理2有，�(B) < 1 �(B) < 1 =)收敛 �(B) < 1，由定理3有I �B非奇异，设x^ = (I �B)�1f ,显然x^ = Bx^+ f (因为(B(I �B)�1+ I)(I�B) = B+ I�B = I，所以B(I �B)�1+ I = (I �B)�1） x(k) � x^ = B(x(k�1) � x^) = � � � = Bk(x(0) � x^) 因为�(B) < 1，由定 理2有x(k) � x^ = B(x(k�1) � x^) = � � � = Bk(x(0) � x^) = 0 2 �为B的任一特征值，�为�所对应的特征向量，即B� = ��. j�j = kB�kk�k � maxy2Rn y 6=0 kByk kyk = kBk < 1) �(B) < 1 注注注 1 �(B) > 1不能说迭代格式关于任何向量都不收敛. 定定定理理理 5 对于任意B 2 Rn�n，有�(B) 6 kBk，k�k是任意一种范数 Pr: �为B的任一特征值，�为对应的特征向量,有B� = ��。显 然存在! 2 Rn，有�; !T为非零矩阵，故B�!T = ��!T。由矩阵范数 有j�j � � !T 6 kBk � � !T 即 j�j 6 kBk 又�为B的任一特征值,所以�(B) 6 kBk 显然，如果�(B) 6 kBk < 1 ) �(B) < 1 )迭代格式x(k) = Bx(k�1) + f收敛，对于任意x(0) 2 Rn. 定定定理理理 6 (误差估计) 如果x = Bx+ f有唯一解x�，且kBk < 1，则 x(k) � x(�) 6 kBk 1�kBk x(k) � x(k�1) x(k) � x(�) 6 kBkk�1 1�kBk x(1) � x(0) 9=; Xiji ng Uni v. 王震:数值分析 -7- 且kBk越小，收敛越快。 Pr: x(�) = Bx(�) + f， x(k) = Bx(k�1) + f x(�) � x(k) = B(x(�) � x(k�1)) 6 kBk x(�) � x(k�1) x(k) � x(k�1) 6 kBk x(k�1) � x(k�2) 6 � � � 6 kBkk�1 x(1) � x(0) 所以 x(k) � x(k�1) = x(�) � x(k�1) � x(�) + x(k) > x(�) � x(k�1) � x(�) � x(k) > x(�) � x(k�1) � kBk x(�) � x(k�1) = (1� kBk) x(�) � x(k�1) 所以 x(�) � x(k) 6 kBk 1�kBk x(k) � x(k�1) 且 x(�) � x(k) 6 kBk 1�kBk x(k) � x(k�1) 6 kBkk�1 1�kBk x(1) � x(0) . 定定定理理理 7 (J-迭代收敛判定定理) ① A按行(按列)严格对角占优或A按 行(按列)弱对角占优且为不可约矩阵，则J-迭代法收敛.② A为对称正定 矩阵，则J-迭代的法收敛的充要条件为矩阵2D � A也是对称正定矩阵. ①Pr: 反证法. 若不收敛，则对BJ有特征值j�j > 1 0 = det(�I �BJ) = det(�I �D�1(L+ U)) = det(D�1) � det(�D � L� U) ) det(�D � L� U) = 0 又�D � L� U = 26666664 �a11 �a12 � � � �a1n �a21 . . . . . . ... ... . . . . . . ... �an1 � � � � � � �ann 37777775 A为按行严格对角占优，则j�j jaiij > j�j � nP j=1 j 6=i aij � nP j=1 j 6=i aij 即,�D � L � U为按行严格对角占优，可由定理1知det(�D � L� U) 6= 0矛盾. Xiji ng Uni v. 西京学院应用数学研究所 -8- A为按行多久占优且不可约=) �D � L� U为按行对角占优且不可约. ②A对称正定=) D� 12 (D � A)D� 12对称且有完全正交特征向量组 U1; � � �un，对应特征值�1; � � ��n.8<: D� 1 2 (D � A)D� 12uk = �kuk uTk ul = �kl 令zk = D� 12uk，由BJ = D�1(D � A)且D�1 = D� 12 �D� 12 .8<: BJzk = D�1(D � A)D� 1 2uk = D � 1 2�kuk = �kzk zTkDzl = �kl z1; � � � zn为BJ的完全特征向量组，�1; � � ��n为对应特征值. J-收敛=) �(BJ) < 1 zTk (2D � A)zl = zTk (D +DD�1(D � A))zl = zTkDzl + z T kD�lzl = (1 + �l)�kl 对任意z 2 Rn，有z = nP k=1 ckzk zTk (2D � A)z = nX k=1 nX l=1 ckclz T k (2D � A)zl = nX k=1 c2k(1 + �k) > 0 即2D � A对称正定. 2D � A对称正定=)设�为BJ的任一特征值.存在z 2 Rn有 BJz = �z = D �1(D � A)z 设� = zTDz > 0; 2� = zT(D � A)z = zT(L + LT )z,故� = zT(D�A)z zTDz = 2� � ,又2D � A与A对称正定.有 zT(D � A)z = � + 2� > 0 zTAz = zT(A�D +D)z = � � 2� > 0 即j�j < 1,所以�(BJ) < 1，所以J-迭代收敛. Xiji ng Uni v. 王震:数值分析 -9- 定定定理理理 8 (GS-迭代收敛判定定理) ① A按行(按列)严格对角占优 或A按行(按列)对角占优且为不可约矩阵，则GS-迭代法收敛.② A对称 正定=)GS-迭代收敛. Pr: ①与定理7类似。②BG = (D � L)�1U = (D � L)�1LT 设�为BG的任一特征值，z 6= 0为其对应特征向量，(D � L)�1Uz = �z (考虑z; �可能为复特征向量与复特征值),即LTz = �(D � L)z， 设zHLz = � + i�; zHDz = � > 0 � = zHLT z zH(D � L)z = �� i� � � �� i� 又A对称正定， 0 < zHAz = zHDz � zHLz � zHLT z = � � 2� � � �� = j�j2 = � 2 + �2 (� � �)2 + �2 = �2 + �2 �2 + �2 + (� � 2�)� < 1 定定定理理理 9 SOR迭代收敛，则松弛因子0 < ! < 2. Pr: 设�1; � � ��n为Bs的n个特征值，则 j�1 � � ��nj = jdet(Bs)j = ��det (D � !L)�1((1� !)D + !U)�� = ��det (D � !L)�1�� jdet((1� !)D + !U)j = j(1� !)nj 又 jdet(Bs)j � (�(Bs))n 故 j1� !j < �(Bs) < 1) 0 < ! < 2 定定定理理理 10 (SOR-迭代收敛判定定理) ①A对称正定，0 < ! < 2 )SOR-迭代收敛. ②A按行(按列)严格对角占优 或A按行(按列)对 角占优且为不可约矩阵，0 < ! � 1)SOR-迭代收敛. Xiji ng Uni v. 西京学院应用数学研究所 -10- Pr: ①Bs = (D � !L)�1((1 � !)D + !U),又A对称正定，Bs = (D � !L)�1((1� !)D + !LT),�为Bs的任一特征值，z 6= 0为�的对应特 征向量. (D � !L)�1((1� !)D + LT!)z = �z zH((1� !)D + LT!)z = �zH(D � !L)z 令zHLz = � + i�zHDz = � > 0 � = zH((1� !)D + LT!)z zH(D � !L)z = (1� !)� + !(�� i�) � � !(� + i�) = (� � !� + !�)� i!� (� � !�)� i!� 故j�j2 = (��!�+!�)2+(!�)2 (��!�)2+(!�)2 ,又0 < zHAz = �zHLz + zHDz � zHLTz = � � 2�,由� < 0; 0 < ! < 2有 (� � !� + !�)2 � (� � !�)2 = ��!(2� !)(� � 2�) < 0 即�(Bs) < 1 )SOR-迭代法收敛. ②反证，若不收敛，则�(Bs) � 1，存 在�有j�j � 1且 0 = det(�I �Bs) = det(�I � (D � !L)�1((1� !)D + !U)) = det((D � !L)�1) det(�(D � !L)� ((1� !)D + !U)) ) det((�� 1 + !)D � �!ZL� !U) = 0 设� = � + i�; j�j � 1; 0 < ! � 1) j�� 1 + !j2 � j�!j2 = ((�� 1 + !)2 + �2)� (!2(�2 + �2)) = (1� !)(!(�2 + �2 � 1) + (�� 1)2 + �2) � 0 即 j�� 1 + !j � j�!j Xiji ng Uni v. 王震:数值分析 -11- 令C = (�� 1 + !)D � �!L� !U ,A按行严格对角占优.故 Cii = j(�� 1 + !)aiij � j�!j jaiij > j�!j ( i�1P j=1 jaijj+ nP j=i+1 jaijj) � i�1P j=1 j�!aijj+ nP j=i+1 j!aijj = nP j=1 j 6=i jCijj 即C = (�� 1 + !)D � �!L� !U为按行严格对角占优.) det(C) 6= 0矛 盾，故�(Bs) < 1，SOR-迭代收敛. 定定定理理理 11 (SSOR-迭代收敛判定定理) A对称正定，0 < ! < 2)SSOR-迭代收敛. Pr: Bss = ( D ! � U)�1D(D ! � L)�1(! � 1 ! D � L)D�1(! � 1 ! D � U) W = D�1=2( D ! � U)P = D1=2(D ! � L)�1(! � 1 ! D � L)D�1=2 WBssW �1 = D�1=2(D ! � U)(D ! � U)�1D(D ! � L)�1(!�1 ! D � L)D�1(!�1 ! D � U)(D ! � U)�1D1=2 = PP T W (I �Bss)W�1 = 2� ! ! (W�1)TAW�1 Bss及(I �Bss)对称正定，则�(Bss) < 1)SSOR-迭代收敛. EX1:运用Jacobi迭代，GS迭代，SOR迭代，SSOR迭代计算，并进行收 敛性分析。 8>>><>>>: 2x1 � 1x2 + x3 = 1 x1 + x2 + x3 = 1 x1 + x2 � 2x3 = 1 ; 8>>><>>>: x1 + 2x2 � 2x3 = �3 x1 + x2 + x3 = 1 2x1 + 2x2 + x3 = 1 Xiji ng Uni v.