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