第六讲线性方程组定常迭代方法通常，求解线性方程组的方法有·直接法：PLU分解，LDLT分解，Cholesky分解等·选代法：·定常（经典，不动点）选代法：Jacobi/Gauss-Seidel，SOR，AOR等·Krylov子空间送代法：CG,MIRES,GMRES,BiCGStab等·快速方法：·基于各种快速变换，如FFT,DCT,DST等·代数多重网格法（Algebraicmultigrid)·快速多极子算法（Fastmultipole)有些方法可能只适用于某类特殊方程，如快速方法.在实际应用中，这些方法可以结合使用，如混合（hy-brid)算法，预处理算法（preconditioning）等本讲以Poisson方程为例，主要介绍下面方法：·定常选代方法：JacobiGS,SOR，AORSSOR·Chebyshev加速·离散Poisson快速方法：快速Sine变换关于Krylov子空间迭代方法（CG，GMRES)，我们将在下一讲介绍更多选代方法可参见[7]6.1离散Poisson方程在本讲中，我们以一个典型的线性方程组为例，逐个介绍各种选代方法，并比较它们之间的性能。这个方程组就是二维Poisson方程经过五点差分离散后得到的线性方程组，6.1.1一维Poisson方程首先介绍一维Poisson方程的离散.考虑如下带Dirichlet边界条件的一维Poisson方程(_u(a) =(a),00<<1da2(6.1)u(0) = a, u(1) = b,其中f(r)是给定的函数,u()是需要计算的未知函数.171
第六讲 线性方程组定常迭代方法 通常, 求解线性方程组的方法有: • 直接法: PLU 分解, LDLT 分解, Cholesky 分解等 • 迭代法: • 定常 (经典, 不动点) 迭代法: Jacobi/Gauss-Seidel, SOR, AOR 等 • Krylov 子空间迭代法: CG, MIRES, GMRES, BiCGStab 等 • 快速方法: • 基于各种快速变换, 如 FFT, DCT, DST 等 • 代数多重网格法 (Algebraic multigrid) • 快速多极子算法 (Fast multipole) 有些方法可能只适用于某类特殊方程, 如快速方法. 在实际应用中, 这些方法可以结合使用, 如混合 (hybrid) 算法, 预处理算法 (preconditioning) 等. 本讲以 Poisson 方程为例, 主要介绍下面方法: • 定常迭代方法: Jacobi, GS, SOR, AOR, SSOR • Chebyshev 加速 • 离散 Poisson 快速方法: 快速 Sine 变换 关于 Krylov 子空间迭代方法 (CG, GMRES), 我们将在下一讲介绍. 更多迭代方法可参见 [7]. 6.1 离散 Poisson 方程 在本讲中, 我们以一个典型的线性方程组为例, 逐个介绍各种迭代方法, 并比较它们之间的性能. 这 个方程组就是二维 Poisson 方程经过五点差分离散后得到的线性方程组. 6.1.1 一维 Poisson 方程 首先介绍一维 Poisson 方程的离散. 考虑如下带 Dirichlet 边界条件的一维 Poisson 方程 − d 2u(x) dx2 = f(x), 0 < x < 1, u(0) = a, u(1) = b, (6.1) 其中 f(x) 是给定的函数, u(x) 是需要计算的未知函数. 171
第六讲线性方程组定常选代方法.172.差分离散取步长h：节点设为工=ih，i=0.1,2....，n+1.我们采用二阶中心差分来近似二阶导n+1数，可得du(r)2u()-u(i-1)-u(i+1) +0(h2.i=1.2.....ndr2h2将其代人(6.1),舍去高阶项后就可得到Poisson方程在工i点的近似离散方程-ui-1 + 2u - ui+1 = h f,其中fi=f(t),ui为ur)的近似.令i=1.2....,n，则可得n个线性方程.写成矩阵形式为Tnu=f,(6.2)其中fi+uoui2-1f22会 tridiag(-1,2,-1),:(6.3)T,fn-1Un-2-1fn + un+1系数矩阵T，的性质易知，T，是不可约弱对角占优的.另外，T，是对称矩阵，因此其特征值都是实数，事实上，我们有下面的结论.引理6.1T，的特征值和对应的特征向量分别为k元入=2-2cosn+i'k2k元2nkrk=1,2,...,nsin,sinn+in+lnt1即Tn=ZZT,其中A=diag(入1,入2,..,入n)是对角矩阵,Z=[21,22,...,zn] 是正交矩阵口证明.直接代人验证即可，由此可知，T，是对称正定的.更一般地，我们有下面的结论，引理6.2设T=tridiag(a,b,c)ERnxn，则T的特征值为kT入=b-2Vaccosk=1,2,...,n,n+i对应的特征向量为Zk，其第i个分量为jknZ:(G) =n+1
· 172 · 第六讲 线性方程组定常迭代方法 差分离散 取步长 h = 1 n + 1 , 节点设为 xi = ih, i = 0, 1, 2, . . . , n + 1. 我们采用二阶中心差分来近似二阶导 数, 可得 − d 2u(x) dx2 xi = 2u(xi) − u(xi−1) − u(xi+1) h 2 + O ( h 2 · d 4u dx4 ∞ ) , i = 1, 2, . . . , n. 将其代入 (6.1), 舍去高阶项后就可得到 Poisson 方程在 xi 点的近似离散方程 −ui−1 + 2ui − ui+1 = h 2 fi , 其中 fi = f(xi), ui 为 u(xi) 的近似. 令 i = 1, 2, . . . , n, 则可得 n 个线性方程. 写成矩阵形式为 Tnu = f, (6.2) 其中 Tn = 2 −1 −1 . . . . . . . . . . . . −1 −1 2 ≜ tridiag(−1, 2, −1), u = u1 u2 . . . un−1 un , f = f1 + u0 f2 . . . fn−1 fn + un+1 . (6.3) 系数矩阵 Tn 的性质 易知, Tn 是不可约弱对角占优的. 另外, Tn 是对称矩阵, 因此其特征值都是实数. 事实上, 我们有下 面的结论. 引理 6.1 Tn 的特征值和对应的特征向量分别为 λk = 2 − 2 cos kπ n + 1 , zk = √ 2 n + 1 · [ sin kπ n + 1 , sin 2kπ n + 1 , . . . , sin nkπ n + 1]⊺ , k = 1, 2, . . . , n, 即 Tn = ZΛZ ⊺ , 其中 Λ = diag(λ1, λ2, . . . , λn) 是对角矩阵, Z = [z1, z2, . . . , zn] 是正交矩阵. 证明. 直接代入验证即可. □ 由此可知, Tn 是对称正定的. 更一般地, 我们有下面的结论. 引理 6.2 设 T = tridiag(a, b, c) ∈ R n×n, 则 T 的特征值为 λk = b − 2 √ ac cos kπ n + 1 , k = 1, 2, . . . , n, 对应的特征向量为 zk, 其第 j 个分量为 zk(j) = (a c )j 2 sin jkπ n + 1
6.1离散Poisson方程.173 -特别地，若a=C=1,则对应的单位特征向量为2k元2k元nkr,sinsinsinZkn+1n+1n+1nt口证明，留作练习，自行验证由引理6.1可知，T，的最大特征值为nTnT4sin2~4,n+12(n + 1)最小特征值为4sin2(n+1)因此，当n很大时，Tn的谱条件数约为4(n + 1)2K2(Tn)~元2十矩阵T，有时也称为二阶差分矩阵，可以写成阶差分矩阵的乘积，即T，二DDT其中-111-1ERnx(n+1)D=-11其中D是一阶差分矩阵，需要注意的是，D不是方阵，因此不能用这个分解来求解线性方程组Tnz=b.6.1.2二维Poisson方程现在考虑二维Poisson方程u(r,y)u(r,y)=f(,y)，-Au(.y)=(r,y)E20r2Oy2(6.4)u(a,y) =uo(r,y),(r,y) e02其中2=[0,1]×[0.1]为求解区域，02表示2的边界五点差分离散为了简单起见，我们在z-方向和y-方向取相同的步长h=，节点设为=ih,i=jh，2+1i,j=0,1,2....,n.在-方向和y-方向同时采用二阶中心差分近似，可得u(r,y)2u(ri,yi)-u(ri-1,yi)-u(ri+1,yj)M0r2h2I(z4.9g)
6.1 离散 Poisson 方程 · 173 · 特别地, 若 a = c = 1, 则对应的单位特征向量为 zk = √ 2 n + 1 · [ sin kπ n + 1 , sin 2kπ n + 1 , . . . , sin nkπ n + 1]⊺ . 证明. 留作练习, 自行验证. □ 由引理 6.1 可知, Tn 的最大特征值为 2 ( 1 − cos nπ n + 1) = 4 sin2 nπ 2(n + 1) ≈ 4, 最小特征值为 2 ( 1 − cos π n + 1) = 4 sin2 π 2(n + 1) ≈ ( π n + 1)2 . 因此, 当 n 很大时, Tn 的谱条件数约为 κ2(Tn) ≈ 4(n + 1)2 π 2 . † 矩阵 Tn 有时也称为 二阶差分矩阵, 可以写成一阶差分矩阵的乘积, 即 Tn = DD⊺ , 其中 D = −1 1 −1 1 . . . . . . −1 1 ∈ R n×(n+1) . 其中 D 是 一阶差分矩阵. 需要注意的是, D 不是方阵, 因此不能用这个分解来求解线性方程组 Tnx = b. 6.1.2 二维 Poisson 方程 现在考虑二维 Poisson 方程 −∆u(x, y) = − ∂ 2u(x, y) ∂x2 − ∂ 2u(x, y) ∂y2 = f(x, y), (x, y) ∈ Ω, u(x, y) = u0(x, y), (x, y) ∈ ∂Ω, (6.4) 其中 Ω = [0, 1] × [0, 1] 为求解区域, ∂Ω 表示 Ω 的边界. 五点差分离散 为了简单起见, 我们在 x-方向和 y-方向取相同的步长 h = 1 n + 1 , 节点设为 xi = ih, yj = jh, i, j = 0, 1, 2, . . . , n. 在 x-方向和 y-方向同时采用二阶中心差分近似, 可得 ∂ 2u(x, y) ∂x2 (xi,yj ) ≈ 2u(xi , yj ) − u(xi−1, yj ) − u(xi+1, yj ) h 2
.174.第六讲线性方程组定常选代方法82u(a,)2u(i,)-u(iyj-1)-(iyj+1)h2Oy?(zt,u)代人（6.4)后，就可以得到二维Poisson方程在（iyi）点的近似离散方程4ui.j-ui-1,j -ui+1,j - uij-1- uij+1= h’ fij其中fü=f(ti，yi)，uis为u(i,y）的近似.这个离散格式可以用下面的图来描述(ri, yj+1)-1(ri-1,yi)(ti,y)(ri+1,y)-141(ri,yj-1)-1写成矩阵形式即为Tu=hf.(6.5)其中TITn+TnI,u=[u1,..,un.1,ui,2,..,un,2,..,u,.n....,un,n].这里表示Kronecker乘积，T，为一维Poisson方程离散后的系数矩阵，即（6.3)十需要注意的是，系数矩阵与网格点的排序有关，不同的排序方式对应不同的系数矩阵。这里是按自然顺序排列的十相类似地，如果对三维poisson方程进行中心差分离散，则对应的系数矩阵为TII+I&T.&I+II&Tn在后面介绍方法时，我们都以二维离散Poisson方程（6.5）为例系数矩阵T的性质由于T=ITn+T,I，根据Kronecker乘积的性质，可以立即得到下面的结论定理6.1设Tn=ZAZT,其中Z=[21,22,..,zn]为正交阵,A=diag(入1,2,..,入n为对角阵，则T的特征值分解为T=(ZZ)(IA+AI)(ZZ)T,即T的特征值为入+入j,对应的特征向量为zi③z(i,j=1,2,.-,n)由于T对称正定，其条件数为nTnTsin21COS4(n + 1)2Amar(T)2(n +1)n+ 1r(T)T元2Amin(T)sin21= cosn+12(n + 1)
· 174 · 第六讲 线性方程组定常迭代方法 ∂ 2u(x, y) ∂y2 (xi,yj ) ≈ 2u(xi , yj ) − u(xi , yj−1) − u(xi , yj+1) h 2 . 代入 (6.4) 后, 就可以得到二维 Poisson 方程在 (xi , yj ) 点的近似离散方程 4ui,j − ui−1,j − ui+1,j − ui,j−1 − ui,j+1 = h 2 fi,j , 其中 fij = f(xi , yj ), ui,j 为 u(xi , yj ) 的近似. 这个离散格式可以用下面的图来描述. (xi−1, yj ) −1 (xi+1, yj ) −1 (xi , yj ) 4 (xi , yj−1) −1 (xi , yj+1) −1 (xi , yj ) 写成矩阵形式即为 T u = h 2 f, (6.5) 其中 T ≜ I ⊗ Tn + Tn ⊗ I, u = [u1,1, . . . , un,1, u1,2, . . . , un,2, . . . , u1,n, . . . , un,n]. 这里 ⊗ 表示 Kronecker 乘积, Tn 为一维 Poisson 方程离散后的系数矩阵, 即 (6.3). † 需要注意的是, 系数矩阵与网格点的排序有关, 不同的排序方式对应不同的系数矩阵. 这里是按 自然顺序排列的. † 相类似地, 如果对三维 poisson 方程进行中心差分离散, 则对应的系数矩阵为 Tn ⊗ I ⊗ I + I ⊗ Tn ⊗ I + I ⊗ I ⊗ Tn. 在后面介绍方法时, 我们都以二维离散 Poisson 方程 (6.5) 为例. 系数矩阵 T 的性质 由于 T = I ⊗ Tn + Tn ⊗ I, 根据 Kronecker 乘积的性质, 可以立即得到下面的结论. 定理 6.1 设 Tn = ZΛZ ⊺ , 其中 Z = [z1, z2, . . . , zn] 为正交阵, Λ = diag(λ1, λ2, . . . , λn) 为对角阵, 则 T 的特征值分解为 T = (Z ⊗ Z)(I ⊗ Λ + Λ ⊗ I)(Z ⊗ Z) ⊺ , 即 T 的特征值为 λi + λj , 对应的特征向量为 zi ⊗ zj (i, j = 1, 2, . . . , n). 由于 T 对称正定, 其条件数为 κ(T) = λmax(T) λmin(T) = 1 − cos nπ n + 1 1 − cos π n + 1 = sin2 nπ 2(n + 1) sin2 π 2(n + 1) ≈ 4(n + 1)2 π 2
6.2定常选代方法.175.故当n越来越大时，K(T)→8o，即T越来越病态求解二维离散Poisson方程的常用方法假定网格部分为n×n,并记N=n?方法串行时间存储空间直接法稠密Cholesky分解O(N3)O(N2)显式求逆O(N2)O(N2)O(N3/2)O(N2)带状Cholesky分解O(N3/2)稀疏Cholesky分解O(N log N)JacobiO(N2)O(N)基本选代法0(N2)Gauss-SeidelO(N)O(N3/2)SORO(N)O(N5/4)O(N)带Chebyshev加速的SSORKrylov子空间迭代CG (共轭梯度法)O(N3/2)O(N)O(N5/4)CG(带修正IC预处理)O(N)快速方法DST（快速Sine变换）O(N log N)O(N)块循环约化O(N)O(N log N)MultigridO(N)O(N)6.2定常选代方法随着矩阵规模的增大，直接法的运算量也随之快速增长，对于大规模的线性方程组，由于运算量太大，直接法一般不再被采用，取而代之的是选代方法当直接求解At=b比较困难时，我们可以求解一个近似等价方程组Mr=b,其中M可以看作是A在某种意义下的近似.设Mz=b的解为(1)易知它与原方程的解工，=A-1b之间的误差满足A(r+- r(1)) = b - Ar(1),如果（1）已经满足精度要求，即非常接近真解工，则可以停止计算，否则需要修正.设修正量为，则△r满足方程A△r=b-Ar(1).但由于直接求解该方程比较困难，因此我们还是通过求解近似方程MA=b-A(1),得到一个近似的修正量△工.于是修正后的近似解为r(2) = r(1) + = (1) + M-1(b - A(1)如果(2）已经满足精度要求，则停止计算，否则继续按以上的方式进行修正
6.2 定常迭代方法 · 175 · 故当 n 越来越大时, κ(T) → ∞, 即 T 越来越病态. 求解二维离散 Poisson 方程的常用方法 假定网格剖分为 n × n, 并记 N = n 2 . 方法 串行时间 存储空间 直接法 稠密 Cholesky 分解 O(N3 ) O(N2 ) 显式求逆 O(N2 ) O(N2 ) 带状 Cholesky 分解 O(N2 ) O(N3/2) 稀疏 Cholesky 分解 O(N3/2) O(N log N) 基本迭代法 Jacobi O(N2 ) O(N) Gauss-Seidel O(N2 ) O(N) SOR O(N3/2) O(N) 带 Chebyshev 加速的 SSOR O(N5/4) O(N) Krylov 子空间迭代 CG (共轭梯度法) O(N3/2) O(N) CG (带修正 IC 预处理) O(N5/4) O(N) 快速方法 DST (快速 Sine 变换) O(N log N) O(N) 块循环约化 O(N log N) O(N) Multigrid O(N) O(N) 6.2 定常迭代方法 随着矩阵规模的增大, 直接法的运算量也随之快速增长. 对于大规模的线性方程组, 由于运算量太 大, 直接法一般不再被采用, 取而代之的是迭代方法. 当直接求解 Ax = b 比较困难时, 我们可以求解一个近似等价方程组 Mx = b, 其中 M 可以看作是 A 在某种意义下的近似. 设 Mx = b 的解为 x (1) . 易知它与原方程的解 x∗ = A−1 b 之间的误差满足 A(x∗ − x (1)) = b − Ax(1) . 如果 x (1) 已经满足精度要求, 即非常接近真解 x∗, 则可以停止计算, 否则需要修正. 设修正量为 ∆x, 则 ∆x 满足方程 A∆x = b − Ax(1) . 但由于直接求解该方程比较困难, 因此我们还是通过求解近似方程 M∆x = b − Ax(1) , 得到一个近似的修正量 ∆˜ x. 于是修正后的近似解为 x (2) = x (1) + ∆˜ x = x (1) + M−1 (b − Ax(1)). 如果 x (2) 已经满足精度要求, 则停止计算, 否则继续按以上的方式进行修正
.176.第六讲线性方程组定常选代方法不断重复以上步骤，于是，我们就得到一个序列(1),(2满足以下递推关系(k+1) = 2(k) + M-1(b - Ar(k)), k = 1, 2, ...这就构成了一个选代方法，由于每次选代的格式是一样的，因此称为定常选代通常，构造一个好的定常选代，需要考虑以下两点：(1）以M为系数矩阵的线性方程组必须要比原线性方程组更容易求解；(2）M应该是A的一个很好的近似，且送代序列【]要收敛常用的定常选代方法是基于矩阵分裂的迭代方法，主要包括：·Jacobi方法·Gauss-Seidel (G-S)方法·超松弛（SOR，SuccessiveOver-Relaxation）方法·对称超松弛（SSOR，SymmetricSOR）方法·加速超松弛（AOR，AcceleratedOver-Relaxation）方法6.2.1矩阵分裂送代方法首先给出矩阵分裂的定义定义6.1（矩阵分裂MatrixSplitting）设AERnxn非奇异，称(6.6)A-M-N为A的一个矩阵分裂,其中 M 非奇异考虑线性方程组Ar =b,(6.7)其中AERnxn非奇异.选代方法的基本思想：给定一个选代初始值r(O)，通过一定的选代格式生成一个选代序列(())=0使得lim 2(K) = ,会 A-1b.给定一个矩阵分裂（6.6)，则原方程组（6.7)就等价于Mr=N+b.于是我们就可以构造出以下的选代格式Mr(+1)= Nr(k) +b, k=0,1,...或(k+1) = M-1N(k) + M-1b ≤ Ga(k) + g, k = 0, 1,..,(6.8)其中G=M-1N称为送代矩阵.这就是基于矩阵分裂（6.6）的选代方法.易知，选取不同的M，就可以构造出不同的选代方法
· 176 · 第六讲 线性方程组定常迭代方法 不断重复以上步骤, 于是, 我们就得到一个序列 x (1), x(2), . . . , , x(k) , . . . . 满足以下递推关系 x (k+1) = x (k) + M−1 (b − Ax(k) ), k = 1, 2, . . . . 这就构成了一个迭代方法. 由于每次迭代的格式是一样的, 因此称为 定常迭代. 通常, 构造一个好的定常迭代, 需要考虑以下两点: (1) 以 M 为系数矩阵的线性方程组必须要比原线性方程组更容易求解; (2) M 应该是 A 的一个很好的近似, 且迭代序列 {xk} 要收敛. 常用的定常迭代方法是基于矩阵分裂的迭代方法, 主要包括: • Jacobi 方法 • Gauss-Seidel (G-S) 方法 • 超松弛 (SOR, Successive Over-Relaxation) 方法 • 对称超松弛 (SSOR, Symmetric SOR) 方法 • 加速超松弛 (AOR, Accelerated Over-Relaxation) 方法 6.2.1 矩阵分裂迭代方法 首先给出 矩阵分裂 的定义. 定义 6.1 (矩阵分裂 Matrix Splitting) 设 A ∈ R n×n 非奇异, 称 A = M − N (6.6) 为 A 的一个矩阵分裂, 其中 M 非奇异. 考虑线性方程组 Ax = b, (6.7) 其中 A ∈ R n×n 非奇异. 迭代方法的基本思想: 给定一个迭代初始值 x (0) , 通过一定的迭代格式生成一 个迭代序列 {x (k)}∞ k=0, 使得 lim k→∞ x (k) = x∗ ≜ A −1 b. 给定一个矩阵分裂 (6.6), 则原方程组 (6.7) 就等价于 Mx = Nx + b. 于是我们就可以构造出以下的迭代 格式 Mx(k+1) = Nx(k) + b, k = 0, 1, . . . , 或 x (k+1) = M−1Nx(k) + M−1 b ≜ Gx(k) + g, k = 0, 1, . . . , (6.8) 其中 G = M−1N 称为 迭代矩阵. 这就是基于矩阵分裂 (6.6) 的迭代方法. 易知, 选取不同的 M, 就可以 构造出不同的迭代方法
6.2定常选代方法.177.6.2.2Jacobi送代将矩阵A分裂为A=D-L-U.其中D为A的对角线部分，-L和-U分别为A的严格下三角和严格上三角部分在矩阵分裂A=M-N中取M=D,N=L+U,则可得Jacobi送代方法(k+1) = D-1(L + U)(k) + D-1b, k = 0, 1,2, ....(6.9)对应的迭代矩阵为G) = D-1(L + U).写成分量形式即为T1(k)(k+1)i=1,2,...,n.aijLiaiij=l.j+i+由于Jacobi 选代中α(k+1)的更新顺序与i无关,即可以按顺序i=1,2,,n计算,也可以按顺序i=n,n.一1,·.,2.1计算，或者乱序计算.因此Jacobi选代非常适合并行计算算法6.1.Jacobi选代1; Given an initial guess r(0)2:while not converge dofor i=l to n do3:_(k+1)4:Ciend for5:6: end while我们有时也将Jacobi迭代格式写为(k+1) = 2(k) + D-1(b - Ar(k) = 2(k) + D-1rk, k = 0, 1, 2, .*-其中rk会b一Az(）是次选代后的残量.这表明，(k+1)是通过对(k）做一个修正得到的下面给出求解二维离散Poisson方程（6.5）的Jacobi选代方法算法6.2.求解二维离散 Poisson方程的 Jacobi 选代方法I: Given an initial guess w(0)2:while not convergedofor i=1 to N do3:for j=l to N do4:+)=(f +1 +1++++1) /45:ui.s6:end for
6.2 定常迭代方法 · 177 · 6.2.2 Jacobi 迭代 将矩阵 A 分裂为 A = D − L − U, 其中 D 为 A 的对角线部分, −L 和 −U 分别为 A 的严格下三角和严格上三角部分. 在矩阵分裂 A = M − N 中取 M = D, N = L + U, 则可得 Jacobi 迭代 方法: x (k+1) = D−1 (L + U)x (k) + D−1 b, k = 0, 1, 2, . . . . (6.9) 对应的迭代矩阵为 GJ = D−1 (L + U). 写成分量形式即为 x (k+1) i = 1 aii bi − ∑n j=1,j̸=i aijx (k) j , i = 1, 2, . . . , n. † 由于 Jacobi 迭代中 x (k+1) i 的更新顺序与 i 无关, 即可以按顺序 i = 1, 2, . . . , n 计算, 也可以按顺序 i = n, n − 1, . . . , 2, 1 计算, 或者乱序计算. 因此 Jacobi 迭代非常适合并行计算. 算法 6.1. Jacobi 迭代 1: Given an initial guess x (0) 2: while not converge do 3: for i = 1 to n do 4: x (k+1) i = ( bi − ∑n j=1,j̸=i aijx (k) j ) / aii 5: end for 6: end while 我们有时也将 Jacobi 迭代格式写为 x (k+1) = x (k) + D−1 (b − Ax(k) ) = x (k) + D−1 rk, k = 0, 1, 2, . . . , 其中 rk ≜ b − Ax(k) 是 k 次迭代后的残量. 这表明, x (k+1) 是通过对 x (k) 做一个修正得到的. 下面给出求解二维离散 Poisson 方程 (6.5) 的 Jacobi 迭代方法. 算法 6.2. 求解二维离散 Poisson 方程的 Jacobi 迭代方法 1: Given an initial guess v (0) 2: while not converge do 3: for i = 1 to N do 4: for j = 1 to N do 5: u (k+1) i,j = ( h 2fi,j + u (k) i+1,j + u (k) i−1,j + u (k) i,j+1 + u (k) i,j−1 ) /4 6: end for
第六讲线性方程组定常选代方法.178.end for7:8:end while6.2.3Gauss-Seidel选代在分裂A=M-N中取M=D-L,N=U,即可得Gauss-Seidel(G-S)送代方法(k+1) = (D - L)-1U(K) + (D - L)-1b)(6.10)对应的选代矩阵为GGs = (D - L)-1U将G-S选代改写为Dr(h+1) = Lr(k+1) + Ur() + b,即可得分量形式i-1n1(k+1)(k+1)(kV7=1.2..b:aijrQijs.naiij=1j=i+1算法6.3.Gauss-Seidel选代1: Given an initial guess r(0)2:while not converge dofori=ltondo3:A(k+1)(k+1)4:Liai=+endfor5:6:end while与Jacobi选代方法类似，我们也可以给出求解问题（6.5)的G-S选代方法G-S方法的主要优点是充分利用了已经获得的最新数据但在G-S方法中，未知量的更新必须按自然顺序进行，因此不适合并行计算下面我们介绍一种适合并行计算的更新顺序：红黑排序，即将二维网格点依次做红黑记号，如右图所示在计算过程中，对未知量的值进行更新时，我们可以先更新红色节点，此时所使用的只是黑色节点的数据，然后再更新黑色节点，这时使用的是红色节点的数据.于是我们得到下面红黑排序G-S选代方法由于在更新红点时，各个点之间是相互独立的，因此可以并行计算，同样，在更新黑点时，各个点之间也是相互独立的，因此也可以并行计算
· 178 · 第六讲 线性方程组定常迭代方法 7: end for 8: end while 6.2.3 Gauss-Seidel 迭代 在分裂 A = M − N 中取 M = D − L, N = U, 即可得 Gauss-Seidel (G-S) 迭代 方法: x (k+1) = (D − L) −1Ux(k) + (D − L) −1 b. (6.10) 对应的迭代矩阵为 GGS = (D − L) −1U. 将 G-S 迭代改写为 Dx(k+1) = Lx(k+1) + Ux(k) + b, 即可得分量形式 x (k+1) i = 1 aii bi − ∑ i−1 j=1 aijx (k+1) j − ∑n j=i+1 aijx (k) j , i = 1, 2, . . . , n. 算法 6.3. Gauss-Seidel 迭代 1: Given an initial guess x (0) 2: while not converge do 3: for i = 1 to n do 4: x (k+1) i = 1 aii ( bi − i ∑−1 j=1 aijx (k+1) j − ∑n j=i+1 aijx (k) j ) 5: end for 6: end while 与 Jacobi 迭代方法类似, 我们也可以给出求解问题 (6.5) 的 G-S 迭代方法. G-S 方法的主要优点是充分利用了已经获得的最新数据. 但在 G-S 方法中, 未知量的更新必须按自然顺序进行, 因此不适合并行计算. 下面我们介绍一种适合并行计算的更新顺序: 红黑排序, 即将 二维网格点依次做红黑记号, 如右图所示. 在计算过程中, 对未知量的值进行更新时, 我们可以先更新红色 节点, 此时所使用的只是黑色节点的数据, 然后再更新黑色节点, 这 时使用的是红色节点的数据. 于是我们得到下面红黑排序 G-S 迭代 方法. 由于在更新红点时, 各个点之间是相互独立的, 因此可以并行 计算. 同样, 在更新黑点时, 各个点之间也是相互独立的, 因此也可 以并行计算
6.2定常选代方法.179 .算法6.4.求解二维离散Poisson方程的红黑排序G-S选代方法1: Given an initial guess w(0)2:while not converge dofor（i，j）为红色节点do3:(k+1)(k)()(k)(h2fi + u+1, + u-1, + u3+1 + u3-1)4:ui.j1end for5:for(i,i）为黑色节点do6:ug+1)(k+1)+ u(*+1)(k+1) + u(k+1)(h2fij +uif)+u+7:-114end for8:9:end while6.2.4SOR送代在G-S方法的基础上，我们可以通过引入一个松弛参数w来加快收敛速度.这就是SOR(SuccessiveOverrelaxation)方法[78].该方法的基本思想是将G-S方法中的第k+1步近似解与第k步近似解做一个加权平均，从而给出一个新的近似解，即(k+1) = (1 - w)() + w (D-1(L(+1) + U(k) + D-1b)(6.11)整理后即为(k+1) = (D - wL)-1 (1 - w)D + wU) r(k) + w(D - wL)-lb,(6.12)其中w称为松弛参数（relaxationparameter)，当w=1时，SOR即为G-S方法，当w1时，称为超松弛（overrelaxation方法.在大多数情况下，当w>1时会取得比较好的收敛效果十SOR方法曾经在很长一段时间内是科学计算中求解线性方程组的首选方法SOR的选代矩阵为GsOR = (D - wL)-1 (1 - w)D +wU),对应的矩阵分裂为1-wp+U.M=N-D-Lww由（6.12)可知SOR迭代的分量形式为i-1nUr(++1) = (1 - w)2(k)(k+1)KbC;aiaiij=1i=i+i-12wr(k)(k+1)kVb,-aiiaijraiij=1j=i
6.2 定常迭代方法 · 179 · 算法 6.4. 求解二维离散 Poisson 方程的红黑排序 G-S 迭代方法 1: Given an initial guess v (0) 2: while not converge do 3: for (i, j) 为红色节点 do 4: u (k+1) i,j = 1 4 ( h 2fi,j + u (k) i+1,j + u (k) i−1,j + u (k) i,j+1 + u (k) i,j−1 ) 5: end for 6: for (i, j) 为黑色节点 do 7: u (k+1) i,j = 1 4 ( h 2fi,j + u (k+1) i+1,j + u (k+1) i−1,j + u (k+1) i,j+1 + u (k+1) i,j−1 ) 8: end for 9: end while 6.2.4 SOR 迭代 在 G-S 方法的基础上, 我们可以通过引入一个松弛参数 ω 来加快收敛速度. 这就是 SOR (Successive Overrelaxation) 方法 [78]. 该方法的基本思想是将 G-S 方法中的第 k + 1 步近似解与第 k 步近似解做一 个加权平均, 从而给出一个新的近似解, 即 x (k+1) = (1 − ω)x (k) + ω ( D−1 (Lx(k+1) + Ux(k) ) + D−1 b ) . (6.11) 整理后即为 x (k+1) = (D − ωL) −1 ((1 − ω)D + ωU) x (k) + ω(D − ωL) −1 b, (6.12) 其中 ω 称为松弛参数 (relaxation parameter) . 当 ω = 1 时, SOR 即为 G-S 方法, 当 ω 1 时, 称为超松弛 (over relaxation) 方法. 在大多数情况下, 当 ω > 1 时会取 得比较好的收敛效果. † SOR 方法曾经在很长一段时间内是科学计算中求解线性方程组的首选方法. SOR 的迭代矩阵为 GSOR = (D − ωL) −1 ((1 − ω)D + ωU), 对应的矩阵分裂为 M = 1 ω D − L, N = 1 − ω ω D + U. 由 (6.12) 可知 SOR 迭代的分量形式为 x (k+1) i = (1 − ω)x (k) i + ω aii bi − ∑ i−1 j=1 aijx (k+1) j − ∑n j=i+1 aijx (k) j = x (k) i + ω aii bi − ∑ i−1 j=1 aijx (k+1) j − ∑n j=i aijx (k) j
第六讲线性方程组定常选代方法.180.算法6.5.求解线性方程组的SOR选代方法1: Given an initial guess r(0) and parameter w2:while not converge dofor i=l to n do3:2(*+1) = (1 - w) (k) +4:aiend for5:6:end while+SOR方法最大的优点是引人了松弛参数w:通过选取适当的w就可以大大提高方法的收敛速度但是SOR方法最大的难点就是如何选取最优的参数6.2.5SSOR送代方法将SOR方法中的L和U相互交换位置，则可得选代格式(k+1) = (D - wU)-1 (1 - w)D + wL) r(k) + w(D - wU)-1b将这个选代格式与SOR相结合，就可以得到下面的两步选代方法(++) = (D -wL)-[(1 -w)D + wUjr(k) + w(D -wL)-1b,(+1) = (D - wU)-[(1 - w)D + wL) r(k+) + w(D - wU)-1b.这就是对称超松弛（SSOR）迭代方法，相当于将L与U同等看待，交替做两次SOR代消去中间选代向量(+),可得(k+1) = GssORz(k) + g,其中送代矩阵GssOR = (D - wU)-1[(1 - w)D +wL](D -wL)-1[(1 -w)D +wU]对应的矩阵分裂为1[D-w(L+U)+w?LD-1]M :w(2-w(D-wL)D-1(D-wU),7w(2-w)u(2-[(1-w)D+wl)D-[(1-w)D+U],N=十对于某些特殊问题，SOR方法不收敛，但仍然可能构造出收敛的SSOR方法
· 180 · 第六讲 线性方程组定常迭代方法 算法 6.5. 求解线性方程组的 SOR 迭代方法 1: Given an initial guess x (0) and parameter ω 2: while not converge do 3: for i = 1 to n do 4: x (k+1) i = (1 − ω)x (k) i + ω aii ( bi − i ∑−1 j=1 aijx (k+1) j − ∑n j=i+1 aijx (k) j ) 5: end for 6: end while † SOR 方法最大的优点是引入了松弛参数 ω: 通过选取适当的 ω 就可以大大提高方法的收敛速度. 但是 SOR 方法最大的难点就是如何选取最优的参数. 6.2.5 SSOR 迭代方法 将 SOR 方法中的 L 和 U 相互交换位置, 则可得迭代格式 x (k+1) = (D − ωU) −1 ((1 − ω)D + ωL) x (k) + ω(D − ωU) −1 b. 将这个迭代格式与 SOR 相结合, 就可以得到下面的两步迭代方法 x (k+ 1 2 ) = (D − ωL) −1 [ (1 − ω)D + ωU] x (k) + ω(D − ωL) −1 b, x (k+1) = (D − ωU) −1 [ (1 − ω)D + ωL] x (k+ 1 2 ) + ω(D − ωU) −1 b. 这就是 对称超松弛 (SSOR ) 迭代方法, 相当于将 L 与 U 同等看待, 交替做两次 SOR 迭代. 消去中间迭代向量 x (k+ 1 2 ) , 可得 x (k+1) = GSSORx (k) + g, 其中迭代矩阵 GSSOR = (D − ωU) −1 [ (1 − ω)D + ωL] (D − ωL) −1 [ (1 − ω)D + ωU] . 对应的矩阵分裂为 M = 1 ω(2 − ω) [ D − ω(L + U) + ω 2LD−1U ] = 1 ω(2 − ω) (D − ωL)D−1 (D − ωU), N = 1 ω(2 − ω) [ (1 − ω)D + ωL] D−1 [ (1 − ω)D + ωU] . † 对于某些特殊问题, SOR 方法不收敛, 但仍然可能构造出收敛的 SSOR 方法