矩阵的LU分解(自编MATLAB)实验报告（完整规范版）.docxVIP

矩阵的LU分解（自编MATLAB实 验报告 1矩阵的LU分解 1.1 LU分解原理 定理：设AC n n，如果A的顺序主子式 ?■???? ? >? >? 少?少? ? ?? ? ? ? ? ?? ?? ? ? ????? ?? ?? ?? ? ?■??? ? ? >? >? 少?少? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ????? ? T工0 ??? -???? ??? -???? ???- ????_?? 则存在唯一的主对角线上元素全为 1的下三角矩阵L与唯一的 上三角矩阵U，使得 A=LU . 证明：对矩阵A的阶数使用数学归纳法. 显然，当n=1时，？???=1 ??????就是唯一的分解式。现假定对n— 阶矩阵，定理的结论成立。对 A进行分块 A=(???-?? ???) A=(?雰？ ????)? 其中？??？??€???「??.由于n-1阶矩阵？??_??的k阶顺序主子式 就是A的k阶主子式（k=1,2,…，n-2）,故它们都不为零?从而由 归纳法假设，？??_??有唯一的LU分解 ???_?? = ???_??????_?? 其中???_??的主对角线上的元素都1?由于 1???-?? = 1????？？？？? ?? ? ? ?? ?? ? 2 ?■ ?■ z-????？？？？? ???? 1???-?? = 1 ??? ?？？？？ ? ? ? ? ? ? ? ? ? ? 2 ?■ ?■ z- ??? ?？？？？ ? ???? ????? | = | ???-?????-??工0 ???- ?????？？ 所以???-??及???_??是n-1阶可逆矩阵 先假设已有A=LU，其中 ???■??L=( ???■?? L=( ???■ ???-?? U= ( ?? ?? ????=( ????=( ■ ■ ■ ■ \ ???-?????-?? ?切?-?? ?? ??€?尹?是待定向量。作乘积 ???_???? ???-?? ??? ????+ ????? =( ?芻？ ？???!=A 则?? ?必须满足 I^n-lP = ” YTUn-l = ^2 > bnn + ytP = ann 注意到及〃-I都是阶可逆矩阵，则由上式可惟一确定 P = Yt =bnn = ann 一 yTp 这就证明了 /的分解的存在性和唯一性. 1.2 LU分解算法 当ZI阶矩阵满足定理的条件时，可以用初等变换的方法求出L 和U. 因为当A=LU时，由于L可逆，故必存在可逆矩阵P使得 PL = I 即PA=PLU=U.也就是说，可以先对A施行行的初等变换得出 上三角矩阵（/,而矩阵P可以通过对单位矩阵/进行相同的行初等变 换得出，即 P(A,I) =(PA,P1) =(U9P) 于是A = P~1U,为保持P为下三角矩阵（从而P"也是下三角矩 阵），在进行行初等变换时，不能进行行的对换，上行的倍数应加到 下行的对应元. 1-3 LU分解用于解方程组 矩阵的三角分解在求解线性方程组时十分方便?如对线性方程组 Ax = ?设4 = 我们先求解方程组。=b.由于E是下三角矩阵， 则解向量y可以通过依次求出其分量ylfy2f-yn而求出，在求解方程 组〃尤=y.解向量兀可以通过该方程组依次求出分量xn,--fx2fx1而快 速得出?于是由两个方程组Ux = y9 Ly = b的求解而给出 LUx = Ly = b = 4尤的解? 1.4程序流程图 1.5 MATLAB 程序 function f=LU_decom(A) [m, n]=size(A) if m~=n fprintf( 'Error:m and n must be equal!m=%d,n=%d\n' ,m,n) end for i=1: n-1 if (det(A(1:i,1:i))==0) fprintf( 'Error:det A(%d,%d)=0!\n' ,i,i) flag='failure' return ; else flag二'ok'; end end L=eye (n); U=zeros (n); for i=1: n U(1,i)=A(1,i); end for r=2: n L(r,1)=A(r,1)/U(1,1); end for i=2: n for j=i:n z=0; for r=1:i-1 z=z+L(i,r)*U(r,j); end U(i,j)=A(i,j)-z; end if abs(U(i,i))<eps flag='failure' return ; end for k=i+1: n m=0; for q=1:i-1 endm=m+L(k,q