王慧勤,雷 剛
(寶雞文理學(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ù)處理的方法來加快迭代法的收斂性。
在用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)式。
定義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)<β。
定理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é)論成立。 例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 利用參數(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.4 數(shù)值例子
5 結(jié) 語