数值分析5LU分解法

1、§3LU分解法Gauss消去法的变形知识预备：1矩阵的初等行变换、初等矩阵及其逆、乘积2 矩阵的乘法3 上三角矩阵的乘积、单位下三角矩阵的乘积4 单位下三角矩阵的逆、可逆的上三角矩阵的逆、Gauss消去法的矩阵解释Gauss消去法实质上是将矩阵A分解为两个三角矩阵相乘。我们知道，矩阵的初等行变换实质就是左乘初等矩阵。第一轮消元：相当于对A左乘矩阵Li,即LiA=A(2)其中a：a(1)a12.a1a1n1a(2)八(2).a22'aa2naai1*一，“，li1一a11q(2).an2'小2)ann一1l211L1=-13101+一*L-ln1001J第二轮消元：对应

2、于一般地k=1,2,，n-1(1)LkA(k)=A(k1)其中一1a(k)a_lik力=k1,k2,nakk1nk整个消元过程为U11U12Ln4Ln/L2L1A=A记U=U22U1nU2nUnn从而A=(Ln/Ln/L2L1)“U=L/lJ；/；式=LU其中L是单位下三角矩阵，即1l211311321lij(j)aij7?r),ajj彳=2,3,，nd=1,，n1(3)Jn11n2【注】消元过程等价于A分解成LU的过程回代过程是解上三角方程组的过程o二、矩阵的三角分解1、若将A分解成LR,即A=L?UI,其中L为单位下三角矩阵,U为非奇异上三角矩阵，则称之为对A的Doolittle分解。当A

3、的顺序主子式都不为零时，消元运算可进行，从而A存在唯一的Doolittle分解。证明：若有两种分解，A=L1Ui,A=L,2U2,则必有Li=L2,Ui=U>°因为LiUIi=L2UI2,而且Li,L2都是单位下三角矩阵，Ull,U2都是可逆上三角矩阵，所以有41L2Li=U2Ui因此L21L1=U2U/=I（单位矩阵）即Li=L2,Ui=U2-2、若L是非奇异下三角矩阵，U是单位上三角矩阵时，A存在唯一的三角分解，A=LU,称其为A的Crout分解（对应于用列变换实施消元）、直接分解（LU分解）算法LU分解算法公式一一按矩阵乘法一1l21A=131u11u12u22第一步：利

4、用A中第一行、第一列元素确定u1nu2nunn-U的第一行、L的第一列元素。由a1j=(1,0,0,0)(u1j,u2j,"uij,0,0)T=Wj(j=1,2,，n)a-=(li1,li2,lii,1,0)(Uii,0,0)T=li1Uii(i=2,3,n)得u砌(j=1,2,n)li1=an/uii(i=2,3,n)第r步：利用A中第r行、第r列剩下的元素确定U的第r行、L的第r列元素(r=2,3,，n).由r-1arj=(lr1,lr2,lrT,1,0,0)(U1j,U2j,Ujj,0,0)TlrkUkjUrj(j=，1,，n)k1得U的第r行元素为r1Urj=arj-、lrk

5、Ukj,j=r,r1,，nk1rJair=(lii,li2,lii1°，°)(Uir,U2r,Urr,0,，0)T=1，likUkr日k(i=r1,r2,n)得r1lir=(ai-likUkr)/Urr(r=2,3,n-1,i=r1,n)k1(4)直接分解的紧凑格式：U1n1U11U12U13l21l31sln1U22U23Jl32ln2UU2n2nnn方程组的三角分解算法(LU分解)对于方程组Ax=b,设A=LU(Doolittle分解)。Ly=b由于Ax=bu)、Ux=y1、求解Ly=b:i1y1=bi,yi=bi-Z鼠丫口。=2,3,，n)(5)kW2、求解Ux=y：

6、nxn=yn/Unn,xi=&-Z5卜*”5,(i=n-1,n-2,，2,1)(6)k=t1LU分解算法步1,输入A,b;步2,对j=1,2，,n求uj:u=a,对i=2,3，,n求l-：l.=2口/u11；步3,对r=2,3,n做(3.1)-(3.2):r-1(3)Urj=arj一£lrkUkj,(j=r,r+1,，n),k4rJ(3.2)iir=(air八likukr)/urr(i=r,1,，n；r=n)；k4i1步4,yi=bi,对i=2,3,n,求yi：yi=Dglikyk；k1n步5，xn=yn,对1=n-1,，1求Xi:Xi=(yi-ZUikXk)/Uiik=t1

7、步6,输出Xj(i=1,2,，n);结束。例子与程序：【例】用LU分解求解方程组|一223%一31477x2=1：-245X3一二7一解：对系数矩阵A进彳TLU分解u11=2,u12=2,u13=3,l21=2,l31=一1u2j=a2jl21u1j,所以u22=3,u23=1l32=(a32-l31u12)/u22=2,u22=2,u33=(a33-l31u13-l32u23)=6因此一11一223A=2131126一先解Ly=b,则y1=3,y2=12y1=5,y3=-7+y1-2y2=6。再解Ux=y,解出x3=1,x2=(5x3)/3=2,x1=(32x2-3x3)/2=2程序：LU_

8、factorization%NotSelectColumnLU_factorizationclearalln=3;a=223;477;-245;b=3;1;-7;%n=3;a=147;258;3611;b=1;1;1;%LU_factorazationfori=2:na(i,1)=a(i,1)/a(1,1);endaforr=2:nforj=r:ns=0.;fork=1:r-1s=s+a(r,k)*a(k,j);enda(r,j)=a(r,j)-s;endfori=r+1:ns=0.;fork=1:r-1s=s+a(i,k)*a(k,r);enda(i,r)=(a(i,r)-s)/a(r,r)

9、;endaend%ExtractLower/UpperTriangularPartl=tril(a);fori=1:nl(i,i)=1;endu=triu(a);lu%LinearLowerTriangularEquationSolutiony=lb%LinearUpperTriangularEquationSolutionx=uy四、列主元LU分解当用LU分解法解方程组时，从第r(r=1,2,n)步分解计算公式可r1知urj二arj-、l*ukj(j=r,r1,n)k1r1lir-(air，likukr)/urr(i=r1,n)k1当urr很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列主元消去法在理论上是等价的），它通过交换A的行实现三角分解PA=LU其中P为置换阵。u1n设第r-1步分解计算己完成，则有Uil21u22Jn1第r步计算时为了避免用绝对值很小的数作除数，引进中间量:,(i=r,n)r4Si-air二.likukrk1则有：u,r=Sr,%；S,(i=r1,n)(1) 选主元：确定ir,使Sir=max§r<<(2) 交换两行：当ir#时，交换A的第r行与第i