matlab编写LDR分解程序

程序运行步骤：

1、将源程序复制到Matlab的Command Window窗口中，如下图所示：

2、回车，则出现如下图所示提示：

3、输入一个矩阵如下图所示：

4、回车，得出如下图所示结果：

A=input('请输入一个矩阵A\n'); %提示输入一个矩阵A b=size(A); %b(1)表示A的行数，b(2)表示A的列数n=b(1);

%判断矩阵阶次，是否方阵

if b(1)~=b(2)

error('不能分解,非方阵')

end

%判断矩阵是否可逆

if n~=rank(A)

error('不能分解，非满秩矩阵')

end

%判断矩阵是否有惟一LDR分解

for k=1:n-1

d=A(1:k,1:k);

D(k)=det(d);

if(D(k)==0)

error('某阶顺序主子式为0')

end

end

%初始化,定义n阶方阵L，R

L=zeros(n,n);

R=zeros(n,n);

%定义L中的主对角线元素均为1

for i=1:n

L(i,i)=1;

end

%开始进行Doolittle分解，按照紧凑格式的公式编程

for k=1:n

for j=k:n

sum=0;

for t=1:k-1

sum=sum+L(k,t)*R(t,j);%求和

end

R(k,j)=A(k,j)-sum; %计算R的第k行元素

end

for i=k+1:n

sum=0;

for t=1:k-1

sum=sum+L(i,t)*R(t,k); %求和

L(i,k)=(A(i,k)-sum)/R(k,k); %计算L的第k列end

end

E=eye(n); %定义一个n阶单位矩阵

L=L*E %在Matlab中输出L矩阵

D=zeros(n,n); %定义D为n阶方阵

%构造D为对角矩阵

for i=1:n

D(i,i)=R(i,i);

end

D=D*E %在Matlab中输出D矩阵

R=inv(D)*R %求矩阵R