数值分析分解法

1、§ 3 LU分解法Gauss消去法的变形知识预备：1矩阵的初等行变换、初等矩阵及其逆、乘积2矩阵的乘法3上三角矩阵的乘积、单位下三角矩阵的乘积4单位下三角矩阵的逆、可逆的上三角矩阵的逆一、Gauss消去法的矩阵解释Gauss消去法实质上是将矩阵 A分解为两个三角矩阵相乘。我们知道，矩阵的初等行变换实质就是左乘初等矩阵。第一轮消元：相当于对A(1)左乘矩阵Li,即L* A(2)其中1LiI211I3101I n10,a(2)0 1aja(2)(2) a22a(2)dn2aa1 n(2) a2nJa(2)nnla5(1)an第二轮消元：对应于L2Aa(3)一般地LkA(k)A(k 1)k

2、 1,2, ,n1( 1)其中Lk1 k 1khk鵜,i k 1,k 2,akk,n1 nk整个消元过程为UnU12U1nLn 1 Ln 2 L2L1 aA(n)记 UU22U2nUnn从而A ( Ln 1 L n 2L2 L1 )1UL11LLn12Ln11U其中L是单位下三角矩阵，即1121L I311 321,hja(j)ijR ,ajj2,3, ,n1, ,n1 n11n2LU的过程【注】消元过程等价于A分解成回代过程是解上三角方程组的过程。二、矩阵的三角分解1、若将A分解成L?U,即A=L?U其中L为单位下三角矩阵，U为非 奇异上三角矩阵，则称之为对A的Doolittle 分解。当A

3、的顺序主子式都不为零时，消元运算可进行，从而A存在唯一的Doolittle 分解。证明：若有两种分解， A=L 1U1，A=L 2U2，则必有L1=L 2，U1=U 2。因为LiUi=L2U2，而且Li, L2都是单位下三角矩阵，Ui,U2都是可逆上三角矩阵 ,所以有11L21L1 U2U11因此L21L1 U2U11 I （单位矩阵）即L1=L2， U1=U2、2、若L是非奇异下三角矩阵，U是单位上三角矩阵时，A存在唯一的三角分解，A=LU,称其为A的Crout分解（对应于用列变换实施消元）三、直接分解（ LU 分解）算法LU分解算法公式按矩阵乘法11u11 u12u1nAl 211u22u

4、2nl31l32111unnln1ln2第一步：步：利用A中第行、第一列元素确定U的第一行、L的第一列元素。由a1j(1,0,0,0) (u1j,u2j, uij ,0, ,0)Tu1j(j1,2, ,n)ai1(li1,li2 ,l ii 1,1, ,0)(u11,0,0)li1 u11(i2,3, ,n)得u1j =a1j( j 1,2,n)li1 =ai1/u 11(i 2,3,n)第r步：利用A中第r行、第r列剩下的元素确定U的第r行、L的第r列元素（r=2 , 3,，n）.由a rj ( 1 r1 , 1 r2 ,1 rr 1 ,1,0,0) (Ulj ,U2j,Ujj ,o.,o)

5、TIrkUkjUrj(j r, r 1,n)k 1得U的第r行元素为Urj arj11 rk u kj ,1j r, r1,nair(1 i1 ,1 i2,1 ii 1 ,1,0,0)(U1r,U2r,U rr ,0,o)T11 ik U kr 1 ir U rr1(ir 1, r2,n)lir(air11 ik Ukr ) / U rr1(r 2,3,n 1,i1,n)直接分解的紧凑格式：1121l 31n1U121n2nl 32l n21fu、nnnU22方程组的三角分解算法（LU分解）对于方程组 Ax=b，设A=LU （Doolittle 分解）。由于1、求解 Ly=b:Ax bLy b

6、Ux yy1b1, yi bii 11ikyk,(i2,3,n)2、求解 Ux=y：Xnyn/Unn,Xj 卜UikXk)/Uii ,(in 1,n2,2,1) (6)k i 1LU分解算法步1，输入A, b;ai j，步 2，对 j=1,2,n 求 u1j : u1 j对 i=2,3,n 求 li1 : li1aM /u11；步 3,对 r=2,3,n 做()-:r 1()u rj a rj1 rk ukj ,( j r,r 3,n),k 1r 1()lir &1ikUkr)/Urr (i 1,；rn)；k 1i 1步 4, y1b1,对i 2,3, ,n，求mbihyk；k 1n步

7、5, Xnyn,对i n 1,1 求xXi(yiUikXk)/uHk i 1步6，输出xi(i 1,2,n);结束。例子与程序：【例】用LU分解求解方程组一223x2477x22 4 5 X3解：对系数矩阵A进行LU分解U112,U122,U133212,131U2j a2j l 21U1j ,所以 u223, U23132(a32l3152)/U222,U222,U33(a33131U13 l 32U23 )6因此A121223132116先解Ly b,则 y3, y21 2y15,y37 y1 2y26 。再解Uxy,解出X31,x2( 5 x3)/32,x1(3 2x2 3x3)/ 22

8、程序： 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);endafor r=2:nfor j=r:ns=0.;for k=1:r-1s=s+a(r,k)*a(k,j);enda(r,j)=a(r,j)-s;endfor i=r+1:ns=0.;for k=1:r-1s=s+a(i,k)*a

9、(k,r);enda(i,r)=(a(i,r)-s)/a(r,r);enda end%Extract Lower/Upper Trian gular Part l=tril(a);for i=1: nl(i,i)=1; end u=triu(a); l u %Lin ear Lower Trian gular Equati on Soluti on y=lb%Linear Upper Triangular Equation Solution x=uy四、列主元LU分解当用LU分解法解方程组时，从第r(r=1,2,n）步分解计算公式可r 1知Urj arjlrkUkj( jk 1r,r 1,n)

10、r 1l ir(airlik ukr ) / urr(i r 1,n)k 1当urr很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列 主元消去法在理论上是等价的），它通过交换A的行实现三角分解PA=LU其中P为置换阵。un设第r-1步分解计算己完成，则有um21l n1第r步计算时为了避免用绝对值很小的数作除数，弓I进中间量：r 1Si airlikUkr, (i r, , n)k 1则有：urrSr,lirSj Sr(i r 1,n)(1) 选主元：确定i r,使Sir max Sir i n(2) 交换两行：当ir r时，交换a的第r行与第ir