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

基于Matlab的n阶方阵的LDU分解实现1.引言矩阵的LDU分解是“矩阵理论与方法”课程中非常重要的一部分。LDU分解在实际工程应用中也非常广泛。LDU分解可以将一个矩阵分解为一个下三角矩阵和一个对角矩阵和一个上三角矩阵的乘积。LDU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。将系数矩阵A转变成等价两个矩阵L和U和对角矩阵的乘积，其中L和U分别是下三角和上三角矩阵，D为对角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LDU。即：Matlab是很好的处理矩阵的工具。它的功能非常强大，包括创建矩阵，对矩阵求逆，转置等操作非常简单，使其成为图像处理，信号分析等领域常用的工具。Matlab官方已经包括了对非奇异矩阵的LU分解函数[L,U]=lu(A)，为了加深对矩阵分解的理解，本文不采用Matlab官方的LU分解函数对矩阵A进行LDU分解，而是根据理论推导和编程实现LDU分解。2.程序设计2.1.输入合法检验LU分解需要被分解矩阵A满足如下条件：1)矩阵A为方阵2)A的顺序主子式故LU分解需先检验A为n阶方阵，然后检验A的n-1个顺序主子式Dk全不为0，才可进行LU分解。而检验主子式Dk可以在n-1次循环LU分解中进行，故先检验矩阵是否为方阵。代码如下2.2.n-1次循环LDU分解LDU分解本质上是高斯消元法。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。LDU分解主要分为两步：1根据高斯消元法对Ai()消元，消元矩阵为Li+1-1；2计算Li+1-1Ai()=Ai+1()以产生下一步迭代的Ai()。2.2.1.根据构造Lj,Lj-1(j=i+1)高斯消元Ai()，使Ai()第i+1列从第i+2行至n行都为0。构造消元矩阵Lj,Lj-1。首先判断是否为0，为0则无法继续分解，退出；否则继续。代码如下2.2.2.计算Li+1-1Ai()=Ai+1()Dk计算Li+1-1Ai()=Ai+1()并储存覆盖Ai()，并计算下次循环的主子式。代码如下2.3.计算并返回L和An-1()为结果通过n-1次循环累乘Li得到L矩阵，并且根据最后得到的An-1()矩阵分解出D和U矩阵。代码如下3.验证结果3.1.可LDU分解的矩阵验证构造可以LDU分解的矩阵ADk调用MyLU函数，分解A矩阵。计算L*D*U验证是否为原矩阵A（方法下同）3.2.输入矩阵A是否为空验证输入矩阵A中没有元素，输入不合法。3.3.输入矩阵A为方阵验证输入矩阵A维度不合法。3.4.顺序主子式出现0的错误验证因为计算011*LA时出现1220a，导致顺序主子式01211220aa。该情况为不合法输入。至此，函数功能和合法性检查全部验证完毕。4.心得体会通过“矩阵理论与方法”的理论指导和Matlab编程的实践经验，我基本掌握了矩阵分解中的LDU分解的推导过程和算法步骤。熟练掌握LDU分解，对今后研究LU分解、Doolittle分解、Crout分解、QR分解等矩阵分解方法的实现有非常大的帮助。对LDU分解的推导过程进行步骤分解和归纳，我将N维矩阵的LDU分解归纳总结为n-1次循环，每步循环进行n-i次元素除法（计算高斯消元系数），2次N维矩阵乘法（计算L和1iA）和1次元素乘法（计算顺序主子式的值）。按乘法计算，即时间复杂度为34121OnninOn。因为计算过程中需要1,jjLL辅助计算矩阵和Li+1-1Ai()=Ai+1()，即空间复杂度为2On。该LDU分解设计还不够快速，占用空间相对较多，是以后改进的方向。5.附录：程序源码function[L,D,U]=MyLU(A)%checkvalidityif(isempty(A))%checkAifisemptyerror('Aisempty!');end[N,D]=size(A);if(N~=D)%checkAifissquareerror('Aisnotasquare!');end%LDUdecomposingL=eye(N);det_k=A(1,1);forstep_n=1:N-1if(det_k==0)%checkifSequenceprincipalminorappears0error('Sequenceprincipalminoris0');endLi=eye(N);Li_inv=eye(N);forstep_row=step_n+1:Nmod=A(step_row,step_n)/A(step_n,step_n);Li(step_row,step_n)=mod;Li_inv(step_row,step_n)=-mod;endA=Li_inv*A;det_k=det_k*A(step_n+1,step_n+1);L=L*Li;endL=L;D=eye(N);U=zeros(N);fori=1:Nif(A(i,i)~=0)D(i,i)=A(i,i);U(i,:)=A(i,:)/D(i,i);elseD(i,i)=0;U(i,i)=1;endendend