数值分析5 LU分解法.docx
数值分析5 LU分解法§3 LU分解法 Gauss消去法的变形 知识预备： 1矩阵的初等行变换、初等矩阵及其逆、乘积 2矩阵的乘法 3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆 一、Gauss消去法的矩阵解释 Gauss消去法实质上是将矩阵A分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。 第一轮消元：相当于对A左乘矩阵L1，即 L(1)1A=A(2) 其中 éê1ùé(1a(1)êa)1112La(1)1nùLê-l211úú1=êê-l3101ú,A(2)=êa(2)22La(2)ú2nêúêLLOOúêLOLú,êúë-lên10L01úûëa(2)La(2)ún2nnúû第二轮消元：对应于 L2)(3)2A(=A 一般地 LkA(k)=A(k+1)k=1,2,L,n-1 其中 1 a(1)li1i1=a(1)11é1êOêê1Lk=ê-lk+1kêêLê-lnkêë整个消元过程为 ùúú(k)úaikú,lik=(k),i=k+1,k+2,L,n1akkúOúú01úûéu11u12Lu1nêu22Lu2nLn-1Ln-2LL2L1A=A(n)记U=êêOLêunnêë从而 ùúú úúúû-1-11-1A=(Ln-1Ln-2LL2L1)-1U=L1L2LL-Ln-2n-1U=L×U 其中L是单位下三角矩阵，即 é1ùêlú121êúú,lij=L=êl31l321êúLLLOêúêëln1ln2LL1úûæi=2,3,L,nö÷,ç÷, a(jjj)çj=1,L,n-1èø(j)aij消元过程等价于A分解成LU的过程 回代过程是解上三角方程组的过程。 二、矩阵的三角分解 1、若将A分解成LU，即A=LU，其中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 2 因此 L2L1=U2U1=I(单位矩阵) 即 L1=L2，U1=U2、 2、若L是非奇异下三角矩阵，U是单位上三角矩阵时，A存在唯一的三角分解，A=LU，称其为A的Crout分解 -1-1三、直接分解算法 LU分解算法公式按矩阵乘法 é1ùêlúéu11u12Lu1n1ê21úêu22Lu2núêA=êl31l321OLêúêLLLOêêúunnêëêúllLL1ën1n2ûùúú úúúû第一步：利用A中第一行、第一列元素确定U的第一行、L的第一列元素。由 a1j=(1,0,0,L,0)×(u1j,u2j,Luij,0,L,0)T=u1jai1=(li1,li2,Llii-1,1,L,0)×(u11,0,L,0)T=li1×u11得 u1j=a1j (j=1,2,L,n) li1=ai1/u11(i=2,3,L,n) (j=1,2,L,n) (i=2,3,L,n) 第r步:利用A中第r行、第r列剩下的元素确定U的第r行、L的第r列元素.由 arj=(lr1,lr2,Llrr-1,1,0,L,0)×(u1j,u2j,Lujj,0,L,0)=ålrkukj+urjTk=1r-1(j=r,r+1,L,n)得U的第r行元素为 urj=arj-ålrkukj,k=1r-1j=r,r+1,L,n 3 由air=(li1,li2,Llii-1,1,0,L,0)×(u1r,u2r,Lurr,0,L,0)=ålikukr+lirurrTk=1r-1(i=r+1,r+2,L,n)得 lir=(air-ålikukr)/urrk=1r-1(r=2,3,L,n-1,i=r+1,L,n)(4) 直接分解的紧凑格式： u11 u12 u13 u1n 1 l21 u22 u23 u2n 2 l31 l32 M M ln1 ln2 unn n 方程组的三角分解算法 对于方程组Ax=b，设A=LU (Doolittle分解)。 由于 Ax=bÛí1、求解Ly=b： ìLy=bîUx=yy1=b1,yi=bi-ålikyk,(i=2,3,L,n) k=1i-12、求解Ux=y： xn=yn/unn,xi=(yi-k=i+1åunikxk)/uii,(i=n-1,n-2,L,2,1) LU分解算法 步1，输入A，b； 步2，对j=1,2,n 求u1j:u1j=a1j, 4 对i=2,3,n 求li1:li1=ai1/u11; 步3,对r=2,3,n 做-(3.2): urj=arj-ålk=1r-1r-1rkukj,(j=r,r+1,L,n), lir=(air-ålk=1ikukr)/urr(i=r+1,L,n;r¹n); i-1步4，y1=b1,对i=2,3,L,n,求yi:yi=bi-ålk=1ikyk; n步5，xn=yn,对i=n-1,L,1求xi:xi=(yi-步6，输出xi(i=1,2,L,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。 5 再解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); 6 l u %Linear Lower Triangular Equation Solution y=lb %Linear Upper Triangular Equation Solution x=uy 四、列主元LU分解 当用LU分解法解方程组时，从第r(r=1,2,n)步分解计算公式可知 urj=arj-r-1ålk=1r-1rkukj(j=r,r+1,L,n) lir=(air-ålikukr)/urrk=1(i=r+1,L,n) 当urr很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法，它通过交换A的行实现三角分解PA=LU，其中P为置换阵。 设第r-1步分解计算己完成，则有 éu11êlê21u22êMêA®êêêêMêêêëln1Lr-1LOur-1,r-1lr,r-1Mln,r-1arrManru1nùMúúMúúú Lur-1,núúLarnúMúúLannúûL第r步计算时为了避免用绝对值很小的数作除数，引进中间量： Si=air-ålikukrk=1,(i=r,L,n) 则有：urr=Sr,lir=SiSr(i=r+1,L,n) 7 选主元：确定ir,使Sir=maxSi r£i£n 交换两行：当ir¹r时,交换A的第r行与第ir行元素 进行分解计算 附：列主元LU分解 a=2 2 3;4 7 7;-2 4 5;b=3;1;-7; l,u=lu(a) y=lb x=uy %x=inv(u)*inv(v)*b l,u,p=lu(a) y=l(p*b) x=uy 8
