矩阵的三角分解法之LU分解之Doolittle分解

function [L,U]=Doolittle(A)
% 矩阵的三角分解法之LU分解之Doolittle分解 A=LU
% Doolittle分解：LU分解中L为单位下三角阵，U为上三角阵
% 说明：n阶(n>=2)矩阵A有惟一杜里特尔分解的充要条件是：A的前n-1个顺序主子式不为0
% 编程思想：使用待定系数法
%
% 作者：野渡无人
% 最后修改日期：2008.4.9
%
% >> A=[2 1 2;4 5 4;6 -3 5]
%
% A =
%
% 2 1 2
% 4 5 4
% 6 -3 5
%
% >> [L,U]=Doolittle(A)
%
% L =
%
% 1 0 0
% 2 1 0
% 3 -2 1
%
%
% U =
%
% 2 1 2
% 0 3 0
% 0 0 -1
%
% >> L*U-A
%
% ans =
%
% 0 0 0
% 0 0 0
% 0 0 0

[n,m]=size(A);
if n~=m
error('请输入方阵');
end
L=eye(n); % 初始化L为单位矩阵，即根据Doolittle分解的特性待定L的对角线上的元素为1
U=zeros(n); % 初始化U为全零矩阵
for i=1:n
U(1,i)=A(1,i);
end
for j=2:n
L(j,1)=A(j,1)/U(1,1);
end
for k=2:n
for j=k:n
U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);
end
for i=k+1:n
L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k))/U(k,k);
end
end