求非对称线性方程组的GMRES和共轭残量法_

pute βij=-(Ar(j+1),Apj)/(Apj,Apj), i=j-k+1,…,j, 7a. pj+1=r(j+1)+∑ji=j-k+1βijpi.
25
吉林大学硕士学位论文
第四章 数值算例
在工程问题的数值模拟中， 无论边界元法还是有限元法 都将形成离散的方程矩阵， 二者不同点在于矩阵系数及性质 上。 有限元法形成的刚度系数矩阵是稀疏矩阵， 边界元积分 方程离散后，一般得到的影响系数矩阵是非对称的稠密矩 阵。 传统边界元法求解线性方程组， 一般采用高斯消去法或 是由高斯消去法派生的相关方法， 解方程所需的计算次数与 方程秩的三次方成比例。 因此当方程组规模变大时， 利用高 斯消去法为基础的直接法变得十分不利， 因为随着方程组规 模的扩大， 直接法所需的求解时间迅速增大， 因此对此类问 题的求解，迭代方法成为了有力的工具。而 GMRES(m)算法 是在 Arnodi 算法的基础上提出的解决具稠密系数矩阵的线 性方程组的有力的数值解法， 其计算次数明显减少， 因此使 用 GMRES(m)算法求解大型边界元法系数方程组是十分有效 的。下面以弹性问题模型为例，来简单介绍 GMRES(m)算法 在求解边界元问题中的应用。 边界元法计算弹性问题时， 首先需要确定物体的边界曲
j=0,…,m-1.
由 r(m)的表示式和 Api 的正交性，即得
α j = (r ( 0) , Ap j ) / ( Ap j , Ap j ) ，
（3-17）得证.假设 x(m-1)已知，现来确定 x（m）.由(3-17),
x ( m ) = x ( m −1) + α m−1 p m −1 ,
（3-18）
证明 近似解和相应的残向量可写成：
23
吉林大学硕士学位论文
m−1 i =0
x ( m) = x ( 0) + ∑ αi pi ，
r
(m)
=r
(0)
+ ∑α i Api .
i =0
m −1
为使||r(m)||2 取极小，必有
∂ ( m) r ∂α j
2 2
= 0 ，即必有正交关系
（r(m),Apj）=0,
Ap j +1 = Ar
( j +1)
+ ∑ βij Api .）每一步算法的计算量约比 GMRES
i =0
j
多 50%. 如同从 FOM 到 IOM 一样， 限制行 7 中求和的项数， 则可 得截断的 GCR，称为 ORTHOMIN（ ORTHOMIN（k）.这只须将算法 10 中的 行 6 和行 7 替换成
21
吉林大学硕士学位论文
6. 7.
β j :=(r
(j+1)
,r(j+1))/(r(j),r(j))
p(j+1):=r(j+1)+ β j p(j)
8. EndDo
这说明 CG 法可视为 FOM 用于对称正定方程组的特例. 同样地， 如果将 GMRES 用于 A 是 Hermite 矩阵 （ AH = A T = A ） 的特殊情况，则可得共轭残量法（CR）.此时，残量是 A-正 交的，即共轭的，而 Api(i=0,1,…)是正交的.对照 CG 算法 8 的结构，并注意到上述正交性条件，可得共轭残量法.由 于残量是共轭的，算法由此得名. 算法 9 共轭残量法（CR）
pute r(0):=b-Ax(0),p0:=r(0) 2.For j=0,1,…until convergence Do： 3. 4. 5. 6. 7. 8.
α j : = (r ( j ) , Ar ( j ) ) / ( Ap j , Ap j )
x(j+1):=x(j)+ α j p j
26
吉林大学硕士学位论文
GMRES（FGMRES）[4]. 为保证预处理方程组系数矩阵的正定性， 一般应要求预 处理矩阵 M 是（对称）正定的.
§3.4 共轭残量法和 ORTHOMIN 算法
将 Arnoldi 方法用于实对称阵 A，则由算法 1 得到的 Hessenberg 阵 Hm 是对称的三对角矩阵.事实上，由（1-4）
α j =(r ,Apj)/(Apj,Apj)
(j)
x(j+1)=x(j)+ α j pj r(j+1)=r(j)- α j Apj
β ij =-(Ar
(j+1)
,Apj)/(Apj,Apj),i=0,1,…j
j
pj+1=r(j+1)+ ∑i =0 βij pi
8.EndDo
（实际计算时，可在行 5 下计算 Ar(j+1),在行 7 下计算
r ( j +1) : = r ( j ) − α j Ap j
β j : = (r ( j +1) , Ar ( j +1) ) / (r ( j ) , Ar ( j ) )
p j +1: = r ( j +1) + β j p j
Compute Apj+1=Ar(j+1)+ β j Apj
22
（Api,Apj）=0,i≠j,
则在仿射子空间 x(0)+κm(A,r(0))上极小化残量范数的近似解 x为
x =x
（m）
(0)
(r ( 0 ) , Api ) ⋅ pi . +∑ i = 0 ( Ap i , Ap i )
m −1
(3-17)
进一步，x(m)可由 x(m-1)算出：
x ( m) = x ( m−1) + (r ( m−1) , Apm−1 ) p . ( Apm−1 , Apm−1 ) m−1
j =0 m− 2
从而， α m−1 = (r ( m−1) , Apm−1 ) / ( Apm−1 , Apm−1 ) . # 于是有下面的 算法 10 广义共轭残量法（GCR）
24
吉林大学硕士学位论文
pute r(0)=b-Ax(0),set po=r(0). 2.For j=0,1,2,…, until covergence Do: 3. 4. 5. 6. 7.
α m −1 = (r ( 0) , Ap m −1 ) /( Ap m −1 , Ap m −1 ),
注意由 r(m)的表示式，有 r(m-1)=r(0)- ∑ α j Ap j .这样，再由正交
j =0
m− 2பைடு நூலகம்
性,有
(r ( m −1) , Ap m −1 ) = (r ( 0 ) , Ap m −1 ) − ∑ α j ( Ap j , Ap m −1 ) = (r ( 0 ) , Ap m −1 ).
吉林大学硕士学位论文
9. EndDo
CR 法的计算量比 CG 法大，而两者的收敛性是类似的， 故一般讲 CG 法优于 CR 法. Krylov 子空间基底的选取，是各种方法的基础.GMRES 用的是正交基；在 CG 算法中，pj 是 A-正交的，即共轭的； 在 CR 算法中，Api 是正交的，即 pi 是 ATA-正交的.下面的结 果可将 CR 算法推广到非对称情形（此时，Hm 不再是三对角 的）. 命题 9 令 p0，p1，…，pm-1 是 Krylov 子空间κm(A,r(0)) 的一组基，它们是 ATA-正交的，即
Hm=VmTAVm,知 Hm 对称，从而它必是三对角的.此时，由相应的
FOM 算法可导出 DIOM（2）.可以证明，算法产生的残向量序 列 r(m)=b-Ax(m)是互为正交的，即（r(i),r(j)）=0,i≠j；而辅 助向量 pi 是 A-正交的，即（Api,pj）=0,i≠j.由此两性质， 可导出共轭梯度法（CG）--算法 8. 算法 8 共轭梯度法(CG)
1. Compute
r(0):=b-Ax(0),
p(0)=r(0)
2 .For j=0,1,…, until convergence Do: 3. 4. 5.
α j :=(r ,r )/(Ap ,p )
(j) (j) (j) (j)
x(j+1):=x(j)+ α j p(j) r(j+1):=r(j)- α j Ap(j)