王慧勤，雷 剛

(寶雞文理學(xué)院數(shù)學(xué)系，陜西 寶雞 721013)

在科學(xué)與工程計(jì)算的數(shù)學(xué)和物理學(xué)許多領(lǐng)域，如流體力學(xué)計(jì)算、有限元分析中的結(jié)構(gòu)與非結(jié)構(gòu)問題、石油地震數(shù)據(jù)處理、數(shù)值天氣預(yù)報(bào)、電力系統(tǒng)優(yōu)化設(shè)計(jì)等問題的數(shù)值解等都?xì)w結(jié)為線性代數(shù)方程組的求解，特別是大型稀疏(即方程中許多未知量的系數(shù)為0)線性方程組的求解處于核心的地位。如在近年來流行的油藏?cái)?shù)值模擬軟件中，其解法器部分設(shè)計(jì)油藏模擬方程離散化后得到的大型稀疏線性方程組的求解，這個(gè)代數(shù)方程的求解占據(jù)了80%以上的計(jì)算量，是整個(gè)問題計(jì)算的瓶頸。大量的數(shù)據(jù)和資料顯示，線性方程組的求解時(shí)間在整個(gè)問題的總計(jì)算時(shí)間中占有非常大的比重，近年來許多學(xué)者提出運(yùn)用預(yù)處理的方法來加快迭代法的收斂性。

1 預(yù)備知識(shí)

在用SOR迭代法求解大型線性方程組Ax=b時(shí)(其中A是n階方陣，x,b是n維向量)，常用文[1]中提到的因子P=(I+C)來加速或改善迭代法的收斂性。也就是對(duì)方程組兩邊分別左乘P(稱為預(yù)處理因子)，轉(zhuǎn)化為

PAx=Pb

(1)

其中I為n階單位矩陣，C為方陣，滿足

a21,a31,…,an1分別是方程組系數(shù)矩陣A對(duì)應(yīng)位置上的元素。通常的處理方法是將矩陣A分裂為對(duì)角矩陣、嚴(yán)格上三角矩陣和嚴(yán)格下三角矩陣三部分[2]，不妨設(shè)A=I-L-U。那么求解方程組的SOR迭代法的迭代矩陣為

T=M-1N

(2)

其中M=(I-ωL),N=[(1-ω)I+ωU]。通常用迭代矩陣T的譜半徑ρ(T)<1來確定迭代法收斂，而且ρ(T)越小，收斂速度就越快。

常見的預(yù)處理形式是在因子P作用下，對(duì)方程組的系數(shù)矩陣PA做類似于系數(shù)矩陣A的分解[2-3]，令PA=I0-L0-U0，這里I0=(I-D1),L0=(L-C+L1),U0=(U+U1)分別是矩陣PA的對(duì)角線部分、嚴(yán)格下三角部分和嚴(yán)格上三角部分，其中D1,L1,U1是CU的對(duì)應(yīng)的三部分分解。那么相應(yīng)地SOR方法的迭代矩陣為

(3)

這里MC=(I-D1)-ω(L-C+L1)],NC=[(1-ω)(I-D1)+ω(U+U1)]。

一般情況下選擇在預(yù)處理因子P的目的是加快迭代法的收斂速度，許多學(xué)者通過選取不同的預(yù)處理因子來加快SOR迭代法的收斂速度，本文在預(yù)處理的基礎(chǔ)上，首先改變?cè)械木仃嚪至研问?，然后引入?yún)?shù)使得分裂形式更加一般化，尋找參數(shù)的取值范圍來加快迭代法的收斂性。

引入?yún)?shù)α和β，利用矩陣分裂理論，將矩陣PA分裂為

和

兩種形式，其中，

MCα=(I-αD1)-ω(L-C+L1),

NCα=[(I-αD1)-ω(I-D1)+ω(U+U1)],

MCβ=β(I-D1)-ω(L-C+L1),

NCβ=[(β-ω)(I-D1)+ω(U+U1)]

則在這兩類分裂下的SOR迭代法的迭代矩陣為

(4)

(5)

容易知道當(dāng)α=1,β=1時(shí)，(4)和(5)式就是常見預(yù)處理后的分裂形式(3)式。

2 定義和引理

定義1[4]設(shè)n×n階方陣A=(aij)的階n≥2，稱A是不可約的，如果對(duì)集合W={1,2,…,n}的任意兩個(gè)非空不相交的子集S和T，S+T=W，都有i和j滿足i∈S,j∈T，使aij≠0，否則稱A為可約的。

定義2[5]如果n×n階矩陣A能表示為A=sI-B，I為n階的單位矩陣，B≥0,當(dāng)s≥ρ(B)時(shí)，稱A為M-矩陣，特別稱n×n階矩陣A為奇異M-矩陣，當(dāng)s=ρ(B)時(shí)；稱n×n階矩陣A為非奇異的M-矩陣，當(dāng)s>ρ(B)時(shí)，其中ρ(B)為B的譜半徑。

定義3[5]所有n×n階的實(shí)矩陣A=(aij)n×n所組成的集合為Rn×n，其中Rn×n的子集

Zn×n={A=(aij)n×nA∈Rn×n,

aij≤0,(?i,j,i≠j)}

當(dāng)A∈Zn×n，并且aii>0(?i)成立時(shí)，稱矩陣A為L(zhǎng)-矩陣。

引理1[6](Perron-Frobenius定理)如果A為n×n階非負(fù)方陣，那么

(i)矩陣A有非負(fù)特征值等于它的譜半徑ρ(A)；

(ii)矩陣A有與譜半徑ρ(A)相對(duì)應(yīng)的非負(fù)特征向量；

(iii)當(dāng)矩陣A的任一元素增加時(shí)，譜半徑ρ(A)不減。

引理2[7]設(shè)n×n階矩陣A為非負(fù)矩陣，則

(i)如果對(duì)常數(shù)α，滿足αx≤Ax對(duì)某一個(gè)非負(fù)向量x且x≠0成立，那么就有α≤ρ(A)；

(ii)如果對(duì)常數(shù)β，Ax≤βx對(duì)某一個(gè)正向量x成立，那么就有ρ(A)≤β，進(jìn)一步，如果矩陣A不可約且有0≠αx≤Ax≤βx,αx≠Ax,Ax≠βx對(duì)某一個(gè)n階非負(fù)向量x成立，那么就有α<ρ(A)<β。

3 主要結(jié)果及證明

定理1 設(shè)線性方程組的系數(shù)矩陣A是非奇異不可約M-矩陣，且0

(ii)當(dāng)0<ω≤β≤1時(shí)，有ρ(TCβ)≤ρ(T)。

證明首先證明對(duì)滿足題設(shè)條件的非奇異不可約M-矩陣A在ω∈(0,1)時(shí)，有

(NCα-ρ(T)MCα)x=

(ρ(T)-1)[αD1+C(1-ω)+ωL1]x

其中x為與ρ(T)對(duì)應(yīng)的正特征向量,α≥1。

由于A為非奇異不可約M-矩陣，易知(4)式對(duì)應(yīng)的SOR迭代矩陣是非負(fù)不可約的，且[(1-ω)I+ωU]x=ρ(T)(I-ωL)x，結(jié)合CL=0，有

ωCUx=[ρ(T)C(1-ω)L-(1-ω)CI]x=

[ρ(T)C-(1-ω)C]x

所以

(NCα-ρ(T)MCα)x=

{[(I-αD1)-ω(I-D1)+ω(U+U1)]-

ρ(T)[(I-αD1)-ω(L-C+LL1)]}x=

[-(α-ω)D1+αρ(T)D1+ωU1-

ρ(T)ωC+ρ(T)ωL1]x=

[α(ρ(T)-1)D1+ω(D1+U1+L1)-

ρ(T)ωC+(ρ(T)-1)ωL1]x=

[α(ρ(T)-1)D1+(ρ(T)-1)C(1-ω)+

(ρ(T)-1)ωL1]x=

(ρ(T)-1)[αD1+C(1-ω)+ωL1]x

從而上述結(jié)論成立。

另由文[8]的引理3知，T是一個(gè)不可約的非負(fù)矩陣，再由引理1知，存在一個(gè)正向量x=(x1,x2,…,xn)，滿足Tx=ρ(T)x，即

[(1-ω)I+ωU]x=(I-ωL)ρ(T)x

[(I-αD1)-ω(L-C+L1)]-1(ρ(T)-1)·

[αD1+C(1-ω)+ωL1]x

所以在ρ(T)<1，ω∈(0,1)時(shí)，有TCα-ρ(T)x≤0，但TCα-λx≠0，因此，由引理2可得ρ(TCα)≤ρ(T)。

在矩陣和迭代法滿足相同的條件下可以證明當(dāng)0<ω≤β≤1時(shí)有

(NCβ-ρ(T)MCβ)x=

(ρ(T)-1)[(1-β)(I-D1)+D1+ωL1]x

類似地有

[β(I-D1)-ω(L-C+L1)]-1(ρ(T)-1)·

[(1-β)(I-D1)+D1+ωL1]x

因此在ρ(T)<1，ω∈(0,1)時(shí)，有TCβ-ρ(T)x≤0，但TCα-λx≠0，因此，由引理2可得ρ(TCβ)≤ρ(T)。

定理2 設(shè)線性方程組的系數(shù)矩陣A是非奇異不可約M-矩陣，且0

(ii)對(duì)任意0<ω≤β≤1，有ρ(TCβ)≥ρ(TC)。

當(dāng)?shù)仃嘥的譜半徑ρ(T)<1時(shí)有

(iv)對(duì)任意0<ω≤β≤1，有ρ(TCβ)<ρ(TC)。

證明根據(jù)參數(shù)1≤α，可知矩陣(I-D1)-(L-C+L1)≥(I-αD1)-(L-C+L1)，則相應(yīng)的逆矩陣就有[(I-D1)-(L-C+L1)]-1≤[(I-αD1)-(L-C+L1)]-1，又由文[9-10]中的引理3知，迭代矩陣TC是非負(fù)矩陣，再結(jié)合引理1可知其譜半徑是它的一個(gè)特征值，且有與之相對(duì)應(yīng)的非負(fù)特征向量x=(x1,x2,…,xn)T，滿足TCx=ρ(TC)x。從而

TCαx-ρ(TC)x=TCαx-TCx=

(TCαx-ρ(T)x)-(TCx-ρ(T)x) ≥

即

TCαx-ρ(TC)x≥[(I-αD1)-ω(L-C+L1)]-1·

(ρ(T)-1)(α-1)D1x

所以，當(dāng)ρ(T)≥1時(shí)，上述不等式右端大于零，從而ρ(TCα)≥ρ(TC)；

另一方面，又由于

TCαx-ρ(TC)x=TCαx-TCx=

(ρ(T)x-TCx)-(ρ(T)x-TCαx)≤

所以，當(dāng)ρ(T)<1時(shí)，上述不等式右端小于零，所以ρ(TCα)<ρ(TC)。

類似地可以證明當(dāng)0<ω≤β≤1時(shí)有上述的結(jié)論成立。

4 數(shù)值例子

例1 如果矩陣A的表達(dá)式如下所示：

計(jì)算可知，以A為線性方程組的系數(shù)矩陣，對(duì)不同的當(dāng)ω,α?xí)r，譜半徑的比較如表1。

表1 不同參數(shù)的迭代法譜半徑Table 1 comparison of the spectral radius in different parameters values

5 結(jié) 語

利用參數(shù)α,β，給出預(yù)處理后SOR迭代法的兩類新的分裂形式，與常見的預(yù)處理后SOR迭代法收斂速度進(jìn)行比較，得到較好結(jié)果，并找到參數(shù)的最優(yōu)取值，對(duì)科學(xué)計(jì)算中遇到的線性方程組的求解問題提供理論支持，同時(shí)也為今后在求解線性方程組時(shí)可以從預(yù)處理因子和分裂形式兩方面選擇來加速迭代法收斂性提供幫助。

參考文獻(xiàn)：

[1] YUN J H. A note on the modified SOR method for Z-matrices [J]. Applied Mathematics and Computation, 2007, 194: 572-576.

[2] NIKI H, HARADA K, MORIMOTO M, et al. The survey of preconditioners used for accelerating the rate of convergence in the Gauss-Seidel method [J]. Journal of Computational and Applied Mathematics, 2004, 165: 587-600.

[3] HUANG T Z, CHENG G H, CHENG X Y. Modified SOR-type iterative method for Z-matrices [J]. Applied Mathematics and Computation, 2006,175: 258-268.

[4] 張謀成, 黎穩(wěn). 非負(fù)矩陣論[M]. 廣州: 廣東高等教育出版社,1995.

[5] YONG D M. Iterative solution of large linear systems [M]. New York: Academic Press, 1971.

[6] VARGA R S. Matrix iterative analysis [M]. Heidelberg: Spring-Verlag, 2000.

[7] LIU Q B, CHEN G L, CAI J. Convergence analysis of the preconditioned Gauss-Seidel method for H-matrices [J]. Computers and Mathematics with Applications, 2008, 56: 2048-2053.

[8] YUN J H. Convergence of SSOR multisplitting method for an H-matrix [J]. Journal of Computational and Applied Mathematics, 2008, 217: 252-258.

[9] WANG X Z, HUANG T Z, FU Y D. Comparison results on preconditioned SOR-type iterative method for Z-matrices linear systems [J]. Journal of Computational and Applied Mathematics, 2007, 206: 726-732.

[10] LI W. Comparison results for solving preconditioned linear systems [J]. Journal of Computational and Applied Mathematics, 2005, 176: 319-329.