gaussxiaoqufa.m

数值分析中的应用程序

function [x,det,biaozhi]=Gaussxiaoqufa(A,b)
% 求线性方程组的列主元Gauss消去法，其中
%   A --- 方程组的系数矩阵
%   b --- 方程组的右端项
%   x --- 方程组的解
%   det --- 系数矩阵A的行列式的值
%   biaozhi --- biaozhi=0表示计算成败；biaozhi=1表示计算成功
format long
[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 dimension of b!');
    return;
end
biaozhi=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-20
        biaozhi=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-20
    biaozhi=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