线性代数方程组数值解法及MATLAB实现综述
线性代数方程组数值解法及MATLAB实现综述廖淑芳20122090 数计学院12计算机科学与技术1班（职教本科）一、分析课题
随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。
直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。
迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。
二、研究课题-线性代数方程组数值解法
一、直接法
1、Gauss消元法
通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x向量。
1.1消元过程
1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。
11121121222212n n n n nn
n a a a b a a a b a a a b ??
? ? ?
???(1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322
(3)(3)(3)3333()()000
00
n n n
n n nn
n a a a a b a a a b a a b a b ?? ?
? ? ? ?
???
步骤如下：
第一步：1
11
1,2,,i a i i n a -?
+=第行第行
11121121222212
n n n n nn
n a a a b a a a b a a a b ??
? ? ?
???1112
11(2)(2)(2)22
22
(2)(2)(2)2
00n
n
n nn
n a a a b a a b a a b ?? ? ? ? ???
第二步：(2)2
(2)222,3,
,i a i i n a -?+=第行第行 111211(2)(2)(2)2222(2)(2)(2)2
00n n
n nn
n a a a b a a b a a b ?? ? ? ? ???11
12
1311(2)(2)(2)(2)222322
(3)(3)(3)33
33(3)(3)(3)3
0000
0n
n n
n nn
n a a a a b a a a b a a b a a b ?? ? ? ? ? ?
???
类似的做下去，我们有：
第k 步：()
()k ,1,
,k ik
k kk
a i i k n a -?+=+第行第行。
n －1步以后，我们可以得到变换后的矩阵为：
11121311(2)(2)(2)(2)222322
(3)(3)(3)3333()()00000
n n n
n n nn
n a a a a b a a a b a a b a b ?? ? ? ? ? ? ??
?
注意到，计算过程中()
k kk a 处在被除的位置，因此整个计算过程要保证它不为0。
所以，Gauss 消元法的可行条件为：()
0k kk
a ≠。 就是要求A 的所有顺序主子式均不为0，即
11
11
det 0,1,,i i ii a a i n a a ?? ?
≠= ? ???
因此，有些有解的问题，不能用Gauss 消元求解。
另外，如果某个()
k kk a 很小的话，会引入大的误差。
例 用Gauss 消去法解方程组：
（1）1234123341518311151111631112x x x x -??????
? ? ?---- ? ? ?= ? ? ? ? ? ?
-?????? （1）对增广矩阵进行初等变换
122133144
3
()21
()121
()4
123
34153715123341505
22218311155321911116044343111277
70
0444E E E E E E E E E +→-+→-+→-??
??-? ?- ? ?---- ? ???????→ ? ? ? ?
?-?? ?-- ?
?
?
233244344
5
()6
77()()6
11
123341233
4151537370505151522222211
29
11
29
00
0111136367357910000
003633E E E E E E E E E +→+→-+→--????
?
? ? ?-- ? ? ? ??????→??????→ ? ?
?
? ?
?
? ??
??
?
得等价方程组
123423434412334153715
522211291136
910
33x x x x x x x x x x -++=??
?-++=??
?
+=??
?=??
回代得40x =，33x =，22x =，11x =。
第一步：将(1)/3使x 1的系数化为1,再将(2)、(3)式中x 1的系数都化为零，即由(2)-2×(1)(1)
得
由(3)-4×(1)(1)得
第二步：将(2)(1)
除以2/3，使x 2系数化为1，得
再将(3)(1)式中x 2系数化为零，由(3)(1)-(-14/3)*(2)(2)，得
第三步：将(3)(2)
除以18/3，使x 3系数化为1，得
经消元后，得到如下三角代数方程组：
)1(321)1( (23)
1
32=++
x x x )1(32)2( (034)
32=+x x )
1(32)3( (63)
10
314-=--x x )
2(32)2(......02=+x x )2(3)3( (63)
18
-=x )
3(3)3(......1-=x
1.2回代过程
由(3)(3)得 x 3=1,将x 3代入(2)(2)得x 2=-2,将x 2 、x 3代入(1)(1)得x 2=1,
所以，本题解为[x]=[1,2,-1]T
1.3高斯消元的公式
综合以上讨论，不难看出，高斯消元法解方程组的公式为 第一步，消元
（1） 令
a ij (1) = a ij , （i,j=1,2,3,…,n ）
b i (1) =b i , （i=1,2,3,…,n ） （2） 对k=1到n-1，若a kk (k)≠0，进行
l ik = a ik (k)
/ a kk (k)
, （i=k+1,k+2,…,n ） a ij (k+1) = a ij (k) - l ik * a kj (k), （i,j= k+1,k+2,…,n ） b i (k+1) = b i (k) - l ik * b k (k), （i= k+1,k+2,…,n ）
第二步，回代
若a nn (n) ≠ 0
x n = b n (n) / a nn (n)
x i = (b i (i) – sgm(a ij (i) * x j )/- a ii (i) ,（i = n-1,n-2,…,1）,( j = i+1,i+2,…,n )
2 、LU 分解法
求解线性代数方程组除了高斯消元法外，还常用LU 分解法（三角形分解法）。LU 分解法的优点是当方程组左端系数矩阵不变，仅仅是方程组右端列向量改变，
即外加激励信号变化时，能够方便地求解方程组。矩阵的三角分解法可分为直接三角分解法,列主元三角分解法,平方根法,三对角方程组的追赶法。下面讨论直接三角分解法。
设n阶线性方程组Ax=b 。假设能将方程组左端系数矩阵A，分解成两个三角阵的乘积，即A=LU ，式中，L为主对角线以上的元素均为零的下三角矩阵, 且主对角线元素均为1的上三角矩阵；U为主对角线以下的元素均为零
所以有，LUx=b令Ux=y，则Ly=b
由A=LU,由矩阵的乘法公式：a
1j = u
1j
, j=1,2,…,n
a i1 = l
i1
u
11
, i=1,2,…,n
推出u
1j = a
1j
, j=1,2,…,n
l i1 = a
i1
/u
11
, i=1,2,…,n
这样就定出了U的第一行元素和L的第一列元素。
设已定出了U的前k-1行和L的前k-1列，现在确定U的第k行和L的第k列。由矩阵乘法：
. -
当r>k 时，l kr =0, 且l kk =1，因为
所以，
同理可推出计算L 的第k 列的公式：
因此得到如下算法——杜利特（Doolittle ）算法:
(1) 将矩阵分解为A=LU,对k=1,2,…,n ；j=k,k+1,…n ； i=k,k+1,…n ；
公式1
(2) 解L y =b
(3) 解U x =y
对大规模稀疏问题，如果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU 分解，可以保障单位下三角矩阵L 及上三角矩阵U 仍为带状结构.
∑==n
r rj
kr kj u l a 1
∑=+=-=n r rj
kr kj kj n
k k j u l a u 1
,...,1,∑=+=n
r rj
kr kj kj u l u a 1
∑=+=-=n
r kk
rj kr ik ik n
k k i u u l a l 1
,...,1,/)(∑-==-=1
1,...,2,1:2k r r
kr k k n
k y l b y 公式∑+=-=-
=n
k r kk
r
kr k k n n k u x u
y x 1
1
,...,1,/)(3：公式???
?
?
?
???
=-=-=∑∑==1/)(u 11
kk kk
n
r rj ky ik jk n
r yj
ki kj ki l u u l a l u l a
3、直接三角分解法
Gauss 消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把Gauss 消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些修正和简化，这些方法可统称为直接三角分解法。
矩阵的三角分解
设A 的顺序主子式0(1,2,
,)k k n ?≠=，则可建立线性方程组Ax b =的Gauss
消去法与矩阵分解的关系，即矩阵A 的LU 分解。这个问题前面已经讲的比较详细了，此处不再赘述。
Doolittle 分解法
首先假设A 的顺序主子式k ?都不为零，则A 可作Doolittle 分解，即A LU =，其中L 是单位下三角阵，有1ii l =，i j <时0ij l =；U 是上三角阵，i j >时0ij u =。仔细写出为
1112
111121212222122
2121211
1n n n n n n nn n n nn a a a u u u a a a l u u A LU a a a l l u ??????
?
??
? ? ???
=== ? ???
?
?????????
（2.11） 在前面逐步推导L 和U 的元素公式都要借助于有关的()
k ij a 来表示。现在强调
指出，只要从给定的A 通过比较（2.11）式的两边就可能逐步地把L 和U 构造出来，而不必利用Gauss 消去法的中间结果，这种方法称为Gauss 消去法的紧凑格式。
根据矩阵的乘法规则，比较（2.11）式的两边可得
min(,)
1
i j ij ik kj k a l u ==
∑
，,1,2,
,i j n = （2.12）
先算出U 的第1行，在（2.12）令1i =，由111l =，10(1)j l j n =<≤可得
11j j u a =，1,2,
,j n =
接着在（2.12）中令1j =，由1111i i a l u =?，从而算出L 的第1列
1111111//i i i l a u a a ==，2,
,i n =
若U 的第1至1k -行和L 的第1至1k -列已经算出，则由（2.12）可计算出
U 的第k 行元素
1
1k kj kj kr rj r u a l u -==-∑，,1,
,j k k n =+ （2.13）
接着再计算出L 的第k 列元素
1
1
k ik ir rk
r ik kk
a l u l u -=-=
∑，1,,j k n =+ （2.14）
解方程组Ax b =就化为求解LUx b =，分两步求解。第一步解Ly b =，其中L 为下三角阵，只要用逐次向前代入的方法；第二步解Ux y =，其中U 为上三角阵，用逐次向后代入方法即可解除x 。
例 用Doolittle 方法求解：
12346
2
1
16241011141510135x x x x -?????? ? ? ?
- ? ? ?
=
? ? ?- ? ? ?
---?????? 解：第1步算出L 和U ：
111311
16
5119
161037L ?? ? ? ? ?=
? ? ?-
- ???，6211102
133337910
1019174U -??
?
? ? ?=-
? ? ?
??
?
第2步解Ly b =得：
231916,3,,574T
y ?
?=- ??
?
第3步解Ux y =得：
(1,1,1,1)T x =--
二、 迭代法
迭代法具有的特点是速度快。与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根。 对方程组：Ax b =做等价变换：x Gx g =+ 如：令A M N =-，则
11()Ax b M N x b Mx b Nx x M Nx M b --=?-=?=+?=+
则，我们可以构造序列(1)() k k x G x g +=+ 若()*k x x →* **x G x g Ax b ?=+?= 则，我们可以构造序列(1)() k k x G x g +=+ 若()*k x x →* **x G x g Ax b ?=+?= 同时：(1)()()**(*)k k k x x Gx Gx G x x +-=-=-
1(0)(*)k G x x +=
=-
所以，序列收敛0k G ?→ 与初值的选取无关 1. Jacobi 迭代
11111
11
n n n nn n n
a x a x
b a x a x b +
+=???
?++=?
11221111
2211233122211 111()
1() 1()
n n n n n n n n n n nn x a x a x b a x a x a x a x b a x a x a x b a ---?=++-??
-?
=+++-?????
-?=++-??
(1)()()11221111(1)()()()2
2112331222(1)()()11 111()1
() 1()k k k n n k k k k n n k k k n n n n n n nn x a x a x b a x a x a x a x b a x a x a x b a +++---?=++-??
-?=+++-?∴?
??
-?=++-??
格式很简单：
1
(1)
()()1
1
1()i n
k k k i
ij j ij
j i j j i ii x a x a x
b a -+==+-=+-∑∑
迭代矩阵
记A D L U =--
1100nn a D a ?? ?=
? ??
?
211
1
00
00n nn a L a a -?? ?- ?
?= ? ? ?--??121100000n n n a a U a ---??
?
? ?= ?
-
?
??
?
易知，Jacobi 迭代有
()D L U x b --= ()Dx L U x b =++ 11()x D L U x D b --=++
11100 B () , D L U I D A f D b ---∴=+=-=，0B 是简单迭代的迭代矩阵。
看上述公式和过程很抽象，我们来举个简单例子：
111122133121122223323113223333
a x a x a x
b a x a x a x b a x a x a x b ++=??
++=??++=? （0ii a ≠） 1
112213311
2221123322
33311322331()1()1
()x b a x a x a x b a x a x a x b a x a x a ?=--???=--???=--??
得
得上述变换的方程后，我们任取一向量(0)(0)(0)(0)12
3(,,)T
x x x x =作初始近似，代入， (1)(1)(1)(1)()()()()123123(,,),,(,,)T
m m m m T
x x x x x x x x ==
假设上述方程的准确解是123(,,)T a a a a =
那么如果lim ,m m x a →∞
=?我们就说上述迭代是收敛的。
2. Gauss －Seidel 迭代
在Gauss －Seidel 迭代中，使用最新计算出的分量值
(1)()()
1
1221111
(1)(1)()()2
2112331222(1)(1)(1)11 111()1
() 1
()k k k n n k k k k n n k k k n
n n n n n nn x a x a x b a x a x a x a x b a x a x a x b a ++++++---?=++-??
-?=+++-?∴?
??
-?=++-?? 1
(1)
(1)()1
1
1()i n
k k k i
ij j ij
j i j j i ii x a x a x
b a -++==+-=+-∑∑
易知，Gauss －Seidel 迭代有
()D L U x b --= ()Dx L U x b =++ 11()x D L U x D b --=++
1110 () , G D L U I D A g D b ---∴=+=-=，0G 是简单迭代的迭代矩阵。
二种方法都存在收敛性问题。 收敛条件
迭代格式收敛的充要条件是B
G
的谱半径<1。
有例子表明：Gauss-Seidel 法收敛时，Jacobi 法可能不收敛；而Jacobi 法收敛时， Gauss-Seidel 法也可能不收敛。
若为严格对角占优的话则是收敛的，如130.10.51
20.0250.3110.030.142
???? ? ?
? ? ? ? ?????
，
假如方程组不满足收敛，有时候对其进行变换，可以改变收敛性。 如求下述方程组的解：
(0)1807109,8,(0,0,0)9117T A b x -????
? ?
=-== ? ? ? ?--????
9117180,71098A b --????
? ?=-= ? ? ? ?-????
格式
(1)()()123(1)
(1)1
1(1)(1)1
11/9(7)1/8(7)1/9(8)
k k k k k k k x x x x x x x +++++?=++?=+??=+? 结果
(1)(2)(3)
(4)(0.7778,0.9722,0.9753)(0.9942,0.9993,0.9994)(0.9999,0.9999,0.9999)
(1.0000,1.0000,1.0000)T T T
T
x x x
x ====
0()1B ρ<，G ρ（）<1也满足收敛。
如122111221A --?? ?=-- ? ?-??
Jacobi ，0B 的特征方程为221
102
2
λ
λ
λ
---=-
解得1230λλλ===，所以用Jacobi 迭代法收敛。
Gauss-Seidel ，0G 的特征方程为22
1022λλ
λλλλ
---=--
解得1230,22λλλ==-=-
0()21B ρ=->, 所以Gauss-Seidel 迭代法是不收敛的。
3. 松弛迭代 记()(1)()k k k x x x +?=- 则(1)()()k k k x x x +=+?
可以看作在前一步上加一个修正量。若在修正量前乘以一个因子ω，有
(1)()()k k k x x x ω+=+? (1)()(1)()()k k k k x x x x ω++=+-
对Gauss －Seidel 迭代格式(1)1(1)()()k k k x D Lx Ux b +-+=++
(1)()1(1)1()1()()k k k k k x x D Lx D Ux D b x ω+-+--=+++- (1)()1(1)1()1(1)()k k k k x x D Lx D Ux D b ωω+-+--=-+++
写成分量形式，有
(1)()()()
111221111(1)
()()()()2
22112331222(1)()()()
11 111(1)()(1)()1(1)()k k k k n n k k k k k n n k k k k n n n n n n n nn x x a x a x b a x x a x a x a x b a x x a x a x b a ωωωωωω+++---?=-+++-??
?=--+++-??
??
-?=-+++-?
?
关于直接法和迭代法的例题： 1用选主元三角分解法求解下列方程组
1234123
1234123433246
2421022431742268x x x x x x x x x x x x x x x +++=??++=??
+++=??+++=? 12343324624201022431742268x x x x ??????
? ? ?
? ? ?
=
? ? ?
? ? ?
??????3324
62420102243174
2268?? ?
?
?
?
??
列主
4226
84
2268422681131363
136222420108111130131112243172233
3
324
633113
100
13422242??
????
?
?-- ? ?
? ?
? ?→→ ?- ?
? ? ? ??? ? ?
--????
这里100042268110003132
6,,1
1810001112333310001014
2L U y ??
???? ?
? ? ?- ? ? ?=== ?
?- ? ? ? ? ?-?? ?????
再通过Ux y =，得x 的解。
123424233x x x x ???? ? ?- ? ?= ? ? ? ?-??
?? 2用Jacobi 迭代法求解下列方程组，并且使精度为310-。
1234
0.240.0880.0930.1590.040.08420x x x -?????? ??? ?-= ???
? ??? ?-??????以4，3，4分别除3方程两边 12310.060.0220.0310.0530.010.0215x x x -?????? ??? ?-= ??? ? ??? ?-?????? 即11223300.060.0220.03
00.0530.010.0205x x x x x x -???????? ? ??? ?
=-+ ? ??? ? ? ??? ?-????????
有Jacobi 迭代式123()
(1)1(1)()2(1)()300.060.0220.03
00.0530.010.0205m m m m m m x x x x x x +++????-???? ? ? ? ?
=-+ ? ? ? ? ? ? ? ? ?-????????
可证明此迭代格式是收敛的，极限值即为解。
40.240.080.060.020.0930.150.030.050.040.0840.010.02λ
λλλλλ--???? ? ?-→- ? ? ? ?--????22
20.06
0.02
0.0018
0.050.0006
0(0.06)(30.010.0009)00
3(0.0018)λλλλ
λλλλλ?? ?- ? ?-+→- ? ?++- ?
?-?
? 解得1230.06,λλλ=-=
=
由此可知，用Jacobi 迭代法是收敛的。 取(0)(2,3,5)T x =为初始值，迭代得，
23(1)1(1)(1)00.060.0222 1.920.03
00.0533 3.190.010.02055 5.04x x x ??-????????
?
??? ? ?=-+= ? ??? ? ? ? ??? ? ? ?-??????????
23(2)1(2)(2)00.060.02 1.922 1.9090.03
00.05 3.193 3.1940.010.020 5.045 5.045x x x ??-???????? ?
??? ? ?=-+= ? ??? ? ? ? ??? ? ? ?-?
???????
??…
其实23(3)1(3)(3)00.060.02 1.9092 1.9090.03
00.05 3.1943 3.1940.010.020 5.0455 5.045x x x ??-???????? ? ??? ? ?
=-+= ? ??? ? ? ? ??? ? ? ?-?????????? 即在310-，(3)(2)x x =，可以停止计算了。所以(3)x 即为近似解。
分析：此题题目可以直接说是解方程组，然后求解过程我们可以先验算Jacobi 和Gauss －Seidel 迭代的收敛情况，再选择算法，一般来说Gauss －Seidel 迭代要比Jacobi 迭代快。
三、用Matlab 软件求解线性方程组 一．高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[]' [n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); 消去 end
x=[0 0 0 0]' 回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
2.列主高斯消去法
*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件
a=[]
d=[] '
[n,n]=size(a);
c=n+1
a(:,c)=d; (增广)
for k=1:n-1
[r,m]=max(abs(a(k:n,k))); 选主
m=m+k-1; (修正操作行的值)
if(a(m,k)~=0)
if(m~=k)
a([k m],:)=a([m k],:); 换行
end
a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去end
end
x=[0 0 0 0]' 回带
x(n)=a(n,c)/a(n,n);
for g=n-1:-1:1
x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)
end
二迭代法
J迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
d=[-12;20;3]
x=[0;0;0]; 初始向量stop=1.0e-4 迭代的精度L=-tril(a,-1)
U=-triu(a,1)
D=inv(diag(diag(a)))
X=D*(L+U)*x+D*d; J迭代公式n=1;
while norm(X-x,inf)>=stop
x=X;
X=D*(L+U)*x+D*d;
n=n+1;
end
X
n
G-S迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; 初始向量
d=[-12;20;3]
stop=1.0e-4
L=-tril(a,-1)
U=-triu(a,1)
D=(diag(diag(a)))
X=inv(D-L)*U*x+inv(D-L)*d; G-S迭代公式
n=1;
while norm(X-x,inf)>= stop 时迭代中止否则继续x=X;
X=inv(D-L)*U*x+inv(D-L)*d;
n=n+1;
end
X
n
SOR迭代公式
a=[5 2 1;-1 4 2;2 -3 10]
x=[0;0;0]; 初始向量
d=[-12;20;3]
stop=1.0e-4