§3 LU 分解法——Gauss 消去法的变形知识预备：1矩阵的初等行变换、初等矩阵及其逆、乘积2矩阵的乘法3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆一、Gauss 消去法的矩阵解释Gauss 消去法实质上是将矩阵A 分解为两个三角矩阵相乘。

我们知道，矩阵的初等行变换实质就是左乘初等矩阵。

第一轮消元：相当于对A（1）左乘矩阵L 1，即)2()1(1A A L =其中)1(11)1(11)2()2(2)2(2)2(22)1(1)1(12)1(11)2(131211,,1001011a a l a a a a a a a Al 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)其中nk k i a a l l l L k kkk ikik nk kk k ,,2,1,,10111)()(1 ++==⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=+整个消元过程为U A A L L L L n n n 记)(1221=-- ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=nn nnu u u u u u 22211211………（2） 从而U L U L L L L U L L L L A n n n n ⋅===---------1112121111221)(其中L 是单位下三角矩阵，即,1,,1,,3,2,,1111)()(21323121⎪⎪⎭⎫⎝⎛-===⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=n j n i l l l l l l L j jj j ija 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因此 )(112112单位矩阵I U U L L ==--即L 1=L 2，U 1=U 2、2、若L 是非奇异下三角矩阵，U 是单位上三角矩阵时，A 存在唯一的三角分解，A=LU ，称其为A 的Crout 分解（对应于用列变换实施消元）三、直接分解（LU 分解）算法LU 分解算法公式——按矩阵乘法⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=nn nnn n u u u u u u l l l l l A 22211211213231211111 第一步：利用A 中第一行、第一列元素确定U 的第一行、L 的第一列元素。

由),,2,1()0,,0,,,()0,,0,0,1(1211n j u u u u a jT ij j j j ==⋅=),,3,2()0,,0,()0,,1,,,(111111211n 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,,,(1121121n r r j u u l u u u l l l a rjkj r k rk Tjj j j rr r r rj +=+=⋅=∑-=-得U 的第r 行元素为∑-=+=-=11,,1,,r k kj rk rj rj n r r j u l a u由),,2,1()0,,0,,,()0,,0,1,,,(1121121n r r i u l u l u u u l l l a rrir kr r k ik Trr r r ii i i ir ++=+=⋅=∑-=-得),,1,1,,3,2(/)(11n r i n r u u l a l rrkr r k ik ir ir +=-=-=∑-=…………(4)直接分解的紧凑格式：方程组的三角分解算法（LU 分解）对于方程组Ax=b ，设A=LU (Doolittle 分解)。

由于 ⎩⎨⎧==⇔=yUx bLy b Ax1、求解Ly=b ：∑-==-==1111),,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 uy x u y x ii ni k k iki 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 做（）-: （）∑-=+=-=11),,,1,(,r k kj rkrj rj n r r j u la u（）);;,,1(/)(11n r n r i u u la l rrr k kr ikir ir ≠+=-=∑-=步4，∑-=-===1111;:,,,3,2,i k k iki i i y lb y y n i b y 求对步5，ii ni k k iki i i n n u x uy 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 alln=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:na(i,1)=a(i,1)/a(1,1); end afor r=2:n for j=r:n s=0.;for k=1:r-1s=s+a(r,k)*a(k,j); enda(r,j)=a(r,j)-s; endfor i=r+1:n s=0.;for k=1:r-1s=s+a(i,k)*a(k,r); enda(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; endu=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)步分解计算公式可知 ∑-=+=-=11),,1,(r k kjrkrj rj n r r j u la u),,1(/)(11n r i u u l a l rrkr 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 lu u A,1,,12221111第r 步计算时为了避免用绝对值很小的数作除数，引进中间量：∑-==-=11),,(,r k krik 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 ni 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。