crout.m

用Crout分解法分解系数矩阵

function [L U]=crout(A) 
n=size(A);
L=zeros(n,n);U=zeros(n,n);
for i=1:n
    U(i,i)=1;%U矩阵的主对角上都是1
end
for i=2:n
    L(i,i-1)=A(i,i-1);%即ri=bi
end
U(1,2)=A(1,2)/A(1,1);
for i=2:n-1
    U(i,i+1)=A(i,i+1)/(A(i,i)-A(i,i-1)*U(i-1,i));
end
L(1,1)=A(1,1);
for i=2:n
    L(i,i)=A(i,i)-A(i,i-1)*U(i-1,i);
end