LU三角分解（matlab程序）

A=[1,2,3;1,3,5;1,3,6];
b=[2,3,4];
x=grout(A,b);




function x=grout(B,c)
n=size(B,1);
L=eye(n);
U=zeros(n);
for i=1:n
    s=0;
    t=0;
    for j=1:i-1
        s=s+L(i,j)*U(j,i:n);
        t=t+L(i+1:n,j)*U(j,i);
    end
   U(i,i:n)= B(i,i:n)-s;
   L(i+1:n,i)=(B(i+1:n,i)-t)/U(i,i);
end
y=grout1(L,c);
x=grout2(U,y);
 


function y=grout1(B,c)
n=size(B,1);
y=zeros(n,1);
for i=1:n
    s=0;
    for j=1:i-1
        s=s+B(i,j)*y(j);
    end
    y(i)=c(i)-s;
end
end




function x=grout2(U,y)
n=size(U,1);
x=zeros(n,1);
for i=n:-1:1
    s=0;
    for j=n:-1:i+1
        s=s+U(i,j)*y(j);
    end
    x(i)=(y(i)-s)/U(i,i);
end
end