function [L,U]=Crout(A)
%  矩陣的三角分解法之LU分解之Crout分解  A=LU
%  Crout分解：LU分解中L為下三角陣，U為單位上三角陣 
%  說明：n階(n>=2)矩陣A有惟一Crout分解的充要條件是：A的前n-1個順序主子式不為0
%  編程思想：使用待定系數法
[n,m]=size(A);
if n~=m    
    error('請輸入方陣');
end
U=eye(n); 
% 初始化U為單位矩陣，即根據Crout分解的特性待定U的對角線上的元素為1 L=zeros(n);
% 初始化L為全零矩陣

for i=1:n     
    L(i,1)=A(i,1); 
end
for j=2:n    
    U(1,j)=A(1,j)/L(1,1); 
end 
for 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);  
    end
end