计算方法 第3章 线性方程组数值解法

第三章
线性方程组的数值解法
线性方程组是应用最为广泛的数学模型，很多复杂问题中都含有线性方程组子问题，因此讨论线性方程组问题的求解很有必要，本章将讨论线性方程组的数值解法。

线性方程组的一般形式：
⎪⎪⎩⎪⎪⎨⎧=+++=+++=+
++
n n
nn n n n n n n b x a x a x a b x a x a x a b x a x a x a 2
2112
2222121112
121
11 (3.1)
如果记：
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭
⎫
⎝⎛=⎪⎪
⎪⎪
⎪
⎭
⎫
⎝⎛=n n nn n n n n b b b b x x x x a a a a a a a a a A 21212
1
222
21
11211
,
,
这里A 称为系数矩阵，x 称为解向量，b 称为右端项，则得线性方程组的矩阵形式： b Ax = (3.2)
求解线性方程组问题(3.1)或(3.2)的数值方法可分为两类：直接解法和迭代解法，其直接
解法是通过有限次初等运算，求得其解，虽然直接解法的推导过程都是无误差的，但是由于计算机的运算都是有舍入误差的，所求解其实是一个有误差的近似解；迭代解法则是从某个初始近似解出发，按照一个确定的迭代公式得到一个更好的近似解，反复迭代，直到求得一个满足精度要求的近似解。

本章首先讨论线性方程组的直接解法然后再介绍迭代解法。

3.1 消去法
1. 顺序Gauss 消去法
首先回顾一下线性代数中所讲的线性方程组消去法过程，然后归纳出消去法的数值算法，请看如下的例子：
例3.1 求解线性方程组
⎪⎩⎪
⎨⎧=-+=+-=-+2
24056242321
321321x x x x x x x x x
解：求解线性方程组的第一阶段称为消元过程，其方法是：第2个方程减去第1个方程的
2
1
倍，第3个方程减去第1个方程的2倍，得 ⎪⎩
⎪
⎨⎧-=+--=+-=-+102736362423232321x x x x x x x 第3个方程减去第2个方程的
3
7
倍，得 ⎪⎩
⎪
⎨⎧-=--=+-=-+3123636242332321x x x x x x 这一过程就是消元过程，即把方程化为等价的上三角方程(对角线下变为0)。

第一个阶段完成后，进入第二个阶段，称为回代过程，其方法是：先由第3个方程解出
3x ，将3x 代入第2个方程解出2x ，再将3x 和2x 代入第1个方程解出1x 也就解出所有的未
知量。

如下就是所求的解：
4
1
,23,41321===
x x x
归纳以上求解方法，求解线性方程组包括两个过程，消元过程和回代过程。

首先给出回代过程的算法，回代过程其实是一个特殊形式的方程组的求解方法，就是一
个上三角线性方程组的求解方法，如：
⎪⎪⎪⎪⎩
⎪
⎪⎪⎪⎨
⎧==+++=+++++=++++++++++++n n nn k
n kn k k k k kk n n k k k k n n k k k k b x a b x a x a x a b x a x a x a x a b x a x a x a x a x a .............................................................................11,2211,222221
111,11212111 (3.3)
第1步：根据第n 个方程解出n x ，得
nn
n
n a b x =
第2步：根据已求出的n x 和第n -1个方程求1-n x ，得
1
,1,11-----=
n n n
n n n n a x a b x
假如已经求出11,,,+-k n n x x x ，代入第k 个方程可以求出k x ，得：
kk n
k j j kj
k k a x a
b x /)(1
∑+=-
=
(3.4)
回代过程就是对1,2,,1, -=n n k 实施这一公式，注意必须从后向前计算方可，所以此过程叫做回代过程，(3.4)式也是求解下三角方程的一般公式，可看作一个独立的算法。

算法3.1(上三角线性方程组的回代算法)
0) [初始化] 设置上三角方程系数矩阵A ，右端项向量b 1) [回代过程] 对于1,,1, -=n n k 循环
1.1) k k b x =
1.2) 对于n k j ,,1 +=循环
j ij k k x a x x -=
1.3) kk k k a x x /= 2) [算法结束]
再看Gauss 消去法的消元过程，对于线性方程组
⎪⎪
⎩⎪⎪⎨
⎧=+++=+++=+++n
n nn n n n n n n b x a x a x a b x a x a x a b x a x a x a (22112)
222212********* 的消去法本质上可看作将其增广矩阵用初等行变换化为梯形矩阵的过程。

为清楚的表示每次消元前后系数矩阵和右端项的状态，通常以)
(k A
和)
(k b
表示系数矩阵
和右端项，其元素分别记作n j i b a k i k ij ,,2,1,,,)
()
( =，其上标)(k 表示在第k 次消元前的状态，其初始增广矩阵为：
⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛)1()
1()1(2)1(1
)1(2)1(2)1(22)1(21)1(1)1(1)1(12)1(11n nn
n n n
n b a a a b a a a b a a a
(3.5)
首先对第1列消元，也就是对增广矩阵做初等行变换将元素)
1(1)1(21,,n a a 约化为0，
希望消元之后形如：
⎪⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛)2()2(2
)2(2)2(3)2(3)2(32)2(2)
2(2)2(22)1(1)
1(1)1(12
)1(110
00n n n n
n n
b a a b a a b a a b a a a
(3.6*)
为此，对n i ,,3,2 =，令消元因子)
1(11
)1(11/a a m i i =，将第i 行减去第1行的1i m 倍，则可将(3.5)式变为(3.6*)式，之后的论述中，在不引起混淆的情况下，将(3.6*)式记为：
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝
⎛n nn
n n n n b a a b a a b a a b a a a
233322222111211
000 (3.6)
第2列消元，对n i ,,4,3 =，令2222/a a m i i =，将第i 行减去第2行的2i m 倍，则有
⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛n nn n n n n
b a a b a a b a a a b a a a a
3333322232211131211
00000 对第3、4、…、n -1列消元，最后得到上三角方程组的增广矩阵为
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝
⎛n nn
n n
n
b a b a a b a a a b a a a a
000033332223221113
1211 以上的消元过程可用一个统一公式表示为：
.1,...,2,1,,...,1,...,1,/-=+=⎪
⎭
⎪⎬⎫⎪⎩⎪
⎨⎧+=-⇐-⇐=n k n k i n k j a m a a b m b b a a m kj ik ij ij k ik i i kk
ik ik (3.7)
通常在编程求解线性方程组时，消元因子),,1(n k i m ik +=存储于矩阵A 对角线以下相应的位置上，因为此时这些位置的元素将被消元为0。

算法3.2 (顺序Gauss 消去法)
0) [初始化] 置系数矩阵A ，右端项向量b ； 1) [消元过程] 对于1,,2,1-=n k 循环
对于n k i ,,1 +=循环 1.1) kk ik ik a a a /= 1.2) k ik i i b a b b -=
1.3) .,,1,n k j a a a a kj ik ij ij +=-= 2) [回代过程] 执行算法3.1，即
1,2,1,/)(1
-=-
=∑+=n n k a x a
b x kk n
k j j kj
k k ，
3) [算法结束]
定理 3.1 顺序Gauss 消去法能够顺利进行的充要条件是系数矩阵A 的顺序主子式
n k D k ,,2,1,0 =≠，且)
()2(22)1(11k kk
k a a a D =，特别有)()2(22)1(11n nn n a a a D A ==。

证：从前述的算法可知，消元过程能够顺利进行的充要条件是步骤1.1)中对角线元素
0)
(≠k kk a ，以下只要证明顺序主子式为)()2(22)1(11
k kk k a a a D =即可。

事实上，顺序消元过程所做行变换并不改变顺序主子式的值，而在完成第k 次消元过程后，增广矩阵被约化成：
⎪⎪⎪⎪⎪
⎪⎪⎪⎭
⎫
⎝
⎛+++++++)()
()(1,)(1)
(,1)
(1
,1)()
()(1,)
()1(1)
1(1)1(1,1)
1(1)1(11k n
k nn
k k n k k k n k k k k k k k kn
k k k k kk n k k b a a b a a b a a a b a a a a
显然第k 个顺序主子式为：
)
()2(22)1(11)
()
2(2)2(22)
1(1)1(12
)1(11
k kk k kk
k
k
k a a a a a a a a a D
==
2．顺序Gauss 算法的计算量
通常我们在度量一个有限步终止的算法的计算量时，主要考虑整个计算过程中使用乘除法运算的总次数，在以后统计算法的计算量时，我们仅统计算法中乘除运算的次数。

1) 消元过程的计算量
6
5236)12)(1()1()
2())(2()2(2
31
1
11
111
n
n n n n n n n k k k n k n k n n k n k n k n k i -
+=--+
-=+=-+-=+-∑∑∑∑-=-=-=+=
2) 回代过程的计算量
222)1()1(21
1
n
n n n k k n n
k n k +=+==+-∑∑== 3) 总的乘除法运算的计算量为
总计算量=)3
1(333332
3n O n n n n =≈-+ 4) 可以估计出总的加减法运算的计算量也约为3
3
n
3. 列主元Gauss 消去法
前节所讲顺序Gauss 消去法，有以下两个问题：
1) 在消元过程中当某个元素对角线为零时，算法无法进行； 2) 如果某个对角线元素的绝对值非常小，将会引起很大的误差。

因此顺序Gauss 消去法不具实用性，但由此方法改进的列主元消去法是最为广泛使用的方法，以下我们通过一个实例体会一下顺序Gauss 消去法的问题和稍作改进之后的效果：
例3.2 用顺序Gauss 消去法求解线性方程组(运算过程中保留4位有效数字)
⎩⎨
⎧=-=+27
.3726.7453.415
.8713.87002.02121x x x x 解：写出方程组的增广矩阵
()⎪⎪⎭
⎫
⎝⎛-=27.3726.7453.415.8713.87002.0b A
消元过程：22275.2226002
.0453
.421≈==
m ，第2行减去第1行的21m 倍，得 ⎪⎪⎭
⎫ ⎝⎛--194100194000015.8713
.87002.0 回代过程得到计算解为：001.11940001941002==
x 35002
.0001
.1*13.8715.871-=-=x ，此方程
的准确解为1,1021==x x ，显然此计算解已经没有任何意义。

如果将两个方程交换一下位置再求解，增广矩阵：
⎪⎪⎭
⎫ ⎝⎛15.8713.87002.072.377.26-4.453 消元过程：0004491.0453
.4002
.021==
m ，第2行减去第1行的21m 倍，得 ⎪⎪⎭
⎫ ⎝⎛13.8713.87072.377.26-4.453 回代过程得到计算解为：00.10453
.4000
.126.727.37,000.113.8713.8712=⨯+===x x ，与精确解完全相同！
第1种解法产生非常大的误差的原因是11a 的绝对值太小，做为分母将产生太大的误差
而导致最后的解面目全非；用第2种方法是通过交换方程的位置使做分母的11a 的绝对值尽可能的大，正如第一章误差分析所述，避免绝对值太小的数做分母，可以提高计算结果的精度。

这一过程称为选主元(所选的对角线元素称为主元)，一般的做法如下。

设进行第k 步计算时增广矩阵如下：
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛-------n nn
nk
k kn kk k n
k k k k k n k k n k k b a a b a a b a a a b a a a a b a a a a a
0000001,1,11,12221,2221111,11211
为了在第k 次消元过程中使消元因子中的分母绝对值尽可能的大，在第k 列中取第k 到第n 个元素中绝对值最大者：
ik n
i k rk a a ≤≤=max
(3.8)
在选定主元之后，考虑如下情况
(1) 如果0=rk a ，则线性方程组的行列式的值为0，其解不定(无解或有无穷多组解，
停止计算)；
(2) 如果k r =，则kk a 的绝对值最大，以kk a 作为主元；
(3) 如果k r >则交换矩阵的第k r ,行，交换之后第k 行对角线元素绝对值变为最大。

而交换增广矩阵的两行相当于交换两个方程的位置，显然不改变原方程的解；
(4) 另一个需要指出的是，由于误差的原因判断0=rk a 的方法采取：对于指定的充分
小的0>ε，如果ε<rk a ，则认为0=rk a 。

增加选主元机制的Gauss 消去法即为列主元Gauss 消去法，列主元Gauss 消去法的计算过程包括：选主元，交换主元行，消元。

以下我们写出适合编程的Gauss 算法。

算法3.3 (列主元Gauss 消去法)
0) [初始化] 置系数矩阵A ，右端项向量b ，取充分小的数0>ε； 1) [消元过程] 对于1,,2,1-=n k 循环
1.1) [选主元] },,,max{n k i a a ik rk == 1.2) 如果ε<rk a ，停止计算，算法失败 1.3) 如果k r >交换k r ,方程
.,,,n k j a a rj kj =⇔， r k b b ⇔
1.4) [消元] 对于n k i ,,1 +=循环
1.4.1) kk ik ik a a a /= 1.4.2) k ik i i b a b b -=
1.4.3) .,,1,n k j a a a a kj ik ij ij +=-=
2) [回代过程] 执行算法3.1，即计算
1,2,,1,,/)(1
-=-
=∑+=n n k a x a
b x kk n
k j j kj
k k
3) [算法结束]
关于列主元Gauss 消去法的计算量，从列主元Gauss 消去法的算法可以看出，仅比顺序Gauss 消去法多找主元部分，这一部分的操作主要是比较运算及必要时交换两个方程，其他的部分算法完全一样，因此影响不大。

4. 全主元Gauss 消去法
全主元Gauss 消去法与列主元Gauss 消去法的区别在于取主元的方法是
},|max{|||n j i k a a ij rq ≤≤=
如果k r >交换增广矩阵的第k r ,行，如果k q >交换增广矩阵的第k q ,列(本质上是交换未知量q x 与r x ，计算结束时应该交换回来)。

显然全主元消去法取主元的方法范围更大，主元的绝对值会更大，计算过程也就更可靠更稳定，但计算量稍大。

实际计算效果表明，列主元Gauss 消去法已经很稳定，全主元法的改进甚微，是没必要的，所以通常使用列主元Gauss 消去法。

5. Gauss-Jordan 消去法
列主元Gauss 消去法在消元过程中是将对角线以下的元素约化为0，将系数矩阵化为上三角矩阵，通过回代得到方程的解；Gauss-Jordan 消去法是将系数矩阵约化成单位矩阵，无需回代就直接得到了方程组的解，也称为无回代过程的消去法。

设线性方程组的增广矩阵为：
()⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛=n nn
n n n n b a a a b a a a b a a a b A
21
222221
111211
第1步计算，设011≠a ，第1行除以11a ，第)1(≠i i 行减去新第1行的1i a 倍，则增广
矩阵变为如下形式
⎪⎪⎪⎪
⎪⎭
⎫
⎝
⎛n nn
n n
n b a a
b a a b a a 222221112
001
第2步计算，设022≠a ，第2行除以22a ，第)2(≠i i 行减去新第2行的2i a 倍，则增
广矩阵变为如下形式
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝
⎛n nn
n n n n
b a a b a a b a a b a a
333332223111300001001
设第k -1步计算结束时增广矩阵为如下形式
⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫
⎝
⎛---n nn
nk k kn kk k n
k k k n k b a a b a a b a a b a a
1,1,11111 1
第k 步计算，设0≠kk a ，第k 行除以kk a ，第≠i i (k )行减去新第k 行的ik a 倍，则增广
矩阵变为
⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫
⎝
⎛++++++n nn
nk k n k k k k n
k k k n k b a a b a a b a a b a a
1,11,1,1,111
,11 1 完成n -1步计算之后，增广矩阵变为
⎪⎪
⎪
⎪
⎪
⎭
⎫
⎝
⎛n b b b 100010001
21
显然最后一列就是方程组的解。

注意到在求解过程中，总是假设对角线上的主元0≠kk a ，若0=kk a 算法将无法进行
下去，如同列主元消去法，在消元过程中也可用选列主元的策略来避免这种情况发生，除非系数矩阵奇异。

算法3.4 (求解线性方程组的Gauss-Jordan 消去法)
0) [初始化] 置系数矩阵A ，右端项向量b ，取充分小的数0>ε； 1) [消元过程] 对于n k ,,2,1 =循环
1.1) [选主元] },,,max{n k i a a ik rk == 1.2) 如果ε<rk a ，停止计算，算法失败！ 1.3) 如果k r >交换k r ,方程
.,,,n k j a a rj kj =⇔ r k b b ⇔
1.4) [规格化主元行] 1.,,1,/=+==kk kk kj kj a n k j a a a
1.5) [消元] 即对于k i ≠循环
1.5.1) k ik i i b a b b -=
1.5.2) n k j a a a a kj ik ij ij ,,, =-=
2) [算法结束]
6. Gauss-Jordan 消去法求逆矩阵及算法设计 首先通过一个求逆过程来写出相应算法。

例3.3 用Gauss-Jordan 消去法求下矩阵的逆阵：
⎪⎪⎪⎭
⎫ ⎝⎛=641321110A
解：在矩阵A 的右边添加单位矩阵进行初等行变换
()⎪⎪⎪
⎭
⎫ ⎝⎛=100641010321001110I A
3
2213
2110320001110010321100641001110010321100641010321001110r
r r r ⇔⇔→⎪⎪⎪⎭
⎫ ⎝⎛-→⎪⎪⎪⎭⎫ ⎝⎛→⎪⎪⎪⎭⎫ ⎝⎛消元
⎪⎪⎪⎭⎫
⎝⎛----→⎪⎪⎪⎭
⎫ ⎝⎛---→⎪⎪⎪⎭⎫ ⎝⎛-112100113010
1200011-0001012
0001001110110320010321212
12
121
2
123消元消元 此时对应初始状态的单位矩阵就是所求的逆矩阵：
⎪⎪⎪⎭
⎫ ⎝⎛----=-112113120
1
A
通常编程实现求逆矩阵时并不添加右边的单位矩阵，而是对原始矩阵直接进行变换，在完成计算时原矩阵就变成了逆矩阵。

其实我们注意一下求逆的过程，有如下两个特征：
首先、从初始状态到最终状态，增广矩阵的每个状态都有n 个列是单位列向量，这也
就意味着每次变换，都会将原矩阵的一列变为单位列向量，而将原单位矩阵的一个列向量变为逆矩阵的一列，每次变换时可让这一对列向量占用共同的列，也就是可将右边新产生的一列直接放在左边将要变为单位列向量的位置上。

其次、在求逆过程中如果没有行交换，则增广矩阵中单位列向量的消失与新增，在左右两部分都是按自然顺序递增的。

从例3.3可以看出，因交换主元行，右边的单位列向量顺序随主元行而变，因此在完成计算后需要换回原位置，这正说明原矩阵的行交换对应逆矩阵的列交换这一性质。

以下过程是通常编程时的实际计算过程：
⎪⎪⎪⎭
⎫ ⎝⎛---→⎪⎪⎪⎭⎫ ⎝⎛----→⎪⎪⎪⎭⎫ ⎝⎛---→⎪⎪⎪⎭⎫
⎝
⎛-→⎪⎪⎪⎭⎫ ⎝⎛-→⎪⎪⎪⎭⎫ ⎝⎛→⎪⎪⎪⎭⎫ ⎝⎛⇔⇔⇔⇔11
21-131********
012-01211032132
13211103216411103216413211102132322111131
132c c c c r r r r 消元消元消元
整个计算有3=n 个消元过程，但是最后一次消元不须选主元，所以主元行标分别是{2,3}，按逆序分别交换列21,32⇔⇔，则得到正确的逆矩阵，此过程直观性太差，仅作为调试程序时核对之用。

算法3.5(矩阵求逆)
0) [初始化] 输入矩阵A ，并以向量p 记录主元行下标(共n -1个主元行下标) ，取充分
小的数0>ε；
1) [行变换] 对n k ,,2,1 =循环
1.1) [选主元] },,,max{n k i a a ik rk ==，r p k =
1.2) 如果ε<rk a ，矩阵不可逆，算法失败；否则计算rk a /1=σ
1.3) [交换主元行] 如果k r >交换k r ,行 n j a a kj rj ,,2,1, =⇔ 1.4) [规格化主元行] k j a a kj kj ≠=,σ，σ=kk a 1.5) [消元] 对n k k i ,1,1,,1+-=循环
1.5.1) n k k j a a a a kj ik ij ij ,,,1,,1, +-=-=
1.5.2) ik ik a a σ-=
2) [交换列] 对1,2,,1 -=n k 循环
如果k p k ≠，交换k p k ,列 n i a a k ip ik ,,2,1, =⇔
3) [算法结束] 以下是计算机求逆验证矩阵
⎪⎪
⎪⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎭
⎫
⎝⎛++-+--
--+-
++=2222
12
2121112
1211212212
1222n n n n n n A
(3.9)
()
j i n A ij
--=-1
，n j i ,,2,1 =，
(3.10)
3.2 LU 分解
1. 消去法的矩阵运算表示形式
从矩阵计算的角度讲Gauss 消去法是通过一系列初等行变换将矩阵约化成上三角矩阵；
Guass-Jordan 消去法是通过一系列初等行变换将矩阵约化成单位矩阵。

矩阵的初等行变换本质上是矩阵左乘一个初等矩阵，因此消去法也可以用矩阵乘积形式进行描述。

Gauss 消去法的矩阵形式：
⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛--+++++++--+----++n nn
k n nk
k n
k k k k k k kn k k kk k n k k k k k k k n k k n nk k k b a a a b a a a b a a a b a a a a b a a a a a m m
1
,1,11,1,11,1,11,1,11,1111,11111
,10
0000011111
⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫ ⎝
⎛=++++++--+----+n nn
k n k n
k k k k kn k k kk k n k k k k k k k n k k n b a a b a a b a a a b a a a a b a a a a a
1
,1,11,11,1,11,1,11,1111,111110
00000000 Gauss-Jordan 消去法的矩阵形式：
⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛----+++++++--+--++-n nn
k n nk
k n
k k k k
k k kn k k kk k n k k k k k n k k nk k k kk k k k b a a a b a a a b a a a b a a a b a a a m m a m m
1
,1,11,1,11,1,11,1,1111,11,1,1,10
0000101
11/111
⎪⎪⎪⎪⎪⎪⎪⎪
⎪⎪⎭
⎫ ⎝
⎛=++++++--+-+n nn
k n k n k k k k kn k k k n k k k n k b a a b a a b a a b a a b a a
1
,1,11,11,1,11,1111,10
00000100010001
其中(3.9)式和(3.10)式中的kk ik ik a a m /=就是我们前述消去法的消元因子，如果记：
1,,2,1,1111
,1-=⎪⎪
⎪⎪⎪⎪
⎪⎪⎭
⎫
⎝⎛=
+n k l l L nk k k k
(3.11)
则Gauss 消去法的整个过程可以表示为：
b L L L Ax L L L Ux n n 121121 --==
注意到形如(3.11)式的单位下三角矩阵具有以下两个性质：逆阵仍然是单位下三角阵；任意两个单位下三角阵的乘积仍然是单位下三角阵，所以上式也可以表示为：
b L Ax L Ux 11--==
因此Gauss 消去法的本质上是将系数矩阵A 进行如下的矩阵分解：
LU A = (3.12)
其中L 是单位下三角矩阵，U 是上三角矩阵，此种分解也称为杜里特尔(Doolittle)分解，如果L 是下三角矩阵，U 是单位上三角矩阵，则称为克洛特(Crout)分解，以下内容讨论的是杜里特尔分解方法。

2. LU 分解
上述(3.12)式也可以写成：
⎪⎪⎪
⎪
⎪
⎭
⎫
⎝
⎛⎪⎪⎪⎪⎪⎭⎫
⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛nn n n n n nn n n n n u u u u u u l l l a a a a a a a a a 222
112112
1
21
2
1
2222111211
11
1
实现矩阵的LU 分解算法既可以仿照消去法设计算法也可以直接从该式推导出计算公
式，以下我们讨论直接LU 分解方法。

首先用L 的第1行乘以U 的各列，则有：
n j u a j j ,,2,1,111 =⨯=，
由此可以得到U 的第1行：
n j a u j j ,,2,1,11 ==
再让L 的各行乘以U 的第1列得到n i u l a i i .,2,1111 ==，由此又可以得到L 的第1列
n i u a l i i ,,3,2,/1111 ==
假如我们已经得到了L 的前k -1列和U 的前k -1行：
⎪⎪
⎪⎪⎪⎪
⎪⎪
⎭⎫
⎝
⎛⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛=--------?????
1???111
,1,11,1111,1111
,1
1,11
,1
n k k
k k k n k k k n n k k k k u u u u u u u l l l l l A
按如下方法确定U 的第k 行和L 的第k 列，将L 的第k 行与U 的第)(k j j ≥各列相乘，有
n k j u u l u l a kj j k k k j k kj ,,,,11,11 =+++=--
从该式可以得到U 的第k 行
n k j u l u l a u j k k k j k kj kj ,,),(,11,11 =++-=--
将L 的第i (i>k )各行与U 的第k 列相乘，有
n k i u l u l u l a kk ik k k k i k i ik ,,1,,11,11 +=+++=--
从该式可以得到L 的第k 列
{}n k i u u l u l a l kk k k k i k i ik ik ,,1,/)(,11,11 +=++-=--
则有如下结果：
⎪⎪
⎪⎪⎪⎪
⎪⎪
⎭⎫
⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫
⎝⎛=--------??1?111,1,11,1111
,1111,11,11,1
kn kk n k k k k k n k k nk k n n k k k k u u u u u u u u u l l l l l l A
按如此方法，从1=k 执行到n -1就完成了矩阵的LU 分解。

综上所述，我们得到直接LU 分解的计算公式：
n k n
k i u u l a l n
k j u l a u k q kk qk iq ik ik k q qj kq kj kj ,...,2,1,,...,1/)(,...,1
11
1=⎪⎪⎩
⎪⎪⎨⎧
+=-==-=∑∑-=-=
(3.13)
定理3.2 设n 阶矩阵A 的顺序主子式均不为0，则A 存在LU 分解，且分解是唯一的。

事实上，如果分解成功，A 的第k 个顺序主子式的值就是kk u u u 2211，因此结论是显然的。

本定理表明LU 分解可以正常进行的前提是顺序主子式不为0，表现在分解过程中就是所有n k u kk ,,2,1,0 =≠。

通常可以将下三角矩阵L 和上三角矩阵U 放在一个矩阵内，即：
⎪⎪
⎪
⎪
⎪
⎭
⎫
⎝⎛nn n n n n u l l u u l u u u LU 2
1
22221
11211
：
(3.14)
例3.4 求矩阵A 的LU 分解
⎪⎪⎪⎭
⎫ ⎝⎛=516234212A
解：设
⎪⎪⎪⎭
⎫
⎝
⎛⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛332322
13121132
31
2111
1
516234212u u u u u u l l l 对于k =1，利用(3.13)式，先计算U 的第1行
2,1,2131211===u u u
再根据3121,a a 的计算公式计算L 的第1列 112121u l a =， 22/4/112121===u a l
113131u l a =，
32/6/113131===u a l
当前的分解为
⎪⎪⎪⎭
⎫
⎝⎛⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛3323223221
213121
516234212u u u l
对于k =2，根据2322,a a 的计算公式计算U 的第2行 22122122u u l a +=， 112312212222=⨯-=-=u l a u
23132123u u l a +=， 222213213223-=⨯-=-=u l a u
再根据32a 的计算公式计算L 的第2列
2232123132u l u l a +=，
()()21/131/2212313232-=⨯-=-=u u l a l
当前的分解为
⎪⎪⎪⎭⎫
⎝
⎛-⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛3321212123121516234212u 对于k =3，根据33a 的计算公式计算33u
332332133133u u l u l a ++=
5)2()2(235233213313333-=-⨯--⨯-=--=u l u l a u
最终得到分解为：
⎪⎪⎪⎭
⎫
⎝⎛--⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛521212123121516234212 将下三角矩阵L 和上三角矩阵U 放在一个矩阵内，如：
⎪⎪⎪⎭
⎫ ⎝⎛---52321221
2 3． 用LU 分解求解线性方程组 考虑线性方程组
b Ax =
如果系数矩阵A 存在LU 分解A=LU ，代入原方程组
b Ux L =)(
令Ux y =，则可先求解
b Ly =
然后再求解
y Ux =
即可得到原线性方程组b Ax =的解。

此处两个方程组，一个是下三角方程，一个是上三角方程。

上三角方程y Ux =的求解可以使用算法3.1，现在考虑下三角方程的求解。

将下三角方程展开
n
n n n n n b y y l y l b y y l b y =+++=+=--11,112
21211
1...
求解上三角线性方程组是按11,,,x x x n n -的次序求解所以称为回代过程，求解此下三角线性方程组时次序正好相反是n y y y ,,,21 所以称为顺代过程，计算公式如下：
n k y l b y k j j kj k k ,...,2,1,
1
1
=-=∑-=
算法3.6 (求解线性方程组的直接LU 分解算法) 0) [初始化] 输入系数矩阵A ，右端项b 1) [LU 分解] 对于n k ,,2,1 =循环
1.1) [计算U 的第k 行]
n k j u l a u k q qj kq kj kj ,,,
1
1
=-=∑-=
1.2) [计算L 的第k 列]
n k i u u l a l k q kk qk iq ik ik ,,1,
/)(1
1
+=-=∑-=
2) [解下三角方程]
n k y l b y k j j kj k k ,...,2,1,
1
1
=-=∑-=
3) [解上三角方程]
1,,1,,
/)(1
-=-
=∑+=n n k a x a
y x kk n
k j j kj
k k
4) [算法结束]
例3.5 用LU 分解方法求解下列线性方程组
⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛1033516234212321x x x
解：先对系数矩阵LU 分解得到：
⎪⎪⎪⎭
⎫
⎝⎛--=⎪⎪⎪⎭⎫ ⎝⎛-=521212,123121U L 解下三角方称：
5
33,1033123121321321-=-==⎪⎪⎪⎭
⎫
⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛-y y y y y y
解上三角方程：
1
11,533521212321321=-==⎪
⎪⎪
⎭
⎫ ⎝⎛--=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--x x x x x x 前述的LU 分解相当于顺序Gauss 消去法，如同顺序Gauss 消去法不具实用性一样，算法3.6也不具实用性，应该使用更稳定的列主元策略。

由线性代数知识可知，交换两行相当于系数矩阵左乘一个交换矩阵，将所有交换矩阵的乘积记为P ，就得到如下定理：
定理3.3 如果矩阵A 非奇异(矩阵可逆)，则存在交换矩阵P ，单位下三角阵L 和上三角阵U ，使得 LU PA =
(3.15)
用列主元实现矩阵的直接LU 分解远不如对算法3.3稍作改进更容易，以下我们给出这
一改进的算法。

算法3.7 (列主元LU 分解算法)
0) [初始化] 设置系数矩阵A ，向量},,2,1{n p =存放主元行下标，取充分小的数
0>ε；
1) [分解过程] 对于n k ,,2,1 =循环
1.1) [选主元] },,,max{n k i a a ik rk == 1.2) 如果ε<rk a ，停止计算，算法失败！ 1.3) 如果k r >交换第k r ,行
r k rj kj p p n j a a ⇔=⇔，,,2,1,
1.4) 对于n k i ,,1 +=循环
1.4.1) kk ik ik a a a /=
1.4.2) .,,1,n k j a a a a kj ik ij ij +=-=
2) [算法结束] 本算法的几点说明：
1) 向量p 记录主元行的下标，必要时可以生成交换矩阵P ，完成矩阵或方程的交换功能，在求解方程时根据主元行下标对右端项作相应的交换；
2) 交换主元行时，要对n 个元素进行交换，这是因为此时矩阵L 也要做相应的交换；
3) 每次循环都对k 行k 列以后的元素进行更新，步骤1.4.1)是对ik l 的最后计算结果，对kj u 的计算则隐含在步骤1.4.2)中；
4) 完成算法之后，下三角矩阵L 除主对角上的1之外其他元素置于矩阵A 对应的位置上，上三对角矩阵U 置于矩阵A 的对应位置上。

3.3 三对角方程组的追赶法
本节考虑如下形式的线性方程组
⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝
⎛-----n n n n n n
n n n d d d d d x x x x x b a c b a c b a c b a c b 132********
333
22211
(3.16)
称上述方程组为三对角方程组，其系数矩阵称为三对角矩阵。

其矩阵形式为：
n n n R d x R A d Ax ∈∈=⨯,,,
(3.17)
其中系数矩阵A 是三对角矩阵。

在一些实际问题中，例如解常微分方程边值问题、偏微分方程的差分方程、解热传导方程、三次样条函数、数值微分等，都会遇到此类线性方程组，求解此类方程组，用Gauss 消去法显然会浪费内存和计算量。

以下我们通过LU 分解，导出求解三对角方程更为简便有效的算法，称其为追赶法。

对三对角矩阵进行求解LU 分解：
⎪⎪
⎪
⎪⎪⎪⎭
⎫ ⎝
⎛⎪⎪⎪⎪⎪⎪⎭⎫
⎝
⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝
⎛---n n n
n n
n n u c u c u c u l l l b a c b a c b a c b 13
2
2
113211
322211
11
1
1
右端两矩阵相乘我们可以得到：
n i b u c l a u l b u i i i i i
i i ,,3,2,1111 =⎭
⎬
⎫
=+==-- 从而得到以下三对角阵LU 分解的递推公式：
n
i c l b u u a l b u i i i i i i i ,...,3,2,/11
1
1=⎭
⎬⎫
-===-- (3.18)
在得到三对角阵的LU 分解之后，先解下三角方程：
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝
⎛--n n n n n
d d d d y y y y l l l 1211213211
11
由此得到y 的计算公式：
n
i y l d y d y i i i i ,...,3,211
1=-==-
再求解上三角方程：
⎪⎪
⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝
⎛----n n n n n n n y y y y x x x x u c u c u c u 12112111
22
11
由此得到x 的计算公式： 1
,...,1/)(/1-=-==+n i u x c y x u y x i
i i i i n
n n
定理3.4设三对角线性方程组d Ax =中系数矩阵A 满足： 1) 11c b > 2) n n a b >
3) 1,,3,2,0,-=≠+≥n i c a c a b i i i i i
则系数矩阵A 非奇异，且有：
1) n i u i ,,2,1,0 =≠
2) 1,,2,1,1/0-=<<n i u c i i
3) 1,,3,2,-=+<<-n i a b u a b i i i i i
证明：注意到结论1)与系数矩阵A 非奇异是等价的。

下面我们用归纳法来证明结论1)和结论2)。

对于1=i ，11b u =，显然结论成立；先假设结论对1-i 成立，则对i 有
011
1
>≥-≥-≥-=---i i i i i i
i i i i i c a b c u a b c l b u
从而可以得到结论1)和结论2)均成立。

考虑结论3)，根据结论2)和i u 的计算公式有：
i i i i i i i i i i i i i i i i a b c u a
b c l b u c u a b a b -<+≤-=≤-
<------11
111
算法3.8 (求解三对角方程的追赶法)
0) [初始化] 输入数据d c b a ,,,，令11b u =。

1) [矩阵分解]
n i c l b u u a l i i i i i i i ,,3,2,,/11 =-==--
2) [解下三角方程]
⎩⎨
⎧=-==-n
i y l d y d y i i i i ,,3,211
1 3) [解上三角方程]
()⎩⎨
⎧-=-==+1
,,1//1 n i u x c y x u y x i
i i i i n n n
4) [算法结束]
例3.6 求解三对角方程(计算过程保留小数点以后3位数字)
3
41241
242244
3
43
2
32
1
21-=+=++-=++=+x x x x x x x x x x
解：首先对三对角矩阵进行LU 分解，结果为：
⎪⎪⎪⎪
⎪
⎭
⎫
⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛417.32
3.42923.5002000.410.29210.28610.25014124124124 先解下三角矩阵
417
.3429.150.1000.2,31
1210.2921286.010.250132214321-==-==⎪⎪⎪⎪
⎪
⎭
⎫
⎝⎛--=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛y y y y y y y y
再解上三角矩阵
1
11
1417.3429.1500.1000.1417.32429..32500.32000.443214321-==-==⎪⎪⎪
⎪
⎪
⎭
⎫
⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛x x x x x x x x ，
用追赶法求解三对角方程组的计算量，也就是乘除法运算次数仅为)(n O ，所以对于三对角方程组通常不使用Guass 消去法而是使用追赶法也就很容易理解了。

有时还会遇到如下形式称为广义三对角方程组(如第4章的三次样条中第三类边界条件形成的三弯矩方程)：
⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-----n n n n n n
n
n n n d d d d d x x x x x b a c c b a c b a c b a a c b 132********
333
222111
(3.19)
对(3.19)式也可以设计出)(n O 计算量的算法，仿照Gauss 消去法，可用如下方法求解。

首先消元第一列：显然仅需要消元n c a 、2，在完成对应的两个方程的消元之后会增加
)2,(),2(n n 、两个位置的非零元素，我们仍然用n c a 、2表示，成为如下形式：
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡-----n n n n n n
n
n n n d d d d d x x x x x b a c c b a c b a a c b a c b 1321132111
133
3
222111
一般情况，对于3,,2,1-=n k ，每列仅两个位置),(),1(k n k k 、+的元素n k c a 、1+应该消元，同时处在)1,(),1(++k n n k 、的位置上会新增两个非零元素，我们仍记为本次被消元的元素符号n k c a 、1+，当完成了3-=n k 的消元过程时是如下形式
1111122
22222
22211111n n n n n n n n n n n
n
n n n b c a x d b c a x d b c a x d a b c x d c a b x d ----------⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦⎣⎦⎣⎦
(3.20)
显然只要再完成(3.20)式对n n n a c a 、、1-的消元就可以将其变为如下形式
11
11122
2222
2
2221
111n n n n n n n n n n n n b c a x d b c a x d b c a x d b c x d b x d ---------⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦⎣⎦⎣⎦
(3.21)
剩下的问题就是一个简单的回代算法。

算法3.8A 求解广义三对角方程算法 0) [初始化] 输入数据d c b a ,,, 1) 对于3,,2,1-=n k 执行
k k k k k c b a b b 111
+++-=，k k k k k d b a d d 111+++-=，k k
k k a b a a 1
1++-=
k k n n n a b c b b -
=，k k n n n d b c d d -=，k k
n n c b c c -= 2) 对于2-=n k 做
22111------
=n n n n n c b a b b ，22111------=n n n n n d b a d d ，22
111------=n n n n n a b a
c c 22---
=n n n n n c b c a a ，22---=n n n n n a b c b b ，22---=n n n n n d b c
d d 3) 对于1-=n k
11---
=n n n n n c b a b b ，11
---=n n n n n d b a
d d 4) [回代过程]
n
n
n b d x =
1
111-----=
n n
n n n b x c d x
1,2,,1,21 --=--=
+n n k b x a x c d x k
n
k k k k k ，
5) [算法结束]
显然容易算出该算法的计算量约为n 11次乘除法运算。

3.4 平方根法
1. 平方根法
本节我们讨论对称正定矩阵的一个更为有效的分解方法—平方根法或Cholesky 分解，首先回顾线性代数中对称正定矩阵的有关结论。

定义3.1 设矩阵⨯
∈nn R
A 满足A A T
=，则称矩阵A 是对称矩阵，如果对任意非零向量
n R x ∈均有
0>Ax x T
则称A 是正定矩阵，如果有0≥Ax x T
则称A 是半正定矩阵。

定理3.5 设A 是对称正定矩阵，则有以下结论成立： 1) 所有特征值为正； 2) 所有顺序主子式为正；
3) 如果对称正定矩阵可逆，则逆矩阵也是对称正定矩阵。

定理3.6 设A A R
A T n
n =∈⨯,，且A 的所有顺序主子式均大于0，则存在唯一的单位
下三角阵L 和对角阵D ，使得
T LDL A =
(3.22)
证明：因为A 的顺序主子式不为零，所以A 存在唯一的LU 分解A=LU ，注意到上三角矩阵U 可以分解为：
DV u u u u u u u u u u u u u u u U n
n nn nn n n ≡⎪⎪⎪⎪⎪
⎪⎪⎭
⎫
⎝
⎛
⎪⎪⎪⎪⎪⎭⎫
⎝
⎛=
⎪⎪⎪⎪⎪⎭⎫
⎝⎛=11122
2111111222
11
22211211
将其带入LU A =则有LDV LU A ==，由对称性得
()()LU DL V LDV A A T
T T
T ====
再根据LU 分解的唯一性，则有T L V =得证。

注意到如果矩阵A 对称正定的，则有n i u ii ,,2,1,0 =>，令：
),,,(22112
1nn u u u diag D =
将(3.22)式可以改写成如下形式：
T T
T
T L L LD LD L D LD LDL A 11212
12
12
1))((====
则有如下结论。

定理3.7 (Cholesky 分解) 设A 是正定对称矩阵，则存在唯一的非奇异下三角阵L ，使得
T LL A =
(3.23)
其中L 的对角元素均为正数，称为Cholesky(乔莱斯基)分解或平方根分解。

至此我们仅给出了Cholesky 分解的存在性条件，以下讨论分解算法。

如同LU 分解，
也考虑直接的分解方法，从下式推导出计算公式：
⎪⎪⎪
⎪
⎪
⎭
⎫
⎝⎛⎪⎪⎪⎪⎪⎭⎫
⎝⎛=
⎪⎪⎪⎪⎪⎭⎫ ⎝⎛nn n n nn n n nn n n n n l l l l l l l l l l l l a a a a a a a a a
222121112
1
22
21
11
2
1
2222111211
(3.24)
我们写出下三角部分的表达式为
n k k i l l l l l l a kk ik k j kj ij n j kj ij ik ,,1,,1
1
1
+=+==∑∑-==
假设已得到L 的第1列到第k -1列，现要计算第k 列，则有
n k n
k i l l l a l l a l kk k j kj ij ik ik k j kj
kk kk ,,2,1,,,1,/11112
=⎪⎪⎩
⎪
⎪⎨⎧+=⎪⎪
⎭⎫ ⎝⎛-=-=∑∑-=-= (3.25)
例3.7 求以下对称矩阵的Cholesky 分解，并求其行列式的值
⎪⎪⎪⎭
⎫ ⎝⎛=653553333A
解：33
3,333,31131311121211111======
==
l a l l a l a l
()
()22
335,23522
213132322
212222=
⨯-
=
-=
=-=-=l l l a l l a l
()()
12362
2
2322313333=--
=--=l l a l
⎪⎪⎪⎪⎭
⎫ ⎝
⎛=12323
3
L ()()61
232
22
==
A
2．改进平方根法
平方根算法是稳定的，缺点是算法中有开平方运算。

一次开平方运算的CPU 时间远远
大于乘除运算(开平方运算一般是迭代法实现的)，因此通常我们应该尽量避免使用平方根运算，改进平方根法就是无开方运算的平方根法。

其实，所谓无开方运算的平方根法就是基于上节所述的(3.22)式进行求解，将(3.22)展开
⎪⎪⎪⎪⎪
⎭
⎫
⎝⎛⎪⎪⎪⎪⎪⎭⎫
⎝
⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=111111
21212
12
1
21
n n n n n l l l d d d l l l A 以LD T =作为中间矩阵，则有j ij ij d l t =，直接对T
TL A =按矩阵乘法规则推导出计算公式
为：
n k l t a d k i d t l l t a t a d k j kj kj kk k i ik ik i j ij kj ik ik ,,3,2,1,,2,1,/111
111
1 =⎪⎪
⎪
⎭
⎪⎪⎪⎬⎫-=-=⎪⎭
⎪⎬⎫
=-==∑∑-=-= (3.26)
如果求解线性方程组Ax=b ，以下只要分别求解下三角方程b Ly =和上三角方程
y D x L T 1-=，即用如下计算公式：
1
,,1,/,,2,11
1
1 -=-
==-=∑∑+=-=n n k x l
d y x n
k y l b y n
k i i
ik
k k k k j j
kj k k
(3.27)
在编程实现改进平方根算法时，对角阵D 可以存放在矩阵L 的对角线上，中间矩阵T 可以借用矩阵L 的上三角部分，改进平方根算法如下：
算法3.9 改进平方根矩阵分解算法
0) [初始化] 输入系数矩阵和右端项1111a l b A =，， 1) [分解] 对于n k ,,3,2 =循环
1.1) [计算L ] 对1,,2,1-=k i
ii
ik ki i j jk
ij ik ik l l l l l a l /1
1
=-=∑-=
1.2) [计算D ]
∑-=-=1
1
k j jk kj kk kk l l a l
2) [算法结束]
注：(1) 在编程时，如果希望将下三角矩阵L 和中间矩阵T 存放在原矩阵A 内，只需要将算法3.9中的所有ij l 改为ij a 即可，其它不须作任何修改；(2) 如果利用改进平方根算法求解线性方程组，可以在算法3.9的基础上按照(3.27)公式继续计算即可以得到方程组的解。

例3.8 用改进平方根算法分解以下矩阵
⎪⎪⎪⎪
⎪⎭
⎫
⎝
⎛=1097
4997477744444
A 解：按算法3.9步骤1)可求得如下结果：
⎪⎪⎪⎪⎪⎭
⎫
⎝
⎛=11
1122113331
4444)\(T L 即：
⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫
⎝
⎛=12
3411
11011100110001
D L ，
3.5 迭代法
以上所述Guass 消去法、LU 分解求解线性方程组的方法均属于直接方法，还有一类求解线性方程组的方法属于迭代法，迭代法在求解大型稀疏问题、病态方程组等特殊情况时是一个非常有效的方法。

求解线性方程组的直接方法和迭代法的主要区别是：直接法可以用有限次初等运算求出其解，比如Gauss 消去法可以用大约3/3n 次乘除法运算和3/3n 次加减法运算求出方程组的解；而迭代法思想如同非线性方程的迭代法一样，先给出一个初始近似解，根据一个迭代公式得到一个更好的近似解，通过反复执行迭代过程逐步逼近到精确解，当近似解达到指定的精度时结束计算。

如同求解非线性方程的迭代法，求解线性方程组迭代法的也是先把所给的线性方程组
b Ax =，化成同解的线性方程组
f Bx x += (3.28)
这里称B 为迭代矩阵，给定初始向量)0(x 并按迭代格式
,2,1,0,)()1(=+=+k f Bx x k k
(3.29)
进行迭代计算。

如果按上述迭代格式所得到的向量序列{})
(k x
收敛于某一向量*
x
，则*
x 就是方程组
b Ax =的解，并称此迭代法收敛；否则称为不收敛或发散。

本课程仅介绍Jacobi 迭代法和Gauss-Seidel 迭代法，并简单讨论其收敛条件。

1. Jacobi 迭代法。