数值分析LU分解法
§3 LU 分解法
——Gauss 消去法的变形
知识预备：
1矩阵的初等行变换、初等矩阵及其逆、乘积
2矩阵的乘法
3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆
一、Gauss 消去法的矩阵解释
Gauss 消去法实质上是将矩阵A 分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。 第一轮消元：相当于对A
（1）
左乘矩阵L 1，即
)2()1(1A A L =
其中
)
1(11
)
1(1
1)2()2(2
)2(2)2(22)
1(1)
1(12)1(11)
2(131211,,1001
011a a l a a a a a a a A
l l l L i i nn n n n n =???
???
???????
?=?????
??
?????????---=
第二轮消元：对应于
)3()2(2A A L =
一般地
1,,2,1)
1()(-==+n k A A L k k k (1)
其中
n
k k i a a l l l L k kk
k ik
ik nk k
k k ,,2,1,,10111)()
(1 ++==???
????
??
?
??
???
????
?--=+整个消元过程为
U A A L L L L n n n 记)(1221=-- ?????
?
?
??
??
?
?
?=nn n
n
u u u u u u 22211211………（2） 从而
U L U L L L L U L L L L A n n n n ?===---------1
112121111221)(
其中L 是单位下三角矩阵，即
,1,,1,,3,2,,11
11)()(21323121???
?
??-===?????
??
?????????=n j n i l l l l l l L j jj j ij
a a ij n n …（3） 【注】消元过程等价于A 分解成LU 的过程
回代过程是解上三角方程组的过程。
二、矩阵的三角分解
1、若将A 分解成L?U ，即A=L?U ，其中L 为单位下三角矩阵,U 为非奇异上三角矩阵,则称之为对A 的Doolittle 分解。
当A 的顺序主子式都不为零时，消元运算可进行，从而A 存在唯一的Doolittle 分解。
证明：若有两种分解，A=L 1U 1，A=L 2U 2，则必有L 1=L 2，U 1=U 2。
因为L 1U 1=L 2U 2，而且L 1，L 2都是单位下三角矩阵，U 1,U 2都是可逆上三角矩阵,所以有
112112--=U U L L
因此 )(1
1211
2单位矩阵I U U L L ==--
即
L 1=L 2，U 1=U 2、
2、若L 是非奇异下三角矩阵，U 是单位上三角矩阵时，A 存在唯一的
三角分解，A=LU ，称其为A 的Crout 分解（对应于用列变换实施消元）
三、直接分解（LU 分解）算法
LU 分解算法公式——按矩阵乘法
??????
?
???
???????????
?????????=nn n
n
n n u u u u u u l l l l l A 222112112132312111
11 第一步：利用A 中第一行、第一列元素确定U 的第一行、L 的第一列
元素。由
),,2,1()0,,0,,,()0,,0,0,1(1211n j u u u u a j
T ij j j j ==?=
),,3,2()0,,0,()0,,1,,,(11
1111211n i u l u l l l a i T ii i i i =?=?=-
得 u 1j =a 1j n j ,,2,1( =) l i1=a i1/u 11n i ,,3,2( =)
第r 步:利用A 中第r 行、第r 列剩下的元素确定U 的第r 行、L 的第r 列元素（r=2，3，…，n ）.由
)
,,1,()0,,0,,,()0,,0,1,,,(1
1
21121n r r j u u l u u u l l l a rj
kj r k rk T
jj j j rr r r rj +=+=?=∑-=-得U 的第r 行元素为
∑-=+=-=1
1
,,1,,
r k kj rk rj rj n r r j u l a u
由
)
,,2,1()0,,0,,,()0,,0,1,,,(1
1
21121n r r i u l u l u u u l l l a rr
ir kr r k ik T
rr r r ii i i ir ++=+=?=∑-=-得
),,1,1,,3,2(/)(1
1
n r i n r u u l a l rr
kr r k ik ir ir +=-=-=∑-=…………
(4)
直接分解的紧凑格式：
方程组的三角分解算法（LU 分解）
对于方程组Ax=b ，设A=LU (Doolittle 分解)。
由于 ???==?=y
Ux b
Ly b Ax
1、求解Ly=b ：
∑-==-==1
1
11),,3,2(,,i k k ik i i n i y l b y b y (5)
2、求解Ux=y ：
)1,2,,2,1(,/)(,/1
--=-
==∑+=n n i u x u
y x u y x ii n
i k k ik
i i nn n n (6)
LU 分解算法
步1，输入A ，b ；
步2，对j=1,2,…,n 求,:111j j j a u u = 对i=2,3,…,n 求;/:11111u a l l i i i = 步3,对r=2,3,…,n 做（3.1）-(3.2): （3.1）∑-=+=-
=1
1),,,1,(,r k kj rk
rj rj n r r j u l
a u
（3.2）);;,,1(/)(1
1
n r n r i u u l
a l rr
r k kr ik
ir ir ≠+=-
=∑-=
步4，∑-=-
===1
1
11;:,,,3,2,i k k ik
i i i y l
b y y n i b y 求对
步5，ii n
i k k ik
i i i n n u x u
y x x n i y x /)(:1,,1,1
∑+=-=-==求对
步6，输出);,,2,1(n i x i =结束。
????
??????-=????????????????????-713542774322322x x x 解：对系数矩阵A 进行LU 分解
6
)(,2,2/)(1
,3,1,2,3,2,22332133133332222123132322322121223121131211=--===-===-=-=====u l u l a u u u u l a l u u u l a u l l u u u j j j 所以因此
????
?
???????
???
?????-=613322121121A 先解627,521,3,213121=-+-=-=-===y y y y y y b Ly 则。 再解
22/)323(,23/)5(,1,321323=--=-=--===x x x x x x y Ux 解出
程序：LU_factorization
%Not Select Column LU_factorization clear all
n=3;a=[2 2 3;4 7 7;-2 4 5];b=[3;1;-7]; %n=3;a=[1 4 7;2 5 8;3 6 11];b=[1;1;1]; %LU_factorazation for i=2:n
a(i,1)=a(i,1)/a(1,1); end a
for r=2:n for j=r:n s=0.;
for k=1:r-1
s=s+a(r,k)*a(k,j); end
a(r,j)=a(r,j)-s; end
for i=r+1:n s=0.;
for k=1:r-1
s=s+a(i,k)*a(k,r); end
a(i,r)=(a(i,r)-s)/a(r,r); end a end
%Extract Lower/Upper Triangular Part l=tril(a); for i=1:n l(i,i)=1; end
u=triu(a); l u
%Linear Lower Triangular Equation Solution y=l\b
%Linear Upper Triangular Equation Solution x=u\y
四、列主元LU 分解
当用LU 分解法解方程组时，从第r(r=1,2,…,n)步分解计算公式可知 ∑-=+=-
=1
1),,1,(r k kj
rk
rj rj n r r j u l
a u
),,1(/)(1
1
n r i u u l a l rr
kr r k ik ir ir +=-=∑-=
当rr u 很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列主元消去法在理论上是等价的），它通过交换A 的行实现三角分解PA=LU ，其中P 为置换阵。
设第r-1步分解计算己完成，则有
??
??
?
?
??????
????????
??????→--
-r n n r r r r n l l l u u u l
u u A
,1,,12221111
第r 步计算时为了避免用绝对值很小的数作除数，引进中间量：
∑-==-=1
1
),,(,r k kr
ik ir i n r i u l a S
则有：),,1(,n r i S S l S u r i ir r rr +===
（1） 选主元：确定i n
i r ir r S S i ≤≤=m ax ,使
（2） 交换两行：,时当r i r ≠交换A 的第r 行与第r i 行元素（相当于先
交换原始矩阵A 第r 行与第r i 行元素后，再进行分解计算得到的结果，且1≤ir l ）
（3） 进行分解计算 附：列主元LU 分解
a=[2 2 3;4 7 7;-2 4 5];b=[3;1;-7]; [l,u]=lu(a) y=l\b x=u\y
%x=inv(u)*inv(v)*b
[l,u,p]=lu(a) y=l\(p*b) x=u\y