基于Matlab-的-n阶非奇异方阵的LU分解实现

基于Matlab的n阶方阵的LDU分解实现1.引言矩阵的LDU分解是“矩阵理论与方法〃课程中非常重要的一局部。LDU分解在实际工程应用中也非常广泛。LDU分解可以将一个矩阵分解为一个下三角矩阵和一个对角矩阵和一个上三角矩阵的乘积。LDU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。将系数矩阵A转变成等价两个矩阵L和U和对角矩阵的乘积，其中L和U分别是下三角和上三角矩阵，D为对角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A二LDU。即：包括创立矩阵，对A=LDU=包括创立矩阵，对Matlab是很好的处理矩阵的工具。它的功能非常强大，矩阵求逆，转置等操作非常简单，使其成为图像处理，信号分析等领域常用的工具。Matlab官方已经包括了对非奇异矩阵的LU分解函数[L,U]=lu(A)，为了加深对矩阵分解的理解，本文不采用Matlab官方的LU分解函数对矩阵A进行LDU分解，而是根据理论推导和编程实现LDU分解。2.程序设计2.1.输入合法检验LU分解需要被分解矩阵A满足如下条件：矩阵A为方阵A的顺序主子式故LU分解需先检验A为n阶方阵，然后检验A的n-1个顺序主子式全k不为0,才可进行LU分解。而检验主子式可以在n-1次循环LU分解中k进行，故先检验矩阵是否为方阵。代码如下%checkvalidityif(isempty(A))%checkAifisemptyerror(*Aisempty!*);end[NfD]=size(A);if(ND)%checkAifissquareerror('Aisnotasquare!1);end2・2・n-1次循环LDU分解LDU分解本质上是高斯消元法。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。LDU分解主要分为两步：1根据高斯消元法对Ai消元，消元矩阵为Li；2计算LiA/A/i以产i1 i1生下一步迭代的Ai。2.2.1.根据 构造L,Li(ji1)JJ高斯消元Ai，使Ai第i+1列从第i+2行至n行都为0。构造消元矩阵L,L1。首先判断是否为0，为0那么无法继续分解，退出；否那么继JJ k续。代码如下%LUdecomposingL=eye(N);det_k二Adjl);forstep^n=1:N-1if(det_k二二0)^scheckifSequenceprincipalminorappears0error!'Sequenceprincipalminoris『)；endLi=eye(N);Li_inv=eye(M);forstep_row=step_n+l:Nnod=A(step_row,step_n)/A(step_n,step_n);Lilstep.roWjStep^n)=mod；Li_inv(step_row,step_n)=-mod;end

2.2.2.计算2.2.2.计算Li11AiAi1计算LiAiAi1并储存覆盖Ai,并计算下次循环的主子式i1代码如下A=Li_inv*A;det_k=det_k*A(step_rn-1,$tep_n+l):2.3.计算并返回L和Ani为结果通过n-1次循环累乘L得到L矩阵，并且根据最后得到的Ani矩阵分解出iD和U矩阵。代码如下0)0)A(i,i);A(i,:)/D(ifi);0;1；=L;=eye(N);=zeros(N);ri=l:Nif(A(ifi)D(ifi)=U(i,:)=elseD(i,i)=U(i’i)二end-end-end3.验证结果3・1・可LDU分解的矩阵验证构造可以LDU分解的矩阵A»A=[52-40;21-2l;-4-250;0102]TOC\o"1-5"\h\z5 2-42 1-2-4 -2 50 10调用MyLU函数，分解A矩阵。»[LfDfU]=MyLU(A)L二1・00000000・40001.000000—0.8000-2・00001.0000005・00002・00001-0000D=5-000000000.200000001.00000000-7・0000U二1.00000.4000-0.8000001.0000-2,00005・0000001.00002.00000001,0000计算L*D*U验证是否为原矩阵A〔方法下同〕»L*D*Uans=5・00002・0000-4・000002.00001.0000-2-00001.0000-4・0000-2>00005・0000001.000002・00003・2・输入矩阵A是否为空验证»A=rand(3,3) (:」：)=[]»：LD丄]=MyLU(A)错误使用MMU(line4)AlsemptyJ输入矩阵A中没有元素，输入不合法。3・3・输入矩阵A为方阵验证>>A=[1,Z5,6]1 2 34 5 6>>[L,艮匸]=MvLU(A)错误使用M¥LU(Line9)Aisnotasquare1输入矩阵A维度不合法。3・4・顺序主子式出现0的错误验证

TOC\o"1-5"\h\zL 3 62 6 127 4 3»丄D「匸二MyLU(A)错误使用MvLU(lij^eL7：iScqu亡neepri.ncipslminori..?.0因为计算L,i*A（o）时出现a（i）=0，导致顺序主子式„=a（o）a（i）=0。该情况12221122为不合法输入。至此，函数功能和合法性检查全部验证完毕。4.心得体会通过“矩阵理论与方法〃的理论指导和Matlab编程的实践经验，我根本掌握了矩阵分解中的LDU分解的推导过程和算法步骤。熟练掌握LDU分解，对今后研究LU分解、Doolittle分解、Crout分解、QR分解等矩阵分解方法的实现有非常大的帮助。对LDU分解的推导过程进行步骤分解和归纳，我将N维矩阵的LDU分解归纳总结为n-1次循环，每步循环进行n-i次元素除法〔计算高斯消元系数〕2次N维矩阵乘法〔计算L和A（i+i）〕和1次元素乘法〔计算顺序主子式的值〕按乘法计算，即时间复杂度为°'-1）（n7…2n3…1））K03因为计算过程中需要Lj辅助计算丄-丄厶…A：矩阵和—AA1，即空间复杂度为°O。该LDU分解设计还不够快速，占用空间相对较多，是以后改进的方向。5.附录：程序源码function[L,D,U]=MyLU(A)%checkvalidityif(isempty(A))%checkAifisemptyerror('Aisempty!');end[N,D]=size(A);if(N~=D)%checkAifissquareprincipalminorappears0error('Sequenceprincipalminoris0');endLi=eye(N);Li_inv=eye(N);forstep_row=step_n+1:Nmod=A(step_row,step_n)/A(step_n,step_n);endendA=Li_inv*A;det_k=det_k*A(step_n+1,step_n+1);L=L*Li;error('Aisnotasquare!');end%LDUdecomposingL=eye(N);det_k=A(1,1);forstep_n=1:N-1Li(step_