关晋瑞,宋儒瑛
(太原师范学院 数学系,山西 晋中 030619)
在科学计算和工程应用的许多领域,例如大地测量、结构分析、网络分析、数据分析、最优化、非线性方程组以及微分方程数值解等,许多问题的求解最终都归结为线性方程组,因此研究线性方程组的求解具有重要的理论意义和应用价值[1-3].
本文我们考虑M-矩阵线性方程组
Ax=b,
(1)
其中A∈Rn×n是非奇异的M-矩阵,b∈Rn是给定的向量.这种类型的线性方程组在微分方程数值解、线性互补问题、Markov链以及物理学,生物学和经济学中都有着广泛的应用[4-6].
对于一般线性方程组的求解,已有很多经典的方法,如Gauss消去法、平方根法、Jacobi迭代法、Gauss-Seidel迭代法、SOR、共轭梯度法和Krylov子空间方法等等[3,5].而对于一些特殊类型的线性方程组,国内外许多学者也进行了深入的研究,并发展出了众多有效的算法[7-15].本文我们利用M-矩阵的特点,针对M-矩阵线性方程组提出了一类迭代法.理论分析和数值实验表明,新方法是可行的,而且在一定情况下也是较为有效的.
我们用Rm×n表示实数域上全体m×n的矩阵,用Rn表示实数域上全体n维列向量,用ρ(A)表示矩阵A的谱半径.n阶单位矩阵记为I.
设A=(aij)∈Rn×n,若A的每个元素都是非负的,则称A为非负矩阵,记为A≥0.若对于任意的i≠j都有aij≤0成立,则称A为Z-矩阵.对于Z-矩阵A,若存在非负矩阵B使得A=sI-B且s≥ρ(B)成立,则称该Z-矩阵为M-矩阵.特别地,若s>ρ(B)则称A为非奇异的M-矩阵,若s=ρ(B)则称A为奇异的M-矩阵.
引理1[16]设A∈Rn×n为Z-矩阵,则A是非奇异的M-矩阵当且仅当存在正向量v>0使得Αv>0.
分裂迭代法是求解线性方程组(1)的最简单的一类迭代法.它的基本思路是,给定系数矩阵A的一个分裂
Α=Μ-Ν,
其中M是可逆的,代入方程组(1)得到迭代法
xk+1=M-1Nxk+M-1b,k=0,1,2,…
(2)
其中M-1N称为分裂迭代法(2)的迭代矩阵.
如果对于任意给定的初始向量x0∈Rn,若由(2)生成的序列{xk}都收敛于方程组(1)的解x*,则称迭代法(2)是收敛的.
引理2[5]分裂迭代法(2)收敛当且仅当迭代矩阵M-1N的谱半径ρ(M-1N)<1.
本节,我们对M-矩阵线性方程组提出一类迭代法,并给出该方法的收敛性分析.
设A=(aij)∈Rn×n为非奇异的M-矩阵,且不妨设A的对角元素都为1,否则在方程组两端左乘D-1即可,其中D为A的对角元组成的对角矩阵.基于分裂
A=I-L-U,
其中
可以得到计算方程组(1)的Jacobi迭代法
xk+1=(L+U)xk+b,k=0,1,2,…
其中x0∈Rn是任意给定的初始向量.
若记B=L+U,则有A=I-B,B≥0,于是可将Jacobi迭代法改写为
xk+1=Bxk+b,k=0,1,2,…
(3)
由于A是非奇异的M-矩阵,可知ρ(B)<1.从而根据引理2,Jacobi迭代法对于方程组(1)是收敛的.
下面我们考虑对Jacobi迭代法进行加速.利用预处理的技巧,选取预处理子
则方程组(1)等价于方程组
PAx=Pb.
(4)
若记P=I+S,其中
(5)
则我们可以得到(4)系数矩阵如下的分裂
于是代入方程组(4)可得如下的迭代法
xk+1=(B-SA)xk+Pb,k=0,1,2,…
(6)
其中x0∈Rn是任意给定的初始向量.
容易发现,由于S的特殊结构,(6)的迭代矩阵B-SA与Jacobi迭代法(3)的迭代矩阵B仅在第一行元素不同.容易验证,由矩阵B构造矩阵B-SA的运算量大约为2n2,因此迭代矩阵B-SA是容易构造的.我们称迭代法(6)为Jacobi-like迭代法.
下面讨论迭代法(6)的收敛性.
引理3设A=(aij)∈Rn×n为对角元素都为1的非奇异M-矩阵,矩阵B和S定义如上,则B-SA是非负矩阵.
β11=a12a21+a13a31+…+a1nan1≥0,
而对于k=2,…,n有
引理4设A=(aij)∈Rn×n为对角元素都为1的非奇异M-矩阵,矩阵P定义如上,则PA也是非奇异的M-矩阵.
证明 由于PA=I-(B-SA),而B-SA是非负矩阵,于是PA为Z-矩阵.
由于A=(aij)∈Rn×n为非奇异的M-矩阵,根据引理1,存在正向量v>0使得Av>0.又P=I+S,S是非负矩阵.从而
PAv=(I+S)Av=Av+SAv>0.
利用引理1,可知PA为非奇异的M-矩阵.证毕.
定理1对于M-矩阵线性方程组(1),迭代法(6)对任意给定的初始向量x0∈Rn都是收敛的.
证明 1)若方程组(1)系数矩阵A的对角元素都是1,则根据引理3和4可知
PA=I-(B-SA),B-SA≥0,
且PA为非奇异的M-矩阵.从而由M-矩阵的定义得ρ(B-SA)<1.再利用引理2可知迭代法(6)是收敛的.
2)否则,考虑方程组(1)的等价形式D-1Ax=D-1b,显然该方程组的对角元素都是1.于是根据上述分析,迭代法(6)收敛.证毕.
给出两个数值例子来检验迭代法(6)的数值效果,并与Jacobi迭代法(3)进行比较.实验中,所有的程序都用Matlab (R2012a)编写,并在个人笔记本电脑(2G CPU和8G内存)上运行.
例1考虑M-矩阵线性方程组(1),其中
直接计算可知
其中矩阵B为Jacobi迭代法(3)的迭代矩阵,B-SA为迭代法(6) 的迭代矩阵.由于
ρ(B-SA)=0.822 1<ρ(B)=0.879 2,
迭代法(6)比Jacobi迭代法(3)收敛要快.
例2考虑M-矩阵线性方程组(1),其中
表1 例2的实验结果
从表1可见,对于参数ε的不同值,(6)所需的迭代次数都比较少.
由以上两个例子可见,本文所提的方法(6)是可行的,而且与Jacobi迭代法相比,新方法的迭代次数较少,因此加速效果是较为明显的.
本文中我们提出了一类Jacobi-like迭代法以计算M-矩阵线性方程组的解,并给出了相应的收敛性分析.数值实验表明,该方法是可行的,而且在一定情况下比Jacobi迭代法所需的迭代次数要少.