数值分析LU分解法

1、 3 LU分解法 Gauss消去法的变形 知识预备： 1矩阵的初等行变换、初等矩阵及其逆、乘积 2矩阵的乘法 3上三角矩阵的乘积、单位下三角矩阵的乘积 4单位下三角矩阵的逆、可逆的上三角矩阵的逆 、Gauss消去法的矩阵解释 Gauss消去法实质上是将矩阵 A分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。 第一轮消元：相当于对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 第二轮消元： 对应于 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分解（对应于用列变换实施消元） 三、直接分解（ 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,设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);结束。 例子与程序: 【例】用LU分解求解方程组 22 3x2 47 7x2 2 4 5X3 解：对系数矩阵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的行实现三角分解 PA=LU其中P为置换阵。 un 设第r-1步分解计算己完成，则有 um 21 l n1 第r步计算时为了避免用绝对值很小的数作除数，弓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