LU分解での連立方程式の解放

clear; close all;clc;
A =  [2 2 1;
    3 6 2;
    4 -2 5];
b = [4 2 8]';
[L, U] = lu_nopivot(A);
y = forward_subs(L,b)
[U2,y2] = forward_elim(U,y); %一度対角成分の係数をゼロにする必要あり
x = backward_subs(U2,y2);

function [L,U] = lu_nopivot(A)
n=size(A,1);
L=eye(n);
for k=1:n
    L(k+1:n,k) = A(k+1:n,k)/A(k,k);
    for l=k+1:n
        A(l,:) = A(l,:)-L(l,k)*A(k,:);
    end
end
U=A;
end

function x=forward_subs(A,b)
n=length(A);
x(1) = b(1);
for i=2:n
    res = 0;
    for j=1:i-1
        res = res + ...
               A(i,j) * x(j);
    end
    x(i) = b(i) - res;
end
end

function [A,b]=forward_elim(A,b)
n=length(A);
for i=1:n
    pivot = A(i,i);
    for j=1:n
        A(i,j) = A(i,j)/pivot;
    end
    b(i) = b(i) / pivot;
    for k=i+1:n
        pivot2 = A(k,i);
        for j=i:n
            A(k,j) = A(k,j)-...
                   A(i,j)*pivot2;
        end
        b(k) = b(k) - b(i)*pivot2;
    end
end
end

function x=backward_subs(A,b)
n=length(A);
x(n) = b(n);
for i=n-1:-1:1
    res = 0;
    for j=i+1:n
        res = res +...
               A(i,j) * x(j);
    end
    x(i) = b(i) - res;
end
end