gauss.m

一个自己写的列主元高斯消去法

function [x,det,index]=Gauss(A,b)
% 求线形方程组的列主元Gauss消去法，其中，
% A为方程组的系数矩阵；
% b为方程组的右端项；
% x为方程组的解；
% det为系数矩阵A的行列式的值；
% index为指标变量，index=0表示计算失败，index=1表示计算成功。
[n,m]=size(A); nb=length(b);
% 当方程组行与列的维数不相等时，停止计算，并输出出错信息。
if n~=m
    error('The rows and columns of matrix A must be equal!');
    return;
end
% 当方程组与右端项的维数不匹配时，停止计算，并输出出错信息
if m~=nb
    error('The columns of A must be equal the length of b!');
    return;
end
% 开始计算，先赋初值
index=1;det=1;x=zeros(n,1);
for k=1:n-1
    % 选主元
    a_max=0;
    for i=k:n
        if abs(A(i,k))>a_max
            a_max=abs(A(i,k));r=i;
        end
    end
    if a_max<1e-10
        index=0;return;
    end
    % 交换两行
    if r>k
        for j=k:n
            z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
        end
        z=b(k);b(k)=b(r);b(r)=z;det=-det;
    end
    % 消元过程
    for i=k+1:n
        m=A(i,k)/A(k,k);
        for j=k+1:n
            A(i,j)=A(i,j)-m*A(k,j);
        end
        b(i)=b(i)-m*b(k);
    end
    det=det*A(k,k);
end
det=det*A(n,n);
% 回代过程
if abs(A(n,n))<1e-10
    index=0;return;
end
for k=n:-1:1
    for j=k+1:n
        b(k)=b(k)-A(k,j)*x(j);
    end
    x(k)=b(k)/A(k,k);
end