LU分解法 | matlab

% LU分解法
% M为输入的增广矩阵
% precision为输入的精度要求，如不输入或输入有误，则默认为10位

if nargin == 2
    try
        digits(precision);
    catch
        disp('输入的精度有误，这里按照10位有效数字计算')
        digits(10);
    end
else
    digits(10);
end

A = vpa(M);
row = size(A, 1);
col = size(A, 2);
if ndims(A) ~= 2 | col - row ~= 1
    disp('矩阵大小有误，不能使用LU分解')
    return
end

if det(M(:,1:row)) == 0
    disp('该方程的系数矩阵行列式为零，不能使用LU分解法')
    return
end

 % 调用LU命令
[L, U, P]=lu(double(A))

% 回代求解
for i = row:-1:1
    temp = U(i, col);
    for k = i+1:row
        temp = vpa(temp - t_solution(k)*U(i, k));
    end
    t_solution(i) = vpa(temp / U(i, i));
end

for i = 1:row
    temp = t_solution(i);
    for k = 1:i-1
        temp = vpa(temp - t_solution(k)*U(i, k));
    end
    solution(i) = temp;
end

A = [0.001, 2.000, 3.000, 1.000; -1.000, 3.712, 4.623, 2.000; -2.000, 1.072, 5.643, 3.000];
Mlu(A, 4)

%=> x = [ -0.4904, -0.05104, 0.3675]