matlab用三角分解法解函数

function x=sanjiao(A,b)
he1=0;
he2=0;

sizes = size(A)
for i=1:sizes(1)
    A(i,1)=A(i,1)/A(1,1)
end


%将矩阵A改成L和U矩阵
%理论上L和U是两个单独的矩阵，为了减少内存，直接在A矩阵上修改，用到的时候直接取A矩阵的值
for r=2:sizes(1)
    for i=r:sizes(2)
        for k=1:r-1
            he1=he1+A(r,k)*A(k,i)
        end
        A(r,i)=A(r,i)-he1;
        he1=0;
    end
    for i=r+1:sizes(2)
        for k=1:r-1
            he2=he2+A(i,k)*A(k,r)
        end
        A(i,r)=(A(i,r)-he2)/A(r,r);
        he2=0;
    end
end

%根据Ly=b,求解矩阵y
y=zeros(size(A,1),1);
x=zeros(size(A,1),1);
y(1)=b(1);
he=0;hehe=0;
for i=2:size(y,1)
    for j=1:i-1
        he=he+A(i,j)*y(j)
    end
    y(i)=b(i)-he
    he=0;
end

%根据Ux=y，求解x矩阵
x(size(A,1))=y(size(A,1))/A((size(A,1)),(size(A,1)))
for i=size(A,1)-1:-1:1
    for j=i+1:size(A,1)
        hehe=hehe+A(i,j)*x(j)
    end
    x(i)=(y(i)-hehe)/A(i,i)
    hehe=0
end


end