电力系统算法

zibuyu_2017

function [L,U,y,x]=LU(A,B)
n=length(A);
L=eye(n);
U=zeros(n,n);
for k=1:n
  for j=k:n
      c=0;
      for i=1:k-1
          c=c+L(k,i)*U(i,j);
      end
      U(k,j)=A(k,j)-c;
  end
  if U(k,k)==0
      for i=k+1:n
          z=i;
          if U(k,i)~=0
              break;end
      end
      for i=1:k
          a=U(i,k);U(i,k)=U(i,z);U(i,z)=a;
      end
      for i=1:n
          a=A(i,k);A(i,k)=A(i,z);A(i,z)=a;
      end
  end
  for i=k+1:n
      c=0;
      for q=1:k-1
          c=c+L(i,q)*U(q,k);
      end
      L(i,k)=(A(i,k)-c)/U(k,k);
  end
end
  for k=1:n
      c=0;
      for i=1:k-1
          c=c+L(k,i)*y(i);
      end
      y(k)=(B(k)-c)/L(k,k);
  end
  for k=n:-1:1
      c=0;
      for i=k+1:n
          c=c+U(k,i)*x(i);
      end
      x(k)=(y(k)-c)/U(k,k);
  end

zibuyu_2017

看到这帖子真是高兴！