col_gauss.m

来自「这是用来解线性方程的」· M 代码 · 共 42 行

function [x,L,U,P]=col_gauss(A,b,EPS)
    % ??GUASS?????¨,?????????ó???í列GUASS消元法，用增广矩阵处理
     if nargin<3
        EPS=50*eps;
     end
     n=length(b);
     P=eye(n);% 
     x=zeros(n,1);
     b=b(:);
     A=[A b];
     for k=1:n-1
     [amax,m]=max(abs(A(k:n,k)));%?ó?÷??求主元
     m=k+m-1;
     if k~=m
        A([k m],:)=A([m k],:);
        P([k m],:)=P([m k],:);%?????ó置换阵
     end 
     if abs(A(k,k))<=EPS
       error('det(A)=0,???????¨????!');
      return;   
     end
     for i=k+1:n
     A(i,k)=A(i,k)/A(k,k); 
     A(i,k+1:n+1)=A(i,k+1:n+1)-A(i,k)*A(k,k+1:n+1);
     end
    %??????????????:上述循环可改为：
    %A(k+1:n,k)=A(k+1:n,k)/A(k,k);
    %A(k+1:n,k+1:n+1)=A(k+1:n,k+1:n+1)-A(k+1:n,k)*A(k,k+1:n+1);
 
    end
      if abs(A(n,n))<=EPS
        error('det(A)=0,???????¨????!');
      end
   x(n)=A(n,n+1)/A(n,n);
   for i=n-1:-1:1
       x(i)=(A(i,n+1)-A(i,i+1:n)*x(i+1:n))/A(i,i);
   end
   A=A(:,1:n);
   L=eye(n)+tril(A,-1);
   U=triu(A);