数值分析5-LU分解法

§3 LU分解法

——Gauss消去法的变形

知识预备：

1矩阵的初等行变换、初等矩阵及其逆、乘积

2矩阵的乘法

3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆

一、Gauss消去法的矩阵解释

Gauss消去法实质上是将矩阵A分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。

第一轮消元：相当于对A（1）

左乘矩阵L1，即

L(1)1AA(2)

其中

1(1)a(1)a1112a(1)1nLl2111l3101,A(2)a(2)22a(2)2n,ln1001a(2)a(2)n2nn

第二轮消元：对应于

L(2)2AA(3)

一般地

LkA(k)A(k1)k1,2,,n1 ……………（1）

其中

推荐精选

a(1)li1i1a(1)11

11Lklk1klnk整个消元过程为

(k)aik,lik(k),ik1,k2,,n1akk01u11u12u1nu22u2nLn1Ln2L2L1AA(n)记Uunn从而

………（2） 1111A(Ln1Ln2L2L1)1UL1L2Ln2Ln1ULU

其中L是单位下三角矩阵，即

1l121,lijLl31l321ln1ln21i2,3,,n,,…（3） a(jjj)j1,,n1(j)aij【注】消元过程等价于A分解成LU的过程

回代过程是解上三角方程组的过程。

二、矩阵的三角分解

1、若将A分解成LּU，即A=LּU，其中L为单位下三角矩阵,U为非奇异上三角矩阵,则称之为对A的Doolittle分解。

当A的顺序主子式都不为零时，消元运算可进行，从而A存在唯一的Doolittle分解。

证明：若有两种分解，A=L1U1，A=L2U2，则必有L1=L2，U1=U2。

因为L1U1=L2U2，而且L1，L2都是单位下三角矩阵，U1,U2都是可逆上三角矩阵,所以有

11L2L1U2U1推荐精选

因此 L2L1U2U1I(单位矩阵) 即

L1=L2，U1=U2、

2、若L是非奇异下三角矩阵，U是单位上三角矩阵时，A存在唯一的三角分解，A=LU，称其为A的Crout分解（对应于用列变换实施消元）

11

三、直接分解（LU分解）算法

LU分解算法公式——按矩阵乘法

1lu11u12u1n121u22u2nAl31l321unnll1n1n2 第一步：利用A中第一行、第一列元素确定U的第一行、L的第一列

元素。由

a1j(1,0,0,,0)(u1j,u2j,uij,0,,0)Tu1jai1(li1,li2,lii1,1,,0)(u11,0,,0)Tli1u11得 u1j=a1j (j1,2,,n) li1=ai1/u11(i2,3,,n)

(j1,2,,n)

(i2,3,,n)

第r步:利用A中第r行、第r列剩下的元素确定U的第r行、L的第r列元素（r=2，3，…，n）.由

arj(lr1,lr2,lrr1,1,0,,0)(u1j,u2j,ujj,0,,0)lrkukjurjTk1r1(jr,r1,,n)得U的第r行元素为

urjarjlrkukj,k1r1jr,r1,,n推荐精选

由

air(li1,li2,lii1,1,0,,0)(u1r,u2r,urr,0,,0)likukrlirurrTk1r1(ir1,r2,,n)得

lir(airlikukr)/urrk1r1(r2,3,,n1,ir1,,n)…………

(4)

直接分解的紧凑格式：

u11 u12 u13 … … u1n 1 l21 u22 u23 … … u2n 2 l31 l32   ln1 ln2 unn n

方程组的三角分解算法（LU分解）

对于方程组Ax=b，设A=LU (Doolittle分解)。

由于 Axb1、求解Ly=b：

Lyb

Uxyy1b1,yibilikyk,(i2,3,,n) …………………（5）

k1i12、求解Ux=y：

xnyn/unn,xi(yiki1unikxk)/uii,(in1,n2,,2,1)…（6）

LU分解算法

步1，输入A，b；

步2，对j=1,2,…,n 求u1j:u1ja1j,推荐精选

对i=2,3,…,n 求li1:li1ai1/u11; 步3,对r=2,3,…,n 做（3.1）-(3.2): （3.1）urjarjlk1r1r1rkukj,(jr,r1,,n),

（3.2）lir(airlk1ikukr)/urr(ir1,,n;rn);

i1步4，y1b1,对i2,3,,n,求yi:yibilk1ikyk;

n步5，xnyn,对in1,,1求xi:xi(yi步6，输出xi(i1,2,,n);结束。

ki1uikxk)/uii

例子与程序： 【例】用LU分解求解方程组 223x23477x1 2245x37解：对系数矩阵A进行LU分解

u112,u122,u133,l212,l311u2ja2jl21u1j,所以u223,u231l32(a32l31u12)/u222,u222,u33(a33l31u13l32u23)6因此

1A2112122331 6先解Lyb,则y13,y212y15,y37y12y26。

推荐精选

再解

Uxy,解出x31,x2(5x3)/32,x1(32x23x3)/22

程序：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)步分解计算公式可知 urjarjr1lk1r1rkukj(jr,r1,,n)

lir(airlikukr)/urrk1(ir1,,n)

当urr很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列主元消去法在理论上是等价的），它通过交换A的行实现三角分解PA=LU，其中P为置换阵。

设第r-1步分解计算己完成，则有

u11l21u22Aln1r1ur1,r1lr,r1ln,r1arranru1n

ur1,narnann第r步计算时为了避免用绝对值很小的数作除数，引进中间量：

Siairlikukrk1,(ir,,n)推荐精选

则有：urrSr,lirSiSr(ir1,,n)

rin（1） 选主元：确定ir,使SirmaxSi

（2） 交换两行：当irr时,交换A的第r行与第ir行元素（相当于先

交换原始矩阵A第r行与第ir行元素后，再进行分解计算得到的结果，且lir1）

（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

（注：可编辑下载，若有不当之处，请指正，谢谢!）

推荐精选