har_linsolver_Crout2.m

function x=har_linsolver_Crout2%(A,b)
%=========================================================
% 用克洛脱分解求解线性方程组
% 测试用例
A=[1 2 3 4; 1 4 9 16; 1 8 27 64; 1 16 81 256 ];
b=[2 10 44 190]';
%=========================================================

n=length(b);

%================分解=====================================
for i=1:n
    %计算L的第i列
    if i==1
        A(:,i)=A(:,i);
    else
        for j=i:n
            A(j,i)=A(j,i)-A(j,1:i-1)*A(1:i-1,i);
        end
    end
    %计算U的第i行
    for j=i+1:n
        A(i,j)=(A(i,j)-A(i,1:i-1)*A(1:i-1,j))/A(i,i);
    end
end
disp(A);
% disp(L);
% disp(U);

%========================回代===============================
%求解中间向量
y=zeros(n,1);
y(1)=b(1)/A(1,1);
for i=2:n
    y(i)=(b(i)-A(i,1:i-1)*y(1:i-1))/A(i,i);
end
%求解向量
x=zeros(n,1);
x(n)=y(n);
for i=n-1:-1:1
    x(i)=y(i)-A(i,i+1:n)*x(i+1:n);
end

end