【Matlab】矩阵三角分解法求解方程组

%直接三角分解法
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