线性代数方程组数值解法：直接三角分解法（LU分解）

本文通过matlab实现

如若需要输出L矩阵和U矩阵在代码后添加disp('L')disp('U')即可。

%线性方程组求解（直接法LU分解）
function linearequation(A,B) %A为系数矩阵,B为数值矩阵
i=size(A);
d=size(B);
L=zeros(i);
U=L; %设置LU矩阵
x=zeros(i(1),1);
y=zeros(i(1),1); %设置xy列向量
for m=1:i(1)
    L(m,m)=1;
end
if i(1)~=i(2)|d(1)~=i(1)|d(2)~=1 %判断A矩阵是否为方阵B是否为列向量
    b='输入矩阵有误';
    disp(b);
else
    %计算LU矩阵

        for o=1:i(1)%计算U的第一行
            U(1,o)=A(1,o);
        end

        for o=1:i(1)%计算L的第一列
            L(o,1)=A(o,1)/U(1,1);
        end

    for m=2:i(1) %计算LU的每一列行

        for r=2:i(1) %计算U的行
            c=0;
            for k=1:m-1
                c=c+L(m,k)*U(k,r);
            end
            U(m,r)=A(m,r)-c;
        end 

        for l=2:i(1) %计算L的列
            c=0;
            for k=1:m-1
                c=c+L(l,k)*U(k,m);
            end
            L(l,m)=(A(l,m)-c)/U(m,m);
        end

    end
    
    %求解下三角方程组Ly=B
    y(1)=B(1);
    for k=2:i(1)
        c=0;
        for t=1:k-1
            c=c+L(k,t)*y(t);
        end
        y(k)=B(k)-c;
    end
    
    %求解上三角方程组Ux=y
    x(i(1))=y(i(1))/U(i(1),i(1));
    for k=(i(1)-1):-1:1
        c=0;
        for t=k+1:i(1)
            c=c+U(k,t)*x(t);
        end
        x(k)=(y(k)-c)/U(k,k);
    end

    disp(x); %输出x
  
end