Crout分解法 | matlab

function [x, L, U]=Crout(A, b)
% Crout分解法求解线性方程组
% A为系数矩阵，b为常数项列
N = size(A);
n = N(1);
L = zeros(n, n); % 下三角矩阵
U = eye(n, n);   % 上三角矩阵
L(1:n, 1) = A(1:n, 1);    % L的第一列
U(1, 1:n) = A(1, 1:n) / L(1, 1);  % U的第一行
for k = 2:n
    for i = k:n
        % L的第k列
        L(i, k) = A(i, k) - L(i, 1:(k-1)) * U(1:(k-1), k);
    end
    for j = (k+1):n
        % U的第k行
        U(k, j) = (A(k, j) - L(k, 1:(k-1)) * U(1:(k-1), j)) / (L(k, k));
    end
end

y = SolveDownTriangle(L ,b);
x = SolveUpTriangle(U, y);

function x=SolveDownTriangle(A, b)
% 求解下三角矩阵Ax=b的解
N = size(A);
n = N(1);
for i = 1:n
    if i>n
        s = A(i, 1:(i-1)) * x(1:(i-1), 1);
    else
        s = 0;
    end 
    x(i, 1) = (b(i) - s) / A(i, i);
end

function x=SolveUpTriangle(A, b)
% 求解上三角矩阵Ax=b的解
N = size(A);
n = N(1);
for i = n:-1:1
    if (i<n)
        s = A(i, (i+1):n) * x((i+1):n, 1);
    else
        s = 0;
    end
    x(i, 1) = (b(i) - s) / A(i, i);
end

A = [3, 2, 0, 0; -1, 3, 2, 0; 0, -1, 3, 2; 0, 0, -1, 3];
b = [7; 11; 15; 9];
[x, L, U] = Crout(A, b)

%=> x = [1.3538, 1.4693, 2.8063, 2.5252]