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]