为了和教材上列主元高斯消元法的求解过程一致，编写了MATLAB程序，供大家参考，互相交流学习。

先写出函数

function x = gsxylzy(A,b)

% format rat;clc

% A = [-3 2 6;10 -7 0; 5 -1 5];

% b = [4; 7 ;6];

M = [A,b]

n = length(b);

%%% 消元过程

for i = 1:n-1

    %交换列主元

    for k = i:n-1     

       ak = max(abs(M(k:n,k))); index = find(M(:,k)==ak);

       if isempty(index)

          index = find(M(:,k)==-ak);

       end

     temp = M(index,:);  M(index,:) = M(k,:); M(k,:) = temp;

    end

    fprintf(['第 ',num2str(i),' 次交换列主元\n'])

    M

     % 消元

    for k = i:n-1

        M(k+1,:) = M(k+1,:)-M(i,:)*M(k+1,i)/M(i,i);

    end

    fprintf(['第 ',num2str(i),' 列消元\n'])

    M

end 

%%% 回代过程

x = zeros(n,1); 

x(n) = M(n,n+1)/M(n,n);

fprintf(['第 ',num2str(n),' 个变量 x(',num2str(n),')=%s\n'],num2str(x(n)));

for k = (n-1):-1:1

    x(k) = (M(k,n+1)-M(k,k+1:n)*x(k+1:n))/M(k,k);

    fprintf(['第 ',num2str(k),' 个变量 x(',num2str(k),')=% s\n'],num2str(x(k)));

end

主程序测试应用：

format rat;clc

A = [-3 2 6;10 -7 0; 5 -1 5];

b = [4; 7 ;6];

x = gsxylzy(A,b);

fprintf('方程组的解为\n')

x

运行结果为

M =

      -3              2              6              4       

      10             -7              0              7       

       5             -1              5              6       

第 1 次交换列主元

M =

      10             -7              0              7       

      -3              2              6              4       

       5             -1              5              6       

第 1 列消元

M =

      10             -7              0              7       

       0             -1/10           6             61/10    

       0              5/2            5              5/2     

第 2 次交换列主元

M =

      10             -7              0              7       

       0              5/2            5              5/2     

       0             -1/10           6             61/10    

第 2 列消元

M =

      10             -7              0              7       

       0              5/2            5              5/2     

       0              0             31/5           31/5     

第 3 个变量 x(3)=1

第 2 个变量 x(2)=-1

第 1 个变量 x(1)=2.6645e-16

方程组的解为

x =

       1/3752999689475414 %这里严格说是0

      -1       

       1       

>> 

如果不想出现这种分数近似，也可以用符号形式，则不会出现近似值.

format rat;clc

A = sym([-3 2 6;10 -7 0; 5 -1 5]);

b = sym([4; 7 ;6]);

x = gsxylzy(A,b);

fprintf('方程组的解为\n')

x

则运行结果为

M =

[ -3,  2, 6, 4]

[ 10, -7, 0, 7]

[  5, -1, 5, 6]

第 1 次交换列主元

M =

[ 10, -7, 0, 7]

[ -3,  2, 6, 4]

[  5, -1, 5, 6]

第 1 列消元

M =

[ 10,    -7, 0,     7]

[  0, -1/10, 6, 61/10]

[  0,   5/2, 5,   5/2]

第 2 次交换列主元

M =

[ 10,    -7, 0,     7]

[  0,   5/2, 5,   5/2]

[  0, -1/10, 6, 61/10]

第 2 列消元

M =

[ 10,  -7,    0,    7]

[  0, 5/2,    5,  5/2]

[  0,   0, 31/5, 31/5]

第 3 个变量 x(3)=1

第 2 个变量 x(2)=-1

第 1 个变量 x(1)=0

方程组的解为

x =

       0       

      -1       

       1       

>> 

这里使用了符号变量命令 sym，结果和教材完全一样。