高斯消元法程序

贡献者： addis

预备知识　线性方程组高斯消元法

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

代码 1：GaussEli.m

% 高斯消元法得到梯形系数矩阵
% 显示每个行变换的操作和结果
% 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

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。