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

1、 3 LU分解法 Gauss消去法的變形 知識預(yù)備： 1矩陣的初等行變換、初等矩陣及其逆、乘積 2矩陣的乘法 3上三角矩陣的乘積、單位下三角矩陣的乘積 4單位下三角矩陣的逆、可逆的上三角矩陣的逆 、Gauss消去法的矩陣解釋 Gauss消去法實質(zhì)上是將矩陣 A分解為兩個三角矩陣相乘。 我們知道，矩陣的初等行變換實質(zhì)就是左乘初等矩陣。 第一輪消元：相當于對A（1）左乘矩陣Li,即 L1Aa(2) 其中 1 小 J1) I211 an ai2 a1 n (2) (2) (1) Li I3101 ,a(2) a22 a2n J lan 1 i1 an c(2) J2) I n10 0 1 an2 a

2、nn 第二輪消元： 對應(yīng)于 L2A a(3) 一般地 LkA (k)A(k 1) k 1,2, ,n 1 ( 1) 10 其中 Lk 1 k 1k hk 鵜,i k 1,k 2, ,n akk 1 nk 整個消元過程為 Un U12 U1n Ln 1 Ln 2 L2L1 a A(n)記 U U22 U2n Unn 從而 A ( Ln 1 L n 2L2 L1 ) 1U L11L Ln12 Ln11U 其中L是單位下三角矩陣，即 1 121 L I31 1 321 ,hj a(j) ij R , ajj 2,3, ,n 1, ,n 1,3) 1 n1 1n2 LU的過程 【注】消元過程等價于A分

3、解成 回代過程是解上三角方程組的過程。 二、矩陣的三角分解 1、若將A分解成L? U,即A=L? U,其中L為單位下三角矩陣,U為 非奇異上三角矩陣,則稱之為對A的Doolittle 分解。 當A的順序主子式都不為零時，消元運算可進行，從而A存在唯一的 Doolittle 分解。 證明：若有兩種分解， A=L 1U1, A=L 2U2,則必有L1=L 2, U1=U 2。 因為LiUi=L2U2，而且Li, L2都是單位下三角矩陣，Ui,U2都是可逆 上三角矩陣 ,所以有 11 L21L1 U2U11 因此L21L1 U2U11 I （單位矩陣） 即 L1=L2， U1=U2、 2、若L是非奇

4、異下三角矩陣，U是單位上三角矩陣時，A存在唯一的 三角分解，A=LU,稱其為A的Crout分解（對應(yīng)于用列變換實施消元） 三、直接分解（ LU 分解）算法 LU 分解算法公式 按矩陣乘法 1 1 u11 u12 u1n A l 21 1 u22 u2n l31 l32 1 1 1 unn ln1 ln2 第一 步： 步： 利用 A中第 行、第一列元素確定 U的第一行、L的第一列 元素。由 a1j (1,0,0, ,0) (u1j,u2j , uij ,0, ,0)Tu1j (j 1,2, ,n) ai1 (li1 ,li2 , l ii 1,1, ,0) (u11,0,0)li1 u11 (i

5、 2,3, ,n) 得 u 1j =a1j ( j 1,2, ,n) l i1 =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)T IrkUkjUrj(j r, r 1,n) k 1 得U的第r行元素為 air (1 i1 ,1 i2, Urj arj 1 ii 1 ,1,0, (i r 1, r 2,n) lir (air 直接分解的 1 1 rk U kj , 1 ,0)

6、1 1 ik Ukr ) / U rr 1 緊湊格式： 11 21 l 31 n1 j r, r 1, ,n (U1r,U2r, (r 2,3, Urr ,0, ,o)T 1 1 ik U kr 1 ir U rr 1 ,n 1,i 1, ,n) U12 1n 2n l 32 l n2 y f u、 nn n U22 方程組的三角分解算法（LU分解） 對于方程組 Ax=b,設(shè)A=LU （Doolittle 分解）。 由于 Ly Ax b Ux b y 1、求解 Ly=b: y1b1, yi i 1 bilik yk ,(i k 1 2,3,n) 2、求解 Ux=y： Xny n / u nn

7、, Xi (yi UikXk)/Uii ,(in 1,n 2, ,2,1) (6) LU分解算法 步1，輸入A, b; 步 2，對 j=1,2,- ,n 求 u1j : u1 ja1 j , 對i=2,3， -,n 求 l i1 : li1 ai1 / u11 步3,對r=2,3， ,n 做(3.1 ) -(3.2): (3.1 ) urjarj r 1 lrkUkj,( j r,r 1, k k 1 r 1 i 1 ,n), (3.2 ) lir (a 1, ,n； r n)； y1 2,3,n，求yi : yi bi lik yk； Xn yn,對 i n 1,1 求xi : xi (yi

8、 n UikXk)/Uii i 1 輸出 xi(i 1,2,n);結(jié)束。 例子與程序: 【例】用LU分解求解方程組 22 3x2 47 7x2 2 4 5X3 解：對系數(shù)矩陣A進行LU分解 U11 2,山2 因此 U2j l32 a2j (a 32 2, U13 l 21u1 j ,所以 u223, U23 2, U222, U33 3, I 212, l 31 l 31U12 ) / U22 (a33131U13 l 32U23 )6 A 1 2 1 223 13 21 1 6 先解 Ly b,則 y 3, y2 1 2y15, y37 y1 2y26 。 再解 Ux y,解出X3 1,x2

9、 ( 5 x3)/3 2,x1(3 2x2 3x3)/ 22 程序： 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

10、; 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 Trian gular Part l=tril(a); for i=1: n l(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分解

11、法解方程組時，從第 r(r=1,2, ,n）步分解計算公式可 r 1 知Urj arjlrkUkj( j k 1 r,r 1, ,n) r 1 l ir(airlik ukr ) / urr (i r 1, ,n) k 1 當urr很小時，可能引起舍入誤差的累積、擴大。因此，可采用與列 主元消去法類似方法，將直接三角分解法修改為列主元三角分解法（與列 主元消去法在理論上是等價的），它通過交換A的行實現(xiàn)三角分解 PA=LU其中P為置換陣。 un 設(shè)第r-1步分解計算己完成，則有 um 21 l n1 第r步計算時為了避免用絕對值很小的數(shù)作除數(shù)，弓I進中間量： r 1 SiairlikUkr , (i r, , n) k 1 則有：urrSr, lir Sj Sr(ir 1,n) (1) 選主元：確定i r,使Sirmax Si r i n (2) 交換兩行：當irr