gauss.m

数值分析中几个基础的matlab程序

function[x,det,flag]=Gauss(A,b)
%求线型方程组的列主元Gauss消去法，其中，
%A为方程组的系数矩阵；
%b为方程组的右端项；
%x为方程组的解；
%det为系数矩阵A的行列式的值；
%flag为指标向量，flag='OK'表示计算成功。
[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
%开始计算，先赋初值
flag='OK';det=1;x=zeros(n,1);
for k=1:n-1
%选主元
    max1=0
    for i=k:n
        if abs(A(i,k))>max1
           max1=abs(A(i,k));r=i;
        end
    end
    if max1<1e-10
       flag='failure';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
    flag='failure';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
x(k)=b(k)/A(k,k);
end