关晋瑞,温瑞萍
(太原师范学院 数学系,山西 晋中 030619)
考虑线性方程组的求解:
Ax=b,
(1)
其中A∈n×n是非奇异的矩阵,b∈n是给定的向量.线性方程组广泛出现在科学计算和工程应用的许多领域中,如大地测量、结构分析、网络分析、数据分析、最优化、非线性方程组以及微分方程数值解中的许多问题最终都归结为线性方程组的求解[1-3].如何快速有效地求解线性方程组是数值代数的核心问题[4-8],也是目前仍在继续研究的重要课题之一[9-12].
线性方程组的数值方法大体上可分为直接法和迭代法两大类.国内外许多学者对此进行了深入的研究,提出了众多的有效方法,如经典的Gauss消元法、平方根法、Jacobi迭代法、Gauss-Seidel迭代法、SOR迭代法、共轭梯度法和Krylov子空间方法[1,5,8],以及近年来著名的HSS迭代法[9]等.这些都是求解线性方程组的有效的数值方法.
本文研究病态线性方程组的求解问题.当线性方程组系数矩阵的条件数较大时或者系数矩阵接近奇异时,方程组比较病态,这时直接求解得到解的精度较低甚至是完全错误的.对此,本文提出了一种迭代改进方法,对原方程组的解进行迭代改进,在一定情况下可以有效提高所求得解的精度.
本节介绍文中的符号记法以及后面将要用到的一些预备知识.
求解线性方程组(1)最简单的一类迭代方法是分裂迭代法,它的基本思路是,给定系数矩阵的一个分裂:
A=M-N,
其中M是可逆的,代入式(1)得到原方程组的等价形式Mx=Nx+b,进而可以得到迭代法:
xk+1=M-1Nxk+M-1b,k=0,1,2,….
(2)
迭代法(2)称为单步线性定常迭代法,其中M-1N为迭代矩阵.
如果对于任意给定的初始向量x0∈n,由式(2)生成的序列{xk}都是适定的且收敛于方程组(1)的解x*,则称(2)是收敛的.下述定理给出了单步线性定常迭代法收敛的充要条件.
定理1[13]单步线性定常迭代法(2)收敛当且仅当其迭代矩阵M-1N的谱半径ρ(M-1N)<1.
本节提出一类迭代法以求解线性方程组(1),并给出新方法的收敛性分析.
设α>0是一个给定的参数,利用分裂A=(αI+A)-αI,可以得到迭代格式:
(αI+A)xk+1=αxk+b,k=0,1,2,….
(3)
或者:
xk+1=(αI+A)-1αxk+(αI+A)-1b,k=0,1,2,…
其中x0∈n是给定的初始向量.
迭代法(3)作为一种迭代方法似乎没有太大的应用价值,这是因为与原方程组(1)相比,迭代法(3)需要求解一系列的方程组才能得到原方程组的解.但是,同原方程组相比,迭代法(3)中系数矩阵更加对角占优,当原方程组很病态或者接近奇异时,适当选取参数可以使得式(3)的系数矩阵病态降低从而易于求解.此外,如果利用某种方法已得到原方程组的一个近似解,那么可以利用迭代(3)进行改进,得到精度更高的一个解,这是迭代法(3)的理论意义及价值.
以下分析迭代法(3)的收敛性.可以证明,当系数矩阵为正稳定矩阵时,迭代法总是无条件收敛的.
定义1[14]设A∈n×n,若A的每个特征值都具有正实部,则称A为正稳定矩阵.
定理2 设线性方程组(1)的系数矩阵A∈n×n为正稳定矩阵,则对任意的参数α>0,迭代法(3)是收敛的.
证明式(3)的迭代矩阵为Lα=(αI+A)-1α,由定理1可知迭代法(3)收敛当且仅当ρ(Lα)<1.因此,只需要验证Lα的每个特征值的模都严格小于1.
正稳定矩阵包含了很多常见的矩阵,如实对称正定矩阵、非对称正定矩阵[9]、非奇异M-矩阵[15]等,因此迭代法(3)的应用范围是较为广泛的.
容易发现,迭代法(3)的数值效果与参数α的选取密切相关.一方面,参数α越大,迭代(3)的系数矩阵越对角占优,从而越容易求解;另一方面,参数α又不能太大,否则迭代(3)会收敛很慢.因此,选取合适的参数很重要.在应用中,一般应该把参数α尽量选择小些.
另外,在迭代法(3)的每步计算中,如何有效求解方程组也是比较重要的,这时可以利用直接法去求解,如LU分解,也可以用迭代法如SOR或者共轭梯度法等去求解.
当线性方程组(1)的系数矩阵为一般的方阵时,可以考虑其对应的正则化方程组:
ATAx=ATb.
应用迭代法(3)于上述正则化方程组可得迭代:
(αI+ATA)xk+1=αxk+ATb,k=0,1,2,…
(4)
由于ATA是实对称正定的,根据定理2上述迭代是收敛的.迭代法(4)的缺陷是增大了条件数,但是扩大了迭代法(3)的应用范围.
本节中,通过两个数值例子来验证所提出的新方法的可行性和有效性,并分别给出利用Matlab中基本命令“x=A”直接求解以及利用新方法计算的结果.实验中,所有的程序都用MATLAB(R2012a)编写并在笔记本(2 G CPU和8 G内存)上运行.
例1 考虑线性方程组Ax=b,其中A为n阶Hilbert矩阵,b为n维列向量,且选取b使得方程组的精确解为x*=(1,1,…,1)T.随着n的增大,系数矩阵的条件数急剧增长,方程组变得非常病态[13].
表1 两种方法求得解的精度Tab.1 The accuracy of solution by the two methods
此外,当n=20时,对于参数α不同的值,表2中给出了迭代法(3)要达到同样精度所需要的迭代次数.从表2中可知,参数α选取太大,迭代法(3)收敛很慢,只有参数选取较小些,方法才是比较有效的.
表2 参数取不同值时迭代法(3)所需的迭代次数Tab.2 The number of iterations by iteration method (3) with different parameters
例2 考虑线性方程组Ax=b,其中:
选取b使得方程组的精确解为x*=(1,1,…,1)T.随着n的增大,系数矩阵的条件数增长很快,方程组变得非常病态[13].
表3 两种方法求得解的精度Ta.3 The accuracy of solution by the two methods
表4 迭代法(4)所得解的精度与迭代次数的关系Tab.4 The precision of the solution by iteration method (4) with the number of iterations
由上述两个数值算例可见,新方法是可行的,而且也是有效的.
本文提出了一类求解病态线性方程组的迭代法.理论分析和数值例子表明,新方法是可行的,而且在一定情况下可以求得较高精度的解.但是新方法中的参数选择比较困难,如何选择合适的参数是一个较为复杂的问题,这有待进一步的研究.