研究生数值分析(8)直接三角分解法

§2 直接三角分解法
三角分解法也是直接法，基本思想是： 将系数矩阵A分解为两个三角形矩阵L和U的 乘积A=LU ，将方程组AX=b的求解问题归结为 两个三角形方程组 LY=b与UX=Y的求解问题。
即：先由LY=b求出Y ，然后由UX=Y求出X ， 从而获得AX=b的解。
(1) A为一般稠密（零元素占很小比例）矩阵 的杜利特尔(Doolittlr)和克劳特(Crout)分解法；
1 0 0
3

7

2
2 2

1 12
7 6

1

0
0
16

3
等价的三角方程组为


12


3 3 2
3 7
2

x1

x2


15

15

2

16


x3

16
an1
bn1
cn1


xn1


fn1

an bn xn fn
其中方程组AX=f 的系数矩阵A的元素满足条件：
b1 c1 0

bi ai ci

bn

an
0
（1） 且 aici 0
(i 2,3, , n 1)
由以上推导过程知，方程组AX=f 有唯一解
由（3）式可得计算 i 的递推公式
1i

c1 ci
/ b1 /(bi

ai i1)
(i 2,3,
（4）
, n 1)
只要计算出 i ，其它待定系数 i 与 i 均可
通过已知数 ai , bi 与 i 表示。
上述解方程组 AX=f 的过程归纳为：
杜利特尔Doolittle分解法
设方程组AX=b的系数矩阵的各阶顺序主子式
det( Ak ) 0 (k 1,2, ,n) ，则存在唯一杜利特尔分解 A LU ，其中
1 l21 1 L l31 l32 1

u11 u12
U
u22


u1n
u2n





ln1 ln2 ln3
1

unn

下面介绍直接根据A的元素计算L、U元素的分解方法 第一步求U的第一行元素和L的第一列元素 第二步求U的第二行元素和L的第二列元素......
对那些明确是1或是0的元素不再求。 由矩阵乘法规则与相等条件,
利用 aij 在上述计算过程中，
导出计算 lij 或 uij 的公式。
li2 (ai2 li1u12 ) / u22 (i 3, 4, , n)
（1） （2）
第k步计算U的第k行L的第k列元素的公式为：
k 1
uki aki lkju ji j 1
(i k, k 1, , n) （3）
k 1
lik (aik liju jk ) / ukk j 1
例如
第一步计算由 ai1 li1u11 得
u1i a1i (i 1,2, ,n)
第二步计算由 a1i u1i 得 li1 ai1 / u11 (i 2,3, ,n)
由 a2i l21u1i u2i 得
u2i a2i l21u1i
(i 2,3, ,n)
由 ai2 li1u12 li2u22 得
它与公式（3）相似。若将b作为增广矩阵的最
后一列元素,那么对增广矩阵作LU分解, b也作相应
运算,仍在最后一列,则分解后的最后一列即为y 。
于是
k 1
ukj aki lkju ji j 1
(i k, k 1, , n, n 1)
例3：将方程组
12x1 3x2 3x3 15 18x1 3x2 x3 15 x1 2x2 x3 6 的系数矩阵作LU分解，并求方程组的解。
根据系数矩阵A的特点，设
1 2 2
1 1


1 2

A

LU

n1 n1

1

n1


n n
1
其中 i , i , i 为待定系数。比较A与LU对应
的元素，有
b1 1, c1 11
ai i
(i 2,3,
bi i i1 i (i 2,3,
ci i i
(i 2,3,
由（1）和（2）可以看出
a1 b1 c1 0
, n) , n) （2）
, n 1)
因此有 1 c1 / a1且0 1 1 由 a2 b2 a21 b2 a2 1 b2 a2 c2 0 有 2 c2 / a2且0 2 1
, 2

2, 7
3

7 26
按递推公式（5）计算 yi (i 1, 2,3, 4)
得
y1


1, 4
y2

1, 14
y3


1 52
,
y4

1 90
按递推公式（6）计算 xi (i 4,3, 2,1) 得
x4

1 90
, x3


1, 45
x2

7, 90
x1

13 45
用追赶法解方程组 AX=f 仅需 5n-4 次乘除 法运算, 计算过程稳定。
x2


4

0

0
1 0
4 1
1 4

x3

x4
4 3
(2) A为三对角的追赶法。
1 Doolittle分解法和Crout 分解法 把一个n阶矩阵A分解成两个三角形矩阵相
乘的形式称为矩阵的三角分解。
矩阵三角分解的常见形式是：
A=LU，其中L为下三角阵，U为上三角阵。
作为特例，若L为单位下三角阵(对角元都是1 的下三角阵), U为上三角阵，则称为杜利特尔 (Doolittle)分解；
练习：1、利用LU分解求解四元线性方程组AX=b，
其中
1
A
2
2

1
2 1 2
5 3 2, 2 3 5
3
2
3

4
b 7

1
0

2、用追赶法求方程组的解
4 1 0 0 x1 3

1
4
1
0



3
用回代法解得
x3 3, x2 2, x1 1
4 追赶法求解三对角线性方程组 在样条函数的计算、微分方程数值求
解中常遇到如下形式的线性代数方程组：
b1 c1 a2 b2 c2
x1 f1

x2


f2


2

1
2

0

0
1
2
1 2 0
0 1 2 2
1
0
0
1 2 2

x1

x2

x3 x4



1 2
0
0
0

解：按递推公式（4）计算 1, 2 , 3 得
1

1 2
若U为单位上三角阵(对角元都是1的上三角阵),
L为单位下三角阵,则称为克劳特(Crout)分解。
下面分析实现矩阵杜利特尔(Doolittle)分解和克 劳特(Crout)分解的条件,讨论这些分解的唯一性。
定理2 (矩阵三角分解基本定理)
设 A Rnn 。若A的顺序主子式
det( Ak ) 0
解：增广矩阵为
12 3 3 15 18 3 1 15 1 2 1 6


12 3 3 15
LU分解的紧凑格式为
3
2
3
7

15

2 2 2

1 12
7 6
16 3

16

所以系数矩阵


的三角分解为

A
1 3
0 0 12 3 3


xk

( yk

n
ukj xj ) / ukk

j k 1
（6）
(k n 1,n 2, ,1)
杜利特尔矩阵分解 求解线性方程组的过程为: 10 实现A=LU分解，即
(a)按计算公式(1),(2)依次计算U的第1行元素
u1i (i 1, 2, , n) 与L的第1列元素 li1 (i 2,3, , n)
(k 1, 2, , n)
则存在唯一的杜利特尔分解A=LU，其中L
为单位下三角阵，U 为非奇异上三角阵。
还可以证明存在唯一的克劳特Crout分解
A LU
这里L为非奇异下三角阵，U为单位上三角阵。
如果A是一般非奇异阵，由列主消元法，A适 当行交换后，可使A的各阶顺序主子式 det( Ak ) 0， 从而实现杜利特尔Doolittle或克劳特Crout分解。
式(6)依次计算 xn , xn1, , x1
为便于记忆，我们给出L、U分解紧凑格式:
u11 u12

l21
u22



A


lk1
lk 2



li1 li2

ln1 ln2
u1k u2k
lk ,k 1
uk 1,k ukk lk 1,k
li,k 1
lik
ln,k 1
(b) 对k+2,3,…,n 按计算公式(3),(4)依次
计算U的第k行元素 uki (i k, k 1, , n) 与L的第
k列元素 lik (i k 1, , n; k n)
20 求解三角形方程组LY=b，即按计算公
式(5)依次计算 y1, y2 , , yn
30 求解三角形方程组UX=Y，即按计算公
(i k 1, , n; k n) （4）
在我们利用杜利特尔矩阵分解解线性方程 组AX=b时,只要实现矩阵分解A=LU,依次解三角 形方程组LY=b与UX=Y即可。
计算公式：


y1

yk
b1 bk

k 1
lkj
yj

j 1
(k 2,3,
（5）
, n)
xn yn / unn
一般地，用归纳法可以证明
ai ci 0 （即0 i 1） (i 1, 2, , n 1)
因此我们从关系式（2）解出待定系数为
1 b1
i ai
(i 2, 3, , n)
ii1bici1a/i
i1 i 1
(i 2, 3,
（3）
, n)
10 实现A=LU分解,按递推公式(4)计算
1, 2 , , n1
20 求解方程组 LY=f ,相应的递推公式是

y1 yi

f1 / b1 ( fi ai
yi1)
/(bi

ai i1)
(i 2,3,
（5）
, n)
30 求解方程组 UX=Y ,相应的递推公式是
xn yn

xi

yi

i xi1
（6）
(i n 1, n 2, , 2,1)
计算 1 2 n1 及 y1 y2 yn 的过程 称为追的过程，计算方程组的解 xn xn1 x1 的过程称为赶的过程，因此上述方法称为追赶法。
例4：用追赶法解方程组
lnk
1步 2步 … k-1步 k步
u1 j u2 j uk 1, j ukj
…
u1n
u2n


uk 1,n ukn Nhomakorabea




unn
n步
在解方程组时,对于右端项b也可不必经 过中间过程而按紧凑格式的方法直接得出y， 因为Ly=b ，所以
y1 b1
k 1
yk bk lkq yq (k 2,3, , n) q1