lucourant.m

来自「这是线性议程组的直接解法」· M 代码 · 共 29 行

function [L,U,x]=LUCourant(A,b)
n=length(A);
U=eye(n);
for k=1:n
    if k==1
        L(k:n,k)=A(k:n);
        U(k,k+1:n)=A(k,k+1:n)/L(k,k);
    else
        L(k:n,k)=A(k:n,k)-L(k:n,1:k-1)*U(1:k-1,k);
        if k~=n
            U(k,k+1:n)=(A(k,k+1:n)-L(k,1:k-1)*U(1:k-1,k+1:n))/L(k,k);
        end
    end
end
for k=1:n
    if k==1
        y(k,1)=b(k,1)/L(k,k);
    else
        y(k,1)=(b(k,1)-L(k,1:k-1)*y(1:k-1,1))/L(k,k);
    end
end
for k=n:-1:1
    if k==n
        x(k,1)=y(k,1);
    else
        x(k,1)=y(k,1)-U(k,k+1:n)*x(k+1:n,1);
    end
end
end