高斯消元法程序

贡献者： addis

　　 这里给出一个简单的程序演示高斯消元法的基本步骤。注意在实际应用中，我们解线性方程组一般使用 Matlab 的反斜杠算符：x=A \ y，其中 $ \boldsymbol{\mathbf{A}} $ 是系数矩阵，$ \boldsymbol{\mathbf{y}} $ 是常数列，$ \boldsymbol{\mathbf{x}} $ 是未知量。

% 高斯消元法得到梯形系数矩阵
% 显示每个行变换的操作和结果
% A 为任意尺寸的矩阵
% 如果 A 是增广矩阵， 输入 augmented = true， 如果 A 是系数矩阵， augmented = false
% 输出的 A 是梯形矩阵， q 是一个行矢量， q(ii) 是第 ii 行第一个不为零的列标
% 如果第 ii 行全为零， 则 q(ii) = 0

function [A, q] = GaussEli(A, augmented)
% 系数矩阵的行数 Nr 和列数 Nc
[Nr, Nc] = size(A);
if (augmented)
    Nc = Nc - 1;
end
q = zeros(1,Nr); q(1) = 1;
disp(A);
for ii = 1:Nr-1
    % 处理第 ii 行
    if (ii > 1)
        q(ii) = q(ii - 1) + 1;
    end
    for jj = q(ii) : Nc
        % 检查第 jj 列
        if (max(abs(A(ii:end, q(ii)))) > 0)
            % 消元使 A(ii+1:end, q(ii)) 全为 0
            A = eli1(A, ii, q(ii));
            break;
        else
            % 不需要消元，检查下一列
            if (q(ii) < Nc)
                q(ii) = q(ii) + 1;
            else
                return;
            end
        end
    end
end
end

% 做行变换使 A(ii+1:end, q) 全为 0
% A(ii:end,q) 不能全为 0
function A = eli1(A, ii, q)
% 交换两行，使 A(ii,q) 在 A(ii:end, q) 中最大
[~, iimax] = max(abs(A(ii:end, q)));
iimax = ii + iimax - 1;
if (iimax > ii)
    disp(['交换: r_', num2str(ii), ' <-> r_', num2str(iimax)]);
    temp = A(ii,:);
    A(ii,:) = A(iimax,:);
    A(iimax,:) = temp;
    disp(A);
end

% 用第 ii 行消去所有 A(ii+1:end, q)
for jj = ii+1 : size(A,1)
    if (A(jj,q) == 0)
        continue;
    end
    k = -A(jj,q)/A(ii,q);
    disp(['消元: ', num2str(k), ' * r_', num2str(ii), ' + r_', num2str(jj)]);
    A(jj,:) = A(ii,:) * k + A(jj,:);
    disp(A);
end
end

　　 与 “高斯消元法” 中的步骤略有不同的是，该程序在处理每一行 A(ii,:) 时都会试图做一个行交换使系数 A(ii,q(ii)) 的绝对值尽可能大。这是为了减小数值误差：试想如果 A(ii,q(ii)) 的解析值为 0，但由于数值误差，计算出来是一个很小的数（例如 1e-15），那么用其消元时就可能需要将第 ii 行乘以一个很大的系数（例如 1e15）再加到另一行上，导致程序不稳定。

　　 用例 2 中的增广矩阵测试程序如下：

>> A = [1 1 -1 1; 2 2 -2 1; 1 1 0 2; 2 2 -1 5];
>> y = [3; 7; 3; 4];
>> [A1, q] =GaussEli([A,y],true); q
     1     1    -1     1     3
     2     2    -2     1     7
     1     1     0     2     3
     2     2    -1     5     4

交换: r_1 <-> r_2
     2     2    -2     1     7
     1     1    -1     1     3
     1     1     0     2     3
     2     2    -1     5     4

消元: -0.5 * r_1 + r_2
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
    1.0000    1.0000         0    2.0000    3.0000
    2.0000    2.0000   -1.0000    5.0000    4.0000

消元: -0.5 * r_1 + r_3
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    1.5000   -0.5000
    2.0000    2.0000   -1.0000    5.0000    4.0000

消元: -1 * r_1 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    1.5000   -0.5000
         0         0    1.0000    4.0000   -3.0000

交换: r_2 <-> r_3
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    0.5000   -0.5000
         0         0    1.0000    4.0000   -3.0000

消元: -1 * r_2 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    0.5000   -0.5000
         0         0         0    2.5000   -2.5000

交换: r_3 <-> r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    2.5000   -2.5000
         0         0         0    0.5000   -0.5000

消元: -0.2 * r_3 + r_4
    2.0000    2.0000   -2.0000    1.0000    7.0000
         0         0    1.0000    1.5000   -0.5000
         0         0         0    2.5000   -2.5000
         0         0         0         0         0
         
q =
     1     3     4     0