第六讲线性方程组定常迭代方法通常，求解线性方程组的方法有·直接法：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) 已经满足精度要求, 则停止计算, 否则继续按以上的方式进行修正