lufact.m

来自「余向阳博士Matlab教程的源代码。 包括：1.非线性方程的解法；2.线性方程」· M 代码 · 共 36 行

function X=lufact(A,b)
%  A为n阶矩阵。
%  b是n维向量。
%  X是所求的AX=b的解。
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
C=zeros(1,N);
R=1:N;
for p=1:N-1
[max1,j]=max(abs(A(p:N,p)));
C=A(p,:);
A(p,:)=A(j+p-1,:);
A(j+p-1,:)=C;
d=R(p);
R(p)=R(j+p-1);
R(j+p-1)=d;
if A(p,p)==0
    'A是奇异阵，方程组无惟一解'
    break
end
for k=p+1:N
   mult=A(k,p)/A(p,p);
   A(k,p)=mult;
   A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);
end
end
Y(1)=b(R(1));
for k=2:N
    Y(k)=b(R(k))-A(k,1:k-1)*Y(1:k-1);
end
X(N)=Y(N)/A(N,N);
for k=N-1:-1:1
    X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end