劉仲云，梅聰輝

(長沙理工大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，湖南 長沙，410114)

文中考慮迭代求解大型稀疏線性方程組

Ax=b

(1)

其中，A∈Cn×n為正定矩陣；x,b∈Cn。

對于任意的矩陣，A∈Cn×n很自然地?fù)碛幸环N Hermitian/skew-Hermitian分裂(HSS)

A=H+S

(2)

這里

(3)

基于上述分裂和古典的交替方向法，BAI等[1]提出了如下求解非 Hermitian 正定線性方程組 (1) 的 HSS 迭代方法。

HSS 迭代方法：給定初始向量x(0)，計算

(4)

直到由式(4)產(chǎn)生的迭代序列{x(k)}收斂，其中，α>0是給定的常數(shù)，k=0,1,2，…。

文獻(xiàn)[1]中指出 HSS 迭代方法無條件收斂到方程組 (1) 的唯一解，并且它的收縮因子的上界只與系數(shù)矩陣 Hermitian 部分的譜有關(guān)，與 skew-Hermitian 部分的譜以及迭代涉及的矩陣的特征向量無關(guān)。

HSS 迭代方法求解方程組 (1) 時，每一步迭代都需要解2個子方程，一個系數(shù)矩陣為αI+H，另一個為αI+S。系數(shù)矩陣為αI+H的子方程可以被直接法 (如：Cholesky 分解) 有效求解，而系數(shù)矩陣為αI+S的子方程可以被迭代法 (如：Krylov 子空間方法) 求解。另外，當(dāng)方程組 (1) 的系數(shù)矩陣是 Hermitian 時，它的收縮因子的上界與共軛梯度法的收縮因子上界相同。因此，HSS 這類迭代方法在古典分裂迭代法和子空間投影法之間架起了聯(lián)系的橋梁。

為了加快 HSS 迭代方法的收斂速度，BENZI[4]提出了一種廣義的 HSS (GHSS) 迭代方法，它主要是將H分裂成兩個 Hermitian 半正定矩陣：H=U+V，其中，V是一種簡單矩陣 (如：對角陣)，然后，將V分配給 skew-Hermitian 部分，這樣就得到一個新的分裂：A=U+(S+V)，進(jìn)而得出如下 GHSS 迭代方法。

GHSS 迭代方法：給定初始向量x(0)，計算

(5)

直到由式(5)產(chǎn)生的迭代序列{x(k)}收斂，其中，α>0是給定的常數(shù)，k=0,1,2，…。

文獻(xiàn)[4]中指出 GHSS 迭代方法收縮因子的上界為

‖(αI-U)(αI+U)-1‖2‖(αI-S-V)(αI+S+V)-1‖2，

如果U和V至少有一個是正定的，那么上式中的兩個矩陣范數(shù)至少有一個嚴(yán)格小于1，另一個至多是1。而 HSS 迭代方法收縮因子的上界為

‖(αI-H)(αI+H)-1‖2‖(αI-S)(αI+S)-1‖2，

其中，‖(αI-H)(αI+H)-1‖2<1，‖(αI-S)(αI+S)-1‖2=1。雖然這還不足以說明GHSS迭代方法的收斂速度一定快于HSS迭代方法，但是文獻(xiàn)[4]用大量的數(shù)值實驗驗證了這一論斷。

對于求解大型稀疏線性方程組 (1)，HSS 迭代方法及其變型可能會破壞原系數(shù)矩陣的稀疏性，進(jìn)而增加存儲單元，使得求解所需的存儲量和計算量增大。

受GHSS迭代思想的啟發(fā)，文中考慮一種式(1)系數(shù)矩陣的新的分裂positive-definite/positive-definite分裂(PPS)：A=P1+P2，其中，P1和P2能保留原系數(shù)矩陣的稀疏結(jié)構(gòu)，且均為正定矩陣，于是便得到如下的 PPS 迭代方法。

PPS迭代方法：給定初始向量x(0)，計算

(6)

直到由式(5)產(chǎn)生的迭代序列{x(k)}收斂，其中，β>0是給定的常數(shù)，k=0,1,2，…。

顯然，該方法較好地保留了原系數(shù)矩陣的稀疏性，同時，收斂速度比 HSS 迭代方法更快。

1 收斂性分析

將式(6)方程組中的兩步合并，便可得到一步定常迭代

x(k+1)=H(β)x(k)+G(β)b

(7)

其中，

(8)

H(β)為 PPS 迭代的迭代矩陣。

為證 PPS 迭代方法的收斂性，先給出以下引理。

引理2.1[3]令

L(β)=(βI-P)(βI+P)-1，

如果P∈Cn×n是一個正定陣，那么

有了上面這個引理，接下來分析 PPS 迭代方法的收斂性。

定理2.2令A(yù)∈Cn×n是一個正定陣，分裂A=P1+P2中的P1和P2也均為正定陣，那么對任意的β>0，迭代式(6)無條件收斂到方程組 (1) 的唯一解。

證明:式(7)中的迭代矩陣H(β)與

相似。顯然，H(β)的譜半徑

根據(jù)引理2.1以及P1和P2均為正定陣，可以得出：不等式右邊的兩個矩陣范數(shù)均小于1，因此，對于任意的β>0均有ρ(H(β))<1，證畢。

下面給出β*取值的一種粗糙估計。

。

于是，由文獻(xiàn)[6]中的引理3可知：

ρ(H(β))=ρ(L1(β)L2(β))≤φ(β)

2 PPS迭代方法在對流擴散問題中的應(yīng)用

本部分討論PPS 迭代方法在對流擴散問題中的應(yīng)用[7]，并與 HSS 迭代方法進(jìn)行了比較，進(jìn)而說明 PPS 迭代方法對于求解線性方程組 (1) 的有效性。數(shù)值實驗在Matalb(R2018b) 環(huán)境中實現(xiàn)，并且初始向量選取為零向量，終止條件為‖r(k)‖2/‖r(0)‖2≤tol，其中,r(k)表示第k次迭代的殘差向量，tol設(shè)置為10-6。

在本應(yīng)用中，求解的線性方程組是來自對(0,1)×(0,1)區(qū)域上帶有狄利克雷邊界條件的二維對流擴散方程

-(uxx+uyy)+η(x,y)ux+ζ(x,y)uy=f(x,y)

(9)

的有限差分離散，其中,系數(shù)η和ζ分別表示沿x和y方向的速度分量，且均為非負(fù)常數(shù)[8]。這里考慮一個均勻的n×n的網(wǎng)格，并且采用標(biāo)準(zhǔn)的五點中心差分近似拉普拉斯算子uxx+uyy，用一階中心差分近似ux和uy，這樣得到線性方程組

Au=v

(10)

其中,u是有限維空間中的解向量且按照自然排序(u11,u21,…,un1,u12,…,unn)。

式(10)的系數(shù)矩陣A是一個分塊三對角矩陣，它的第k行分別包含了次對角、主對角和超對角塊，階數(shù)均為n，

Ak,k-1=dIn,Ak,k=tridiag(b,a,c),Ak,k+1=eIn

(11)

(12)

在第(j,k)個網(wǎng)格點，式(10)右端向量滿足vjk=h2fjk，其中，fjk=f(jh,kh)。

當(dāng)η=ζ時，即b=d，c=e方程組 (10) 的系數(shù)矩陣A可以寫成A=In?B+B?In，其中，B=tridiag(b,a/2,c)，定義f(x,y)=ex+y，于是得到下列數(shù)值結(jié)果，見表1。

表1 HSS和PPS的IT和CPUTable 1 IT and CUP of HSS and PPS

表1中，展示了HSS和PPS迭代方法求解方程組(10)所需的IT和CPU對于不同取值的η和h，數(shù)值結(jié)果表明:無論是從收斂速度還是計算復(fù)雜度的角度考慮，PPS迭代方法求解方程組(10) 的效果都比HSS迭代方法要好，其中,無論η=2還是η=10，隨著h的減小(n的增大)，HSS迭代方法求解 (10) 的迭代步數(shù)始終超過 PPS 迭代方法的兩倍，并且在求解時間上，HSS迭代方法也接近甚至超過PPS 迭代方法的兩倍。