数值分析5LU分解法

1、§ 3 LU分解法Gauss消去法的变形知识预备：1矩阵的初等行变换、初等矩阵及其逆、乘积2矩阵的乘法3 上三角矩阵的乘积、单位下三角矩阵的乘积4 单位下三角矩阵的逆、可逆的上三角矩阵的逆一、Gauss消去法的矩阵解释Gauss消去法实质上是将矩阵 A分解为两个三角矩阵相乘。 我们知道，矩阵的初等行变换实质就是左乘初等矩阵。L1,即第一轮消元：相当于对A(1 左乘矩阵L1A二A其中-1121a(2)J1)a12a1n 1J2)a22畀a2na(1)1 一ai1+,打1 一a11a(2)an2j2)ann12第二轮消元：对应于L2A=A般地LkA(kA(k 1)k =1,2, ,n -

2、1其中,na(k)1 ik令,i 二 k 1,k 2,akk-1 nk整个消元过程为U11UmLnLnd L2LjA =A记 U -U22U2n从而UnnA =(LnLn L2L1) U其中L是单位下三角矩阵，即一11211 311 321J【注】消元过程等价于A分解成回代过程是解上三角方程组的过程。n11 n2LU的过程a( j) aij丽，jjU =2,3,，nJ = h，n -1二、矩阵的三角分解1、若将A分解成L?J,即A=L?U,其中L为单位下三角矩阵，U为非奇 异上三角矩阵，则称之为对 A的Doolittle 分解。当A的顺序主子式都不为零时，消元运算可进行，从而A存在唯一的Doo

3、little 分解。证明：若有两种分解，A=L1U， A=LsU2，则必有L1=L2，U=U。因为L1Ui=L2Lb，而且L1，L2都是单位下三角矩阵，U,U2都是可逆上三角矩阵，所以有4-1L；L1 75因此L21L1 二 U2UF = I (单位矩阵)即Li=I_2, Ui=L2-2、若L是非奇异下三角矩阵，U是单位上三角矩阵时，A存在唯一的 三角分解，A=LU,称其为A的Crout分解（对应于用列变换实施消元）、直接分解（LU分解）算法LU分解算法公式按矩阵乘法一1l21 1A =1 311 32U11 U12U22UmU2nUnn第一步：利用A中第一行、第一列元素确定U的第一行、l的第

4、一列元素。由a1 j = (1,0,0,0) (u1 j , u2 j , uij ,0,0)= u1j(j = 1,2,n)玄祁=(1 订，1 i2,1 ii1,0)(U11Q,0)= 1i1 un (i 二 2,3,n)得 u1j =a1j(j =1,2, n )1 i1 =a1 /u 11 (i =2,3, n)第r步：利用A中第r行、第r列剩下的元素确定U的第r行、L的第r列元素(r=2 , 3,，n).由r A.arj =(1r1,lr2, lrr,1,0,0) (U1j,U2j,Ujj,0,0)IrkUkj Urj (j ="1,，n)得U的第r行元素为j r,r -1,

5、 nr -4Urj arj 二 1 rk Ukj ,k=irair =(lii,li2,hi 打,0,0) (Uir,U2r, Urr ,0,，0)丁 八 likUkr lirUrrk_4(i = r 1, r 2, n)r丄lir =(a-'TikUkr)/Urr (r = 2,3厂，n _1,i = r 1,n)k-1直接分解的紧凑格式：U11U12U13 U in 121U22U23U 2nl 31n1l32nnn方程组的三角分解算法(LU分解) 对于方程组 Ax=b，设A=LU (Doolittle 分解)。禹工ly =b由于Ax = b =丿lx = y1、求解 Ly=b:i

6、 -4y1 二bi,yi-bi二hkyk,(i -2,3,n) (5)2、求解 Ux=y：nXn =yn /Unn,Xi 二卜-' UikXk)/® ,(i = n -1, n - 2,2,1) 一( 6) k=t+LU分解算法步1，输入A, b;步 2，对 j=1,2,n求 u1 j : U1j = a1j,对 i=2,3,n 求l i1 ： li = ai1 /u11;步 3,对 r=2,3,n 做( 3.1 ) -(3.2):r -1(3)Urj 二aj -'TrkUkj,(j 二1,n),k 1 r J(3.2 ) lir = (air - 7 likukr)

7、/ urr (i = r T,，n； r = n)；k 4iy1 = b1 ,对 i = 2,3, n,求 yi : yi - bi.1 l ik yk ;k4nXn =yn,对 i =n -1,1 求Xi : Xi =(yi -為 UikXk)/Uiik=t 1输出xi（i =12，n）;结束。例子与程序:_223 14一2【例】用LU分解求解方程组75 一"7 一解：对系数矩阵A进行LU分解Un =2,比2 =2,比3 = 321 = 2,打1 二-1二 3331U13 一32口23)= 6U2j 二 a2j21U1j,所以 U22 =3,U23 =1*2 =(a32 1 丨 3

8、1U12 ) / U22 = 2,U 22 = 2, U33因此A =2 13 1L-1 2 1一16 一2 231=1-2力二 _5, y3_1先解 Ly = b,则 yi = 3, y2一7 y1再解Ux = y,解出 x3 = 1, x2 = ( 一5 x3) /3 = -2, x = (3 2x2 3x3) /2 = 2 程序：LU_factorizati on%Not Select Colu mn 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;%

9、LU_factorazati onfor 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(k,r);enda(i,r)=(a(i,r)-s)/a(r,r);endaend%Extract Lower/Upper Trian gular Partl=tril(a);for i=1: nl(i,i)=1;endu=triu(a);u%Lin ear Lo

10、wer Trian gular Equati on Soluti on y=lb%Linear Upper Triangular Equation Solution x=uy四、列主元LU分解当用LU分解法解方程组时，从第r（r=1,2,n）步分解计算公式可r 二知 Urj =aj 一、 IrkUkj（j =r,r T,，n）kr Jl ir =（alikUkr）/Urr（i = r，n）k=1当U“很小时，可能引起舍入误差的累积、扩大。因此，可采用与列主元消去法类似方法，将直接三角分解法修改为列主元三角分解法（与列 主元消去法在理论上是等价的），它通过交换A的行实现三角分解PA=LU其中P为置换阵。设第r-1步分解计算己完成，则有Uml 21 u22J n1,（i 二r, ,n）第r步计算时为了避免用绝对值很小的数作除数，弓I进中间量:则有：urr = Sr ,hr = Si Sr （i = r 竄，n）r 4Si = airlik u krk 1(1) 选主元：确定ir,使Sir = maxSjr<<(2) 交换两行：当ir =r时,交换A