线性方程组的解法.doc

线性方程组的解法
1 引言
在科学研究和大型工程设计中出现了越来越多的数学问题，而这些问题往往需要求数值解。在进行数值求解时，经离散后，常常归结为求解形如Ax= b的大型线性方程组。而如插值公式，拟合公式等的建立，微分方程差分格式的构造等,均可归结为求解线性方程组的问题．在工程技术的科学计算中，线性方程组的求解也是最基本的工作之一．因此,线性 方程组的解法一直是科学和工程计算中研究最为普遍的问题,它在数值分析中占有极其重要的地位。20世纪50年代至70年 代,由于电子计算机的发展，人们开始考虑和研究在计算机上用迭代法求线性方程组Ax =b的近似解,用某种极限过程去逐渐逼近精确解，并发展了许多非常有效的迭代方法,迭代法具有需要计算机存储单元少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点。例如Jacobi方法、Gauss—Seidel 方法、SOR方法、SSOR方法 ,这几种迭代方法是最常用的一阶线性定常迭代法。
2 主要算法
20世纪50年代至70年代，人们开始考虑和研究用迭代法求解线性方程组。
Ax = b (1)
的近似解，发展了许多有效的方法，其中有Jacobi方法、Gauss—Seidel方法,SOR方法、SSOR方法，这几种迭代方法均属一阶线性定常迭代法，即若系数矩阵A的一个分裂：A =M—N ；M 为可逆矩阵,线性方程组（1）化为 ：
(M－N）X =b；
→M X = NX + b；
→X= M —1NX+ M-1b
得到迭代方法的一般公式:
X（k+1）=HX（k)+d （2）
其中：H =MN-1，d=M-1b，对任意初始向量X（0)
一阶定常迭代法收敛的充分必要条件是： 迭代矩H的谱半径小于1，即ρ（H) < 1；又因为对于任何矩阵范数恒有ρ(H)≤‖H‖，故又可得到收敛的一个充分条件为：‖H‖〈 1。
2。1 Jacobi迭代法
若D为A的对角素构成的对角矩阵，且对角线元素全不为零。系数矩阵A的一个分解:A = D－(L +U）； 这里D为A的对角素构成的对 角矩阵， 为严格下三角阵，U为严格上三角阵.
Jacobi迭代的矩阵形式为 ：
X（k+1）=HJX（k）+dJ (3)
(3) 式 中： dJ= D—1b； HJ=I－D—1A, 称HJ为Jacobi迭代矩阵。 其计算公式为：
迭代矩阵HJ的谱半径ρ（H) < 1，则对于任意迭代初值X（0),Jacobi迭代法收敛。
2。2 Gauss-Seidei迭代法
对于非奇异方程组，若D为A的对角素构成的对角矩阵,且对角线元素全不为零;系数矩阵A的一个分解:
A = （D－L)－U （5）
Gauss—Seide迭代矩阵形式为： X（k+1）=HGSX（k）+dGS
上式中： HGS=(D－L）—1U;dGS=(D－L)-1L
称HGS为G—S迭代矩阵.
计算公式为：

i=1,2，3,…n k=0,1，2，… (6）
若A为严格或不可约对角占优矩阵，或A为对称正定阵，则对于任意初值X（0），Gauss—Seide迭代法收敛.
2.3 SOR（successive over relaxation)迭代法
对于非奇异方程组,若D为A的对角素构成的 对角矩阵，且对角线元素全不为零;系数矩阵A的一个分解:
(7）
这里D为A的对角素构成的对角矩阵,为严格下三角阵，为严格上三角阵。
SOR迭代法的矩阵形式为 X（k+1）=HSOR X(k）+dSOR
式中 :
HSOR=（D—L)-1 [ (1-） D + U]
d SOR= ( D-L ）-1b。
计算公式为：

i=1,2，3，…n k=0,1,2，… （8)
然后按相反的次序(i= n，n-1，…1）用向后的SOR方法逐点计算

i= n，n-1，…1；k=0,1，… (11）
若AX= b中A是正定的，且0﹤﹤2，则SSOR法收敛。
2。4 消元法
消元法是初等数学中求解低阶多元线性方程组的方法．此时线性方组必须是适定方程组.一般是用于二元一次 或三元一次方程组．当未知元增多时.计算效率低甚至无法求解.
2.5 克拉默法则
当系数行列式不为零时．适定方程组有惟一 解 ．其解如（12)式所示 ：xi= Di /D i=1，2,…，n (12)
其中D是系数行列式，Di是在系数行列式基础之上结合方程组右边常数形成的新行列式。在此法则中。行列式的计算显得非常重要。用行列式的性质计算行列式最为有效。 对于二 、三阶行列式可以利用对角线法则计算。
克拉默法则克服了消元法计算效率低甚至无法计算多元一次方程组的缺点．但是对于系数行列式等于零以及欠定或者超定方程组的情况，它是无能为力的.事实上，当未知元数过多时 (如未知元数≥5）。克拉默法则的计算效率就很低。
2.6 逆阵乘积法
对于适定方程组．可以把它表达成矩阵方程的形式：AX=b
解矩阵可以利用逆阵乘积法求出：
A-1 b = X
矩阵运算的实质是把矩阵当作一个“量 ”来运算。使 普通数 的运算有很大的简化。但是该方法的前提是A可逆。本质上仍然是系数行列式｜A|≠0．对于阶数比较高的系数矩阵．直接求解逆阵也是比较 困难的（利用初等变换可以降低求解难度)。当|A｜=0时,或者对于欠定或者超定方程组，逆阵乘积法仍然是无能为力的或不适合的。
2.6 初等变换法
对于欠定或超定方程组的求解。初等变换法是最直接、 最简单的方法。同时该方法也能用于适定方程组的求解。因此，初等变换法是一种求解线性方程组的普适方法和最 基本方法 。秩是矩 阵的本征参 数．利用系数矩阵的秩和增广矩阵秩的关系,可以判断线性方程组解的情况：无解，惟一解和无穷多解。对增广矩阵进行初等变换．可以得到增广矩阵的等价矩阵．从而得到了原来方程组的等价方程组。由于等价矩阵已变换成阶梯形矩阵或最简形矩阵,所以等价方程组变 得非常简单。可以方便地求出解。初等变换也是求逆阵的一种直观方法．所以可以和阵乘积法配合运用。
2。7 利用向量空间概念求解线性方程组
如果把齐次线性方程组的解集看作是一个向量空间 。由于一个向量组 (向量空间)与它的一个最大无关组等价。则只需求出解集的一基础解系即可确定齐次线性方程组的解。求解的方法仍然采用初等变换法．解的形式可以用通解来表达．更能说明解的本质。尤其是无穷多解。基础解系中线性无关的向量形成解向量空间的一个最大无关组 。在非齐次方程组求解中．只要求出对应齐次方程组的通解。然后加上非齐次方程组的一个特解即可．而且这种通解和特解可以在一次初等变换中求出。
2。7 Krylov子空间方法
Krylov子空间方法是解决大型线性方程组的一类十分有效的方法 ，其主要代表是对称正定线性方程组的共轭梯度法和非对称线性方程组的GMRES方法。20世 纪 50年代 初 期 由 Hestenes和 Stiefel首先提出的共轭梯度法 （或简称 CG法 )。近20年来有关的研究取得了前所未有的发展，目前有关的方法和理论已经相当成熟，并且成为求解大型稀疏线性方程组十分有效的一类方法。
通过对经典迭代法的总结 ，发现 SOR迭代法在松弛因子取最优松弛因子opt的情况下,迭代收敛很快。但是，计算opt还需要求得对应的 Jacobi迭代矩阵的谱半径,这常常是比较困难的。
考虑线性方程组： Ax= b （13）
的求解问题，其中A是给定的n阶对称正定阵,b是给定的n维向量， x是待求的n维向量。为此定义二次泛函数
Ψ( x） = xTA x -2bTx （14）
则可以证明求方程组 （13)的解 等价于求二次泛函(14）的极小点。
由此,给定了初始向量x（0），按某一方向去求式（14)的极小值点x(1），就得到下一个迭代值x(1),再由x（1）出发 ，求x（2)等等。这样来逼近精确解x*.若取求最小值的方向为Ψ在x（k-1)（k =1， 2，⋯）处的负梯度方向，就是所谓的最速降法。 然而理论和实际计算表明这个方法的收敛速度较慢，共轭梯度法则是在x（k-1)处的梯度方向r(k-1)和
这一步的修正方向p(k—1)所构成的二维平面内，寻找使Ψ减少最快的方向作为下一步的修正方向，即求极小值的方向p（k) （其第一步仍取负梯度方向）。计算公式为p（1）=r（0）=b— x（0）
再逐次计算
qk=（ r(k-1）,r（k-1））/（A p(k)， p（k)）,(（x，y）表示向量xy的内积
x（k）= x(k-1）+ qk p(k）
r（k）= r(k—1）- qk p（k）
ek=(r(k)，r(k）)/（r（k—1)，r（k-1)）
p(k+1）= r（k)+ ek p（k）（k=1,2，…）
可以证明当i≠j时，（p(i)，A p（j)）=0,(ri,rj)=0
从而p（1）,p(2)…．形成一共轭向量组；r（0），r（1）…形成一正交向量组。后者说明若没有舍入误差的 ，至多n次迭代就可以得到线性方 程组(14）的精确解。然而，在实际计算中一般都有舍入误差，所以r（0）,r（1）…并不是真正相交，而x(0），x（1)…等也只是逐步逼近式（14）的真解，故我们将共轭梯度法作为迭代法来使用。
3 数值实验
下面通过一个数值实例来比较Jacobi方法、Gauss-Seidel方法、SOR方法的收敛速度。
数值实验：我们取一个系数矩阵为64阶的线性方程组：A X= b.其中：
精确解:X=（1 1 1 1… 1）T1×64; X（0) （0 0 0 0… 1）T1×64； 精度为:0.0001.
执行C语言程序interrator-method. C得到result.txt文件.
得到三种迭代法的对比表1。
表1 三种方法的迭代性能

Method
IT CPU Time／ms
‖x(k）—x∞‖∞
Info
Jacobi
18 3.00000
8．19628×10—5
几乎比 Jacobi收敛速度快一倍
Gauss-Seidel
10 2.00000
1．68O24×10—5
SOR
with=1。0
10 2.00000
1．68024×10-5
退化成 Gauss
_ Seidel迭代
使 SOR迭 代
达到最快
的取值于SOR迭代的收敛速度
影响很大
SOR with=opt
6 1.00000
5．47552×10-6
SOR
with=1。5
17 3。00000
2．09836×10-5

4 优点
克拉默法则、逆阵乘积法只能求解系数行列式不为零的适定方程组 ；初等变换法可以直观地解决所有类型的超定、欠定、适定方程组,是一种普适的方法；利用向量空间概念求解线性方程组，更能从本质上把握线性方程组的解的性质。应用MATL B语言编程可以轻松实现这些求解方法。
实用共轭梯度法主要有以下优点 :
1）算法中，系数矩阵A的作用仅仅是用来由已知向量P产生向量 =A,这就充分利用了A的稀疏性。而且，对某些提供矩阵A较为困难而由已知向量产生向量=A又十分方便的应用问题是很有益的 。
2）不需要预先估计任何参数就可以计算,这与SOR方法等不同。
3)每次迭代所需的计算,主要是向量之间的计算，便于并行运算
5 附录
MATLAB语言是一种以矩阵运算为基础的计算语言。对于实现线方程组的求解非常方便。对一个四元一次方程组的求解．可以用克拉默法则和逆阵乘积法来实现,程序如下 ：
tic；
D=[1 1 1 1；1 2 -1 4；2 -3 -1 —5;3 1 2 11］
det(D）
b=［5—2—2 0］；
Dl=[5 1 l l；—2 2 -1 4；-2 -3 —1 —5；0 1 2 11］；
D2=[1 5 1 l; 1 —2 —1 4；2 -2 -1 —5；3 0 2 11］；
D3=[1 1 5 1;1 2—2 4;2 -3 —2 -5;3 1 0 11］；
D4=[1 1 1 5;1 2 -1 —2；2 -3 —1 -2；3 l 2 0];
X1= det (D1)／det(D）;
X2= det (D2)／det（D）;
X3= det (D3）／det（D）;
X4= det (D4)／det（D)；
X=inv（D)＊b；
Toc
其中克拉默法则用行列式除法( Xi=det（Di)／det（D)）来实现;逆阵乘积法用(X=inv(D）＊b）来实现;det（D)是系数矩阵D的行列式运算；inv(D）是D的逆阵运算.
上例中．系数矩阵D不为零，可以用克拉默法则和逆阵乘积法来求解。当系数行列式为零时，只能用初等变换来求解。对于初等变换．利用阶梯生成函数命rreff也可以轻松地实现，举例如下 ：
tic；
A=［3—4 3 2 -1;0 —6 0 —3 -3；4 -3 4 2 —2；1 1 1 0 -1；—2 6 -2 1 3］；
b=［2 -3 2 0 1］;
B=［A,b］,[UB，ip］=rref（B)
U0=UB([1：5］,［1：5]），d=UB（:,6)
toc
UB为经过初等变换以后的行阶梯矩阵.可以轻松地求出方程组的解。
tic、toc为计时函数,CPU-pentium1.7GHz、512MB内存、MATLAB6.5的运算环境中，两个程序的运行时间分别为0。031s和0。015s。可见,MATLAB语言实现线性方程组的求解具有程序简单、直观的特 点。同时还具有计算效率高的优点。在实际计算中摆脱了系数矩阵阶数、未知元数等的限制 。
参考文献
[1］同济大学应用数学系．数学-—线 性代数［M］．4版.北京 :高等教育出版社．2003．
[2]陈怀琛，龚杰民．线性代数实践及MATLABA人民[M］．北京：电子工业出版社．2005．