%直接三角分解法
function my_LU(a, b)
n = length(a);
l = zeros(n, n);%初始化
u = zeros(n, n);
for i=1:n
l(i,i) = 1;
end
u(1,1:n) = a(1,1:n);
l(2:n, 1) = a(2:n, 1) ./ u(1,1);
for r=2:n
for i=r:n
u(r, i) = a(r, i) - sum(l(r,1:r-1) .* (u(1:r-1,i))');
end
for i=r+1:n
if (r~=n)
l(i, r) = (a(i, r) - sum(l(i,1:r-1) .* (u(1:r-1,r)))')./u(r,r);
end
end
end
L=l,U=u%输出LU矩阵
y(1) = b(1);
for i=2:n
y(i) = b(i) - sum(l(i, 1:i-1).*y(1:i-1));
end
x(n) = y(n)/u(n,n);
for i=n-1:-1:1
x(i) = (y(i) - sum(u(i,i+1:n).*x(i+1:n)))./u(i,i);
end
x=x'
end