數(shù)值分析LU分解法

1、§ 3 LU 分解法 Gauss 消去法的變形知識(shí)預(yù)備：1 矩陣的初等行變換、初等矩陣及其逆、乘積2 矩陣的乘法3 上三角矩陣的乘積、單位下三角矩陣的乘積4 單位下三角矩陣的逆、可逆的上三角矩陣的逆一、 Gauss 消去法的矩陣解釋Gauss 消去法實(shí)質(zhì)上是將矩陣A 分解為兩個(gè)三角矩陣相乘。我們知道，矩陣的初等行變換實(shí)質(zhì)就是 左乘初等矩陣。第一輪消元：相當(dāng)于對(duì)A（1）左乘矩陣L1，即L1 A(1)A( 2)其中1a11(1)a12(1)a1(1n)l 211a22(2)a2(n2 )ai(11)L1l 3101,A( 2),l i1a11(1)l n100 1an(22)ann(2

2、)第二輪消元： 對(duì)應(yīng)于L2 A(2)A(3)一般地Lk A(k )A(k1)k1,2, , n1 （ 1）1其中11(k )Lk,likaik(k ),ik 1, k 2, , nl k 1k1akkl nk0 1整個(gè) 消元過(guò)程 為u11u12u1nLn 1 Ln 2L2L1 AA( n ) 記 Uu22u2n （ 2）unn從而A ( Ln 1 L n 2L2 L1 ) 1U L11L21Ln1 2 Ln11U L U其中 L 是單位下三角矩陣，即1Ll 2111, lijaij( j ) ,i2,3, n, （ 3）l 31l 32( j )a jjj1, n1l n1l n21【注】 消

3、元過(guò)程等價(jià)于A 分解成 LU 的過(guò)程回代過(guò)程是解上三角方程組的過(guò)程。二、矩陣的三角分解1、若將 A 分解成 L? U，即 A=L? U，其中 L 為單位下三角矩陣 ,U 為非奇異上三角矩陣 , 則稱之為對(duì) A 的 Doolittle 分解。當(dāng) A 的順序主子式都不為零時(shí)，消元運(yùn)算可進(jìn)行，從而A 存在 唯一 的Doolittle分解。證明： 若有兩種分解，A=L 1 U 1， A=L 2 U 2，則必有L1=L 2， U 1=U 2 。2因?yàn)?L1U1 =L 2U 2 ，而且 L1， L2 都是單位下三角矩陣，U 1,U2 都是可逆上三角矩陣 ,所以有L21L1U 2U11因此L21L1 U 2

4、U 11I (單位矩陣 )即L1=L 2 ， U 1=U 2、2、若 L 是非奇異下三角矩陣， U是單位上三角矩陣時(shí)， A 存在唯一的三角分解， A=LU，稱其為 A 的 Crout 分解（對(duì)應(yīng)于用列變換實(shí)施消元）三、直接分解（ LU 分解）算法LU 分解算法公式按矩陣乘法1l 211u11u12u1nu22u2nAl 31l321l n1l n 21unn第一步： 利用 A 中第一行、第一列元素確定U 的第一行、 L 的第一列元素。由a1 j(1,0,0,0)(u1 j , u2 j , uij ,0,0)Tu1 j( j1,2, n)ai1(l i1 , l i 2 ,l ii1,1,0)

5、(u11,0,0) Tl i1u11(i2,3, n)得u1j=a( j1,2, n )1jli1 =ai1 /u 11 (i2,3, n)第 r 步: 利用 A 中第 r 行、第 r 列剩下的元素確定U的第 r 行、 L 的第 r 列元素（ r=2 ，3， ， n） . 由3r1arj(l r1 ,l r 2 ,l rr1 ,1,0,0) (u1 j ,u2 j ,ujj ,0,0)Tl rk ukjurj ( j r , r 1, , n)k 1Urr1urjarjl rk ukj ,jr , r1, nk1,0) Tr1air(l i1 ,l i 2 ,l ii 1 ,1,0,0)(u1

6、r , u2r ,urr ,0,l ik ukrlir urrk1(ir1, r2, n)r1l ir(airlikukr ) / urr( r 2,3, n1,ir1, n)k1(4)u11u12u13ul21u22u23ul 31l 321n12n2ln1l n2unn n方程組的三角分解算法（LU分解）Ax=bA=LU (Doolittle)AxLybbyUx1Ly=bi 1y1 b1 , yibil ik yk ,(i2,3, , n)5k 12Ux=y4nxn yn / unn , xi ( yiuik xk ) / uii ,(i n 1, n 2, ,2,1) （ 6）ki 1L

7、U 分解算法步 1，輸入 A， b；步 2，對(duì) j=1,2,n求 u1 j : u1 ja1 j ,對(duì) i=2,3,n求 l i1 : li 1ai1 / u11 ;步 3, 對(duì) r=2,3,n做（ 3.1 ） -(3.2):r1（ 3.1 ） urjarjl rk ukj, ( jr , r1,n),k 1r1（ 3.2 ） l ir(airl ik ukr) / urr(ir1, n; rn);k1i1步 4， y1b1 ,對(duì) i2,3, , n,求 yi : yibil ik yk ;k1步 5， xnyn , 對(duì) in 1, ,1 求 xi : xi( yinuik xk ) / ui

8、ik i1步 6，輸出 xi (i1,2, n); 結(jié)束。例子與程序：【例】用LU 分解求解方程組223x23477x21245x37解： 對(duì)系數(shù)矩陣A 進(jìn)行 LU 分解u112,u122,u13 3,l 212, l 311u2 ja2 jl 21u1 j ,所以 u223,u23 1l32(a32l 31u12 ) / u22 2,u222,u33( a33 l 31u13 l 32u23 ) 6因此51223A21311216先解 Lyb, 則 y13, y21 2 y15, y37 y1 2 y2 6 。再解Uxy,解出 x31, x2( 5x3 ) / 32, x1(32x23x3

9、) / 22程序： LU_factorization%Not Select Column LU_factorizationclear 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_factorazationfor 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)

10、*a(k,r);enda(i,r)=(a(i,r)-s)/a(r,r);endaend6%Extract Lower/Upper Triangular Partl=tril(a);for i=1:nl(i,i)=1;endu=triu(a);lu%Linear Lower Triangular Equation Solutiony=lb%Linear Upper Triangular Equation Solutionx=uy四、列主元 LU 分解當(dāng)用 LU 分解法解方程組時(shí)，從第r(r=1,2,n) 步分解計(jì)算公式可r 1知urjarjl rk ukj( jr , r1, n)k 1r1l

11、ir(airlik ukr ) / urr(ir 1, n)k 1當(dāng) urr 很小時(shí)，可能引起舍入誤差的累積、擴(kuò)大。因此，可采用與列主元消去法類似方法，將直接三角分解法修改為列主元三角分解法（與列主元消去法在理論上是等價(jià)的），它通過(guò)交換 A的行實(shí)現(xiàn)三角分解PA=LU，其中 P 為置換陣。設(shè)第 r-1步分解計(jì)算己完成，則有u11u1nl 21u22Aur1, r 1ur 1,nl r ,r1arrarnl n1l n ,r 1 anrann7第 r 步計(jì)算時(shí)為了避免用絕對(duì)值很小的數(shù)作除數(shù)，引進(jìn)中間量：r1Siairl ik ukr, (ir , n)k 1則有： urrSr , lirSi Sr (ir1, , n)（ 1）選主元 ：確定 i r , 使 Sirmax Siri n（ 2）交換兩行： 當(dāng) i rr時(shí) , 交換 A 的