crout.m

来自「自己编写的程序」· M 代码 · 共 30 行

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