MATLAB - 高斯-约旦消元法

高斯-约旦消元法（Gauss-Jordan elimination）是线性代数中用于求解线性方程组和求矩阵逆的一个基本算法。该方法系统地应用初等行运算将矩阵转换为其 简化行阶梯形（RREF）。虽然理解此算法对于掌握线性代数概念至关重要，但 MATLAB 提供了高度优化、数值稳定且高效的工具来完成这些任务，在实际应用中应优先选用。

本教程将涵盖两种方法：

MATLAB 的方式： 使用内置的 rref 函数和反斜杠运算符（\）来获得快速可靠的解。这是实际应用的最佳实践。
算法深入探讨： 逐步实现高斯-约旦算法以理解其底层机制。

问题：一个线性方程组

让我们考虑以下线性方程组：

$$\mathrm{\begin{cases} x:+:y:+:2z:=:9 \ 2x:+:4y:-:3z:=:1 \ 3x:+:6y:-:5z:=:0\end{cases}}$$

我们可以将此系统表示为矩阵形式 Ax = b，其中 A 是系数矩阵，x 是变量向量，b 是常数向量。

% 系数矩阵 A
A = [1 1 2; 2 4 -3; 3 6 -5];

% 常数向量 b
b = [9; 1; 0];

% 增广矩阵结合了 A 和 b
augmentedMatrix = [A, b];
% 即 [1 1 2 9; 2 4 -3 1; 3 6 -5 0]

第一部分：MATLAB 的方式（最佳实践）

对于实际问题解决，始终使用 MATLAB 内置的、高度优化的函数。

方法 1：使用 `rref` 函数

rref 函数直接计算矩阵的简化行阶梯形。这是在 MATLAB 中执行高斯-约旦消元法最直接的方式。

% 定义增广矩阵
augmentedMatrix = [1 1 2 9; 2 4 -3 1; 3 6 -5 0];

% 计算简化行阶梯形
R = rref(augmentedMatrix);

% R 的最后一列是解
solution = R(:, end);

disp('简化行阶梯形为：');
disp(R);

disp('解向量 [x; y; z] 为：');
disp(solution);

输出将显示左侧是一个 3x3 的单位矩阵，最后一列是解向量 [1; 2; 3]，这意味着 x=1，y=2，z=3。

方法 2：使用反斜杠运算符（`\`）

当您只需要 Ax = b 的解而不是 RREF 时，反斜杠运算符（\），也称为 mldivide，是 MATLAB 中最有效且数值鲁棒的方法。它使用一组不同的、更高级的算法（如 LU 分解），这些算法经过优化以提高速度和准确性。

% 定义系数矩阵和常数向量
A = [1 1 2; 2 4 -3; 3 6 -5];
b = [9; 1; 0];

% 使用反斜杠运算符求解 x
solution = A \ b;

disp('解向量 [x; y; z] 为：');
disp(solution);

第二部分：算法深入探讨

现在，让我们从头开始实现高斯-约旦算法，以了解其工作原理。这仅用于教育目的，不应在生产代码中使用，因为 rref 或 \ 更为适用。

该过程涉及两个主要阶段：前向消元（在对角线下方创建零）和后向消元（在对角线上方创建零），以及用于数值稳定性的主元选择。

function solution = solveWithGaussJordan(A, b)
    % SOLVEWITHGAUSSJORDAN 使用手动实现的高斯-约旦消元法求解线性方程组 Ax=b。
    % 它包含用于数值稳定性的部分主元选择。
    % 返回解向量 x。

    % 构建增广矩阵
    [numRows, numCols] = size(A);
    if numRows ~= numCols
        error('系数矩阵 A 必须是方阵。');
    end
    augmentedMatrix = [A, b];
    n = numRows;

    % --- 主要的高斯-约旦消元循环 ---
    for i = 1:n
        % 1. 部分主元选择
        % 在当前列中找到具有最大主元元素的行
        [~, maxRowIndex] = max(abs(augmentedMatrix(i:n, i)));
        maxRowIndex = maxRowIndex + i - 1; % 调整索引使其成为绝对值

        % 将当前行与包含最大主元的行进行交换
        if maxRowIndex ~= i
            augmentedMatrix([i, maxRowIndex], :) = augmentedMatrix([maxRowIndex, i], :);
        end

        % 检查奇异性
        if abs(augmentedMatrix(i, i)) < 1e-10 % 使用容差进行浮点比较
            error('矩阵奇异或接近奇异。无法求解系统。');
        end

        % 2. 规范化主元行
        % 将主元行除以主元元素，使主元元素变为 1
        augmentedMatrix(i, :) = augmentedMatrix(i, :) / augmentedMatrix(i, i);

        % 3. 消除当前列中的其他项
        for j = 1:n
            if i ~= j
                % 从所有其他行中减去主元行的某个倍数
                factor = augmentedMatrix(j, i);
                augmentedMatrix(j, :) = augmentedMatrix(j, :) - factor * augmentedMatrix(i, :);
            end
        end
    end

    % 从最后一列提取解
    solution = augmentedMatrix(:, end);

    % --- 验证 (最佳实践) ---
    residual = A * solution - b;
    if norm(residual) > 1e-6
        warning('解可能不准确。残差范数: %e', norm(residual));
    end
end

% --- 如何使用该函数 ---
% A = [1 1 2; 2 4 -3; 3 6 -5];
% b = [9; 1; 0];
% mySolution = solveWithGaussJordan(A, b);

代码解释

函数定义： 代码被封装在 solveWithGaussJordan 函数中，以实现复用性和清晰度。
输入验证： 它首先检查矩阵 A 是否为方阵，这是此简单求解器的要求。
部分主元选择： 对于每一列 i，它会找到其下方具有最大绝对值的行。该行与当前行 i 进行交换。这种称为部分主元选择（partial pivoting）的技术对于最小化浮点误差和避免除以零至关重要。
奇异性检查： 在主元选择后，它会检查主元元素是否接近零。如果是，则矩阵是奇异的（或接近奇异），不存在唯一解。程序会抛出错误。
规范化： 主元行除以主元元素 A(i,i)，使新的主元元素等于 1。
消元： 对于所有其他行 j，它会减去（新规范化的）主元行的某个倍数，以使元素 A(j,i) 变为零。
解的提取： 循环完成后，矩阵的系数部分是一个单位矩阵，最后一列包含解向量。
验证： 一个好的实践是检查解。残差（A*x - b）应该是一个零向量。我们检查它的 范数（大小）以确定其是否可接受地小。

结论与最佳实践

尽管实现高斯-约旦消元法等算法是一个很好的学习练习，但在专业或学术环境中，为工作选择合适的工具至关重要。

求解 Ax = b： 始终首选反斜杠运算符：x = A \ b。
查找 RREF： 使用内置函数：R = rref(augmentedMatrix)。
数值稳定性： MATLAB 的内置方法经过专业开发和测试，能够比简单的脚本更好地处理各种数值挑战。
性能： 内置函数以低级语言（如 C++ 或 Fortran）实现，并且比 for 循环中的解释性 MATLAB 代码快得多。

MATLAB - 高斯-约旦消元法

MATLAB - 高斯-约旦消元法

问题：一个线性方程组

第一部分：MATLAB 的方式（最佳实践）

方法 1：使用 rref 函数

方法 2：使用反斜杠运算符（\）

第二部分：算法深入探讨

代码解释

结论与最佳实践

方法 1：使用 `rref` 函数

方法 2：使用反斜杠运算符（`\`）