预处理P=(I+C)后双分裂下的SOR 迭代法收敛性*

王慧勤，雷 刚

(宝鸡文理学院数学系，陕西 宝鸡 721013)

在科学与工程计算的数学和物理学许多领域，如流体力学计算、有限元分析中的结构与非结构问题、石油地震数据处理、数值天气预报、电力系统优化设计等问题的数值解等都归结为线性代数方程组的求解，特别是大型稀疏(即方程中许多未知量的系数为0)线性方程组的求解处于核心的地位。如在近年来流行的油藏数值模拟软件中，其解法器部分设计油藏模拟方程离散化后得到的大型稀疏线性方程组的求解，这个代数方程的求解占据了80%以上的计算量，是整个问题计算的瓶颈。大量的数据和资料显示，线性方程组的求解时间在整个问题的总计算时间中占有非常大的比重，近年来许多学者提出运用预处理的方法来加快迭代法的收敛性。

1 预备知识

在用SOR迭代法求解大型线性方程组Ax=b时(其中A是n阶方阵，x,b是n维向量)，常用文[1]中提到的因子P=(I+C)来加速或改善迭代法的收敛性。也就是对方程组两边分别左乘P(称为预处理因子)，转化为

PAx=Pb

(1)

其中I为n阶单位矩阵，C为方阵，满足

a21,a31,…,an1分别是方程组系数矩阵A对应位置上的元素。通常的处理方法是将矩阵A分裂为对角矩阵、严格上三角矩阵和严格下三角矩阵三部分[2]，不妨设A=I-L-U。那么求解方程组的SOR迭代法的迭代矩阵为

T=M-1N

(2)

其中M=(I-ωL),N=[(1-ω)I+ωU]。通常用迭代矩阵T的谱半径ρ(T)<1来确定迭代法收敛，而且ρ(T)越小，收敛速度就越快。

常见的预处理形式是在因子P作用下，对方程组的系数矩阵PA做类似于系数矩阵A的分解[2-3]，令PA=I0-L0-U0，这里I0=(I-D1),L0=(L-C+L1),U0=(U+U1)分别是矩阵PA的对角线部分、严格下三角部分和严格上三角部分，其中D1,L1,U1是CU的对应的三部分分解。那么相应地SOR方法的迭代矩阵为

(3)

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

一般情况下选择在预处理因子P的目的是加快迭代法的收敛速度，许多学者通过选取不同的预处理因子来加快SOR迭代法的收敛速度，本文在预处理的基础上，首先改变原有的矩阵分裂形式，然后引入参数使得分裂形式更加一般化，寻找参数的取值范围来加快迭代法的收敛性。

引入参数α和β，利用矩阵分裂理论，将矩阵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)

容易知道当α=1,β=1时，(4)和(5)式就是常见预处理后的分裂形式(3)式。

2 定义和引理

定义1[4]设n×n阶方阵A=(aij)的阶n≥2，称A是不可约的，如果对集合W={1,2,…,n}的任意两个非空不相交的子集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,当s≥ρ(B)时，称A为M-矩阵，特别称n×n阶矩阵A为奇异M-矩阵，当s=ρ(B)时；称n×n阶矩阵A为非奇异的M-矩阵，当s>ρ(B)时，其中ρ(B)为B的谱半径。

定义3[5]所有n×n阶的实矩阵A=(aij)n×n所组成的集合为Rn×n，其中Rn×n的子集

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

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

当A∈Zn×n，并且aii>0(∀i)成立时，称矩阵A为L-矩阵。

引理1[6](Perron-Frobenius定理)如果A为n×n阶非负方阵，那么

(i)矩阵A有非负特征值等于它的谱半径ρ(A)；

(ii)矩阵A有与谱半径ρ(A)相对应的非负特征向量；

(iii)当矩阵A的任一元素增加时，谱半径ρ(A)不减。

引理2[7]设n×n阶矩阵A为非负矩阵，则

(i)如果对常数α，满足αx≤Ax对某一个非负向量x且x≠0成立，那么就有α≤ρ(A)；

(ii)如果对常数β，Ax≤βx对某一个正向量x成立，那么就有ρ(A)≤β，进一步，如果矩阵A不可约且有0≠αx≤Ax≤βx,αx≠Ax,Ax≠βx对某一个n阶非负向量x成立，那么就有α<ρ(A)<β。

3 主要结果及证明

定理1 设线性方程组的系数矩阵A是非奇异不可约M-矩阵，且0

(ii)当0<ω≤β≤1时，有ρ(TCβ)≤ρ(T)。

证明首先证明对满足题设条件的非奇异不可约M-矩阵A在ω∈(0,1)时，有

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

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

其中x为与ρ(T)对应的正特征向量,α≥1。

由于A为非奇异不可约M-矩阵，易知(4)式对应的SOR迭代矩阵是非负不可约的，且[(1-ω)I+ωU]x=ρ(T)(I-ωL)x，结合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

从而上述结论成立。

另由文[8]的引理3知，T是一个不可约的非负矩阵，再由引理1知，存在一个正向量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)时，有TCα-ρ(T)x≤0，但TCα-λx≠0，因此，由引理2可得ρ(TCα)≤ρ(T)。

在矩阵和迭代法满足相同的条件下可以证明当0<ω≤β≤1时有

(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)时，有TCβ-ρ(T)x≤0，但TCα-λx≠0，因此，由引理2可得ρ(TCβ)≤ρ(T)。

定理2 设线性方程组的系数矩阵A是非奇异不可约M-矩阵，且0

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

当迭代矩阵T的谱半径ρ(T)<1时有

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

证明根据参数1≤α，可知矩阵(I-D1)-(L-C+L1)≥(I-αD1)-(L-C+L1)，则相应的逆矩阵就有[(I-D1)-(L-C+L1)]-1≤[(I-αD1)-(L-C+L1)]-1，又由文[9-10]中的引理3知，迭代矩阵TC是非负矩阵，再结合引理1可知其谱半径是它的一个特征值，且有与之相对应的非负特征向量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

所以，当ρ(T)≥1时，上述不等式右端大于零，从而ρ(TCα)≥ρ(TC)；

另一方面，又由于

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

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

所以，当ρ(T)<1时，上述不等式右端小于零，所以ρ(TCα)<ρ(TC)。

类似地可以证明当0<ω≤β≤1时有上述的结论成立。

4 数值例子

例1 如果矩阵A的表达式如下所示：

计算可知，以A为线性方程组的系数矩阵，对不同的当ω,α时，谱半径的比较如表1。

表1 不同参数的迭代法谱半径Table 1 comparison of the spectral radius in different parameters values

5 结 语

利用参数α,β，给出预处理后SOR迭代法的两类新的分裂形式，与常见的预处理后SOR迭代法收敛速度进行比较，得到较好结果，并找到参数的最优取值，对科学计算中遇到的线性方程组的求解问题提供理论支持，同时也为今后在求解线性方程组时可以从预处理因子和分裂形式两方面选择来加速迭代法收敛性提供帮助。

参考文献：

[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] 张谋成, 黎稳. 非负矩阵论[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.