于春肖，苑潤浩

(燕山大學(xué) 理學(xué)院，河北 秦皇島 066004)

許多工程、物理應(yīng)用問題求解的關(guān)鍵是如何高效地求解大型稀疏矩陣方程組.直接解法由于在進(jìn)行矩陣分解時常引入大量的填充元，導(dǎo)致計算過程中存儲量及計算量一般很大，而且當(dāng)方程組系數(shù)矩陣呈現(xiàn)病態(tài)時，直接解法穩(wěn)定性差，致使任何中間舍入誤差均可能引起最終計算結(jié)果的失真.基于上述考慮，近些年來迭代解法越來越受到重視，成為求解線性代數(shù)方程組的主要方法.可是，當(dāng)系數(shù)矩陣條件數(shù)很大時，迭代法存在不收斂與收斂速度慢2個潛在問題，文獻(xiàn)[1-2]中詳細(xì)論述了預(yù)處理器的構(gòu)造可使矩陣的特征值分布更集中，降低條件數(shù).

不完全分解是一種較為通用的預(yù)處理方法，對于對稱正定矩陣，可用不完全Cholesky因子分解預(yù)處理.Meijerink[3]、Van[4]及胡家贛[5]對此作了探索，指出在實際計算時，無填充的不完全分解(簡記IC(0))使用最為廣泛，但是這樣分解的精度還不足以產(chǎn)生合適的收斂率.隨后有不少學(xué)者就此進(jìn)一步提出了相關(guān)的改進(jìn)方法，如吳建平等[6]提出了帶門檻不完全Cholesky改進(jìn)方法、張永杰等[7]提出的大型稀疏線性方程組的改進(jìn)ICCG方法.

基于此，本文提出了新的不完全Cholesky改進(jìn)分解方法即在對稱逐次超松弛(SSOR)預(yù)處理分解及其改進(jìn)分解的基礎(chǔ)上所產(chǎn)生的一類新的預(yù)處理共軛梯度法(簡稱SSOR-ICCG方法).

1 理論準(zhǔn)備(ICCG法)

針對大型稀疏線性方程組大多采用迭代解法，結(jié)合適當(dāng)?shù)念A(yù)處理技術(shù)，方程組良態(tài)時能夠在較少的迭代次數(shù)內(nèi)收斂到問題的解，并達(dá)到精度要求；方程組病態(tài)時可以很好的降低條件數(shù)，從而使方程組能夠求解.不妨取線性方程組為

Ax=b，

(1)

其中A為對稱正定方陣，x,b∈Rn，且x為未知向量，b為已知向量.

對系數(shù)矩陣A作不完全Cholesky分解，即A=LLT-R,其中L為下三角矩陣，R為剩余矩陣，而M=LLT為預(yù)處理矩陣.此時就得到所謂的不完全Cholesky分解共軛梯度法(ICCG).

由M=LLT≈A，可得到式(1)的等價方程組

Fy=g，

(2)

此處F=L-1AL-T，y=LTx,b=L-1b，于是具體迭代算法如下:

1)對?x(0)∈Rn，計算r(0)=b-Ax(0)，r～(0)=L-1r(0)，p(0)=L-Tr～(0)；

2)對于k=0,1,2,…，計算

在方程組(2)中因為A對稱正定，F(xiàn)亦然，故而M=LLT越接近于A，F(xiàn)就越近似于單位陣I，F(xiàn)的條件數(shù)就越接近最小值1，此時ICCG法的收斂速度就越快.也就是說，對于預(yù)處理矩陣M分解的好壞直接影響到對應(yīng)的ICCG法收斂性的好壞.

2 一般的預(yù)處理矩陣

實際生活中采用最多的是無填充不完全Cholesky分解，簡記ILLT(0)(或IC(0))，即取L的非零結(jié)構(gòu)和A的下三角部分的非零結(jié)構(gòu)完全一樣，如下面所示:

ILLT(0)因子分解[8]

(3)

此外文獻(xiàn)[6-7]指出，雖然上述分解保持了和A一致的稀疏性，但當(dāng)矩陣A弱對角占優(yōu)時，ILLT(0)因子分解的有效性會大大降低.文獻(xiàn)[10]中給出了其修正形式，即通過適當(dāng)?shù)膶窃拚呗?，使預(yù)處理矩陣得以完善進(jìn)而得到一種新的ICCG方法.

3 新預(yù)處理ICCG法

3.1 SSOR預(yù)處理矩陣及其改進(jìn)形式

引理1[11]若矩陣A=AT∈Rn×n正定，則A=(aij)n的對角線元素均大于零，即aii>0.

引理2[11]若矩陣A∈Rn×n非奇異，則ATA為對稱正定矩陣.

上述引理保證對任意非奇異方陣都可轉(zhuǎn)換為對稱正定陣，故假設(shè)本文所討論的矩陣A對稱正定.此外，文獻(xiàn)[5]中胡家贛先生證明了當(dāng)矩陣A對稱正定時，經(jīng)SSOR迭代預(yù)處理矩陣M=LLT的預(yù)處理后，矩陣F=L-1AL-T的條件數(shù)約等于系數(shù)矩陣A的條件數(shù)的平方根，從而使ICCG方法加速收斂.因此現(xiàn)給出以下預(yù)處理矩陣：

(4)

其中LA=-tril(A,-1)為A的嚴(yán)格下三角部分反號，D=diag(diag(A))=diag(a11,a22,…,ann)且

(5)

其中L1=LA+0.5(I-D),L=I-wL1.同時，指出當(dāng)矩陣A的對角優(yōu)勢變?nèi)鯐r，該算法的有效性仍然保持.

(6)

3.2 SSORICCG算法及其改進(jìn)算法

算法1基于前文給出的ICCG算法，引進(jìn)SSOR預(yù)處理技術(shù)，就得到了本文的SSOR-ICCG算法.具體迭代過程如下：

Step1：輸入矩陣A、預(yù)處理因子w、初始值x(0)、誤差精度ε1和最大迭代次數(shù)kmax，置k=0，并計算

Step3：給定迭代終止條件‖x(k+1)-x(k)‖2<ε1以及k>kmax，如果滿足，則算法終止，轉(zhuǎn)到step4；否則轉(zhuǎn)入step2繼續(xù)計算；

Step4：輸出數(shù)值解x(k+1)和迭代次數(shù)k.

上述過程不僅沒有改變或影響ICCG算法自身的收斂性，反而經(jīng)過預(yù)處理降低了矩陣A的條件數(shù)使得算法穩(wěn)定性加強(qiáng)，收斂速度加快.

算法2若采用3)中的預(yù)處理矩陣替換算法1中的SSOR預(yù)處理矩陣，就得到了SSOR-ICCG改進(jìn)算法.由于其迭代過程和算法1相似，故只需將對應(yīng)的矩陣L替換即可，此處就不再贅述.

為了說明上述算法的收斂性，給出如下定理：

定理[8]設(shè)A為一對稱正定矩陣且A有近似分解式M=LLT≈A，取F=L-1AL-T，則求解方程組(2)時，對于任意給定的初值x(0)，預(yù)處理ICCG方法收斂且有以下誤差估計式

(7)

4 數(shù)值算例與分析

為驗證上面理論分析、說明新預(yù)處理ICCG算法的有效可行性，進(jìn)行了如下數(shù)值模擬仿真實驗，并且數(shù)值實驗結(jié)果表明新預(yù)處理ICCG法是比較合理有效的.

文中采用相同的迭代終止條件‖x(k+1)-x(k)‖2<10-6和k≤3 000，符號k表示迭代次數(shù)，‖ε*(x)‖2表示相對誤差，w表示預(yù)處理松弛因子(其中文獻(xiàn)[9]算法w=1.2，本文算法1,2中w=0.85).

算例1：所考查的大型稀疏線性方程組為單位正方形區(qū)域Ω={(x,y)|0

為更直觀的說明文中算法的優(yōu)劣性，針對數(shù)值模擬仿真實驗所得數(shù)據(jù)可從迭代次數(shù)和相對誤差2方面給出計算結(jié)果分析，如下圖1、圖2所示.

此處須說明：1) 圖1、圖2中實際網(wǎng)格內(nèi)點數(shù)n為圖示中的平方數(shù).2) 上述數(shù)值計算及繪圖都是在Intel(R) Core(TM) i3CPU型號的微機(jī)上用Matlab 2010b實現(xiàn)的，并且隨著CPU和Matlab版本不同，數(shù)據(jù)會略有變動.

圖1 不同階數(shù)下算法收斂速度比較 圖2 不同階數(shù)下算法求解精度比較 Fig.1 Algorithm convergence speed comparison Fig.2 Algorithm accuracy comparison under under different order number different order number

由圖1、圖2易得：本文算法(即新預(yù)處理ICCG法)對求解大型稀疏線性方程組是可行的，且比一般的無填充不完全I(xiàn)CCG法(即IC(0)法)有效；此外在相同的網(wǎng)格劃分下，算法1與文獻(xiàn)[9-10]中Evans和張永杰給出的有效方法相比，求解精度相差無幾具有一定的可行性；其改進(jìn)算法2相對上述算法可在較少的迭代次數(shù)下求得高精度的數(shù)值解，而且網(wǎng)格劃分越細(xì)其求解的有效性、優(yōu)越性愈明顯.

算例2：考查大型病態(tài)線性方程組，系數(shù)矩陣為經(jīng)典病態(tài)陣Hilbert矩陣，此時方程組為

Hx=b，

易得，上述方程的精確解為x*=(1,1,…,1)T，不妨取迭代初值x(0)=(0,0,…,0)T.

同樣的根據(jù)數(shù)值仿真實驗所得數(shù)據(jù)，可從迭代次數(shù)和相對誤差2方面給出計算結(jié)果分析，如圖3、圖4所示.

圖3 不同階數(shù)下算法收斂速度比較 圖4 不同階數(shù)下算法求解精度比較 Fig.3 Algorithm convergence speed comparison Fig.4 Algorithm accuracy comparison under under different order number different order number

由上面的圖示可以得出以下結(jié)論：

1)本文算法對求解大型病態(tài)線性方程組很有效，且算法2較算法1更優(yōu)越；

2)算法1較文獻(xiàn)[9-10]解法，迭代次數(shù)略多但求解精度方面卻相差無幾，具有一定的可行性；

3)算法2與上述各種算法相比，可以通過較少次數(shù)迭代求得較高精度的數(shù)值解，且矩陣階數(shù)越大，方程組越病態(tài)其求解有效性越突出.

5 結(jié)束語

文中探討了SSOR預(yù)處理分解及其改進(jìn)分解，并基于此種分解分別給出了對應(yīng)的預(yù)處理矩陣M和新預(yù)處理ICCG法，而后從理論角度說明了算法收斂性.最后，通過數(shù)值仿真實驗進(jìn)一步證實本文算法確實具有很高的可行性和實用價值.

參 考 文 獻(xiàn):

[1] MANTEUFFEL T A. An incomplete factorization technique for positive definite linesr systems[J]. Math Comp, 1980, 34:473-497.

[2] BEAUMENS ROBERT. Iterative solution methods[J]. Applied Numerical Mathematics, 2004,51(6):437-450.

[3] MEIJERINK J A , H A VAN DER VORST. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix[J]. Mathematics of Compytation, 1977, 31(137): 148-162.

[4] HENK A，VAN DER VORST, KEES DEKKER. Conjugate gradient type methods and preconditioning[J]. J Comput Appl Math, 1988, 24(1-2): 73-87.

[5] 胡家贛.線性代數(shù)方程組的迭代解法[M].北京:科學(xué)出版社,1999.

[6] 吳建平,王正華.帶門檻不完全Cholesky分解存在的問題與改進(jìn)[J].數(shù)值計算與計算機(jī)應(yīng)用,2003,3:208-210.

WU Jianping, WANG Zhenghua. Problems and improvements to the incomplete cholesky decomposition with thresholds[J]. Journal of Numerical Methods and Computer Applications, 2003, 3: 208-210.

[7] 張永杰,孫秦,李江海.大型稀疏線性方程組的改進(jìn)ICCG方法[J].計算物理,2007,24(5):582-583.

ZHANG Yongjie, SUN Qin, LI Jianghai. An improved ICCG method for large scale sparse linear equations [J]. Chinese Journal of Computational Physics, 2007, 24(5): 582-583.

[8] 《現(xiàn)代應(yīng)用數(shù)學(xué)手冊》編委會編. 計算與數(shù)值分析卷[M].北京: 清華大學(xué)出版社, 2005.

[9] EVANS D J. The analysis and application of sparse matrix algorithms in the finite element method[J]. J Inst Math Appl,1973, 9: 427-446.

[10] 張永杰,孫秦.大型稀疏線性方程組新的ICCG方法[J].數(shù)值計算與計算機(jī)應(yīng)用,2007,28(2):135-136.

ZHANG Yongjie, SUN Qin. A new ICCG method of large scale sparse linear equations[J]. Journal on Numerical Methods and Computer Applications, 2007, 28(2): 135-136.

[11] 蔡大用,白峰杉.高等數(shù)值分析[M].北京:清華大學(xué)出版社,2000.