function [L,U]=Crout(A)% 矩阵的三角分解法之LU分解之Crout分解 A=LU% Crout分解：LU分解中L为下三角阵，U为单位上三角阵% 说明：n阶(n>=2)矩阵A有惟一Crout分解的充要条件是：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]=Crout(A)%% L =%% 2 0 0% 4 3 0% 6 -6 -1%%% U =%% 1.0000 0.5000 1.0000% 0 1.0000 0% 0 0 1.0000%% >> L*U-A%% ans =%% 0 0 0% 0 0 0% 0 0 0
[n,m]=size(A);if n~=m error('请输入方阵');endU=eye(n); % 初始化U为单位矩阵，即根据Crout分解的特性待定U的对角线上的元素为1L=zeros(n); % 初始化L为全零矩阵
for i=1:n L(i,1)=A(i,1);endfor j=2:n U(1,j)=A(1,j)/L(1,1);endfor k=2:n for j=k:n L(j,k)=A(j,k)-L(j,1:k-1)*U(1:k-1,k); end for i=k:n U(k,i)=(A(k,i)-L(k,1:k-1)*U(1:k-1,i))/L(k,k); endend