MATLAB - 高斯-约旦消元法

高斯-约旦消元法是一种系统的方法，用于求解线性方程组、求矩阵的逆以及执行其他与矩阵相关的计算。该技术是高斯消元法的扩展，旨在通过一系列行操作将给定矩阵转换为其简化行阶梯形式 (RREF)。

高斯-约旦消元法广泛应用于工程、物理、计算机科学和经济学等各个领域。其应用包括：

求解线性方程组。
计算矩阵的逆。
确定矩阵的秩。
求解线性规划问题。

在 MATLAB 中，可以使用内置函数或手动应用必要的行操作来实现高斯-约旦消元法，以达到所需的格式。主要步骤包括：

**向前消元** - 通过使用行操作在枢轴位置下方创建零来将矩阵转换为上三角形式。
**向后消元** - 执行进一步的行操作以在枢轴位置上方创建零，从而得到简化行阶梯形式。

当每个非零行的首项为 1，并且包含首项的列中的所有元素都为零（首项本身除外）时，矩阵被称为简化行阶梯形式。

以下是 MATLAB 中高斯-约旦消元过程的分步细分。

**初始化矩阵** - 定义增广矩阵，该矩阵将系数矩阵和线性方程中常数结合在一起。
**执行行操作** - 使用 MATLAB 命令执行行交换、缩放和加法，以便在适当的位置创建零。
**转换为 RREF** - 继续此过程，直到矩阵处于简化行阶梯形式。
**提取解** - 最终矩阵形式将直接提供线性方程组的解。

高斯-约旦消元法用于求解线性方程组并求矩阵的逆。该方法使用初等行操作将矩阵转换为其简化行阶梯形式 (RREF)。

让我们使用高斯-约旦消元法求解一个线性方程组。考虑以下系统：

$$\mathrm{\begin{cases} x\:+\:2y\:+\:3z\:=\:9 \\ 2x\:+\:3y\:+\:4z\:=\:13 \\ 3x\:+\:4y\:+\:5z\:=\:17\end{cases}}$$

以下是使用高斯-约旦消元法求解线性方程的步骤。

步骤 1 - 形成增广矩阵。

因此，我们从上述线性方程得到的增广矩阵如下所示。

$$\mathrm{A \: = \: \begin{bmatrix} 1 & 2 & 3 & | & 9 \\ 2 & 3 & 4 & | & 13 \\ 3 & 4 & 5 & | & 17\end{bmatrix}}$$

步骤 2 - 执行初等行操作以使矩阵处于 RREF。

让我们在 MATLAB 中实现这一点。

% Define the augmented matrix
A = [1 2 3 9; 2 3 4 13; 3 4 5 17];

% Get the number of rows
n = size(A, 1);

% Gauss-Jordan Elimination with partial pivoting
for i = 1:n
    % Pivoting: Swap rows if necessary to make the largest absolute value in the
    % current column the pivot element (for numerical stability)
    [~, pivotRow] = max(abs(A(i:n, i)));
    pivotRow = pivotRow + i - 1;
    if pivotRow ~= i
        A([i, pivotRow], :) = A([pivotRow, i], :);
    end
    
    % Make the diagonal element 1
    A(i, :) = A(i, :) / A(i, i);
    
    % Make the other elements in the current column 0
    for j = 1:n
        if i ~= j
            A(j, :) = A(j, :) - A(j, i) * A(i, :);
        end
    end
end

% Extract the solution
solution = A(:, end);

% Display the result
disp('The solution is:');
disp(solution);

以下是 Matlab 代码的解释。

步骤 1 - 定义增广矩阵。

A = [1 2 3 9; 2 3 4 13; 3 4 5 17];

此矩阵将变量的系数和方程右侧的常数结合在一起。

步骤 2 - 获取行数。

n = size(A, 1);

我们确定行数（也等于方程数）。

步骤 3 - 高斯-约旦消元循环。

% Gauss-Jordan Elimination with partial pivoting
for i = 1:n
    % Pivoting: Swap rows if necessary to make the largest absolute value in the
    % current column the pivot element (for numerical stability)
    [~, pivotRow] = max(abs(A(i:n, i)));
    pivotRow = pivotRow + i - 1;
    if pivotRow ~= i
        A([i, pivotRow], :) = A([pivotRow, i], :);
    end
    
    % Make the diagonal element 1
    A(i, :) = A(i, :) / A(i, i);
    
    % Make the other elements in the current column 0
    for j = 1:n
        if i ~= j
            A(j, :) = A(j, :) - A(j, i) * A(i, :);
        end
    end
end

代码：for i = 1:n，此处循环迭代矩阵中的每一行 i。目标是使系数矩阵的对角元素（即元素 A(i, i)）等于 1，并将当前列 i 中的其他元素设为 0。

[~, pivotRow] = max(abs(A(i:n, i)));
pivotRow = pivotRow + i - 1;
if pivotRow ~= i
    A([i, pivotRow], :) = A([pivotRow, i], :);
end

此代码块查找当前列 i 中（从第 i 行到最后一行）绝对值最大的行。这样做是为了减少数值不稳定性。绝对值最大的行成为枢轴行。

如果最大元素不在当前行，则将其与 pivotRow 交换。这确保了使用最大元素作为枢轴。

A(i, :) = A(i, :) / A(i, i);

这将第 i 行的整个内容除以对角元素 A(i, i)，使此对角元素为 1。这是将矩阵简化为行阶梯形式的重要一步。

for j = 1:n
    if i ~= j
        A(j, :) = A(j, :) - A(j, i) * A(i, :);
    end
end

此循环迭代所有行 j，除了当前行 i。对于每一行 j，它从行 j 中减去行 i 的倍数，以使元素 A(j, i) 为零。此操作将矩阵转换为行简化阶梯形式，其中当前列中的所有元素（枢轴（对角元素）除外）都变为零。

步骤 4 - 提取解

solution = A(:, end);

高斯-约旦消元完成后，矩阵 A 的最后一列包含方程组的解。此步骤将最后一列提取为解向量。

步骤 5 - 显示结果。

disp('The solution is:');
disp(solution);

最后，显示解。这将显示满足线性方程组的 x、y 和 z 的值。

当我们在 matlab 命令窗口中执行代码时，得到的输出为：

>> % Define the augmented matrix
A = [1 2 3 9; 2 3 4 13; 3 4 5 17];

% Get the number of rows
n = size(A, 1);

% Gauss-Jordan Elimination with partial pivoting
for i = 1:n
    % Pivoting: Swap rows if necessary to make the largest absolute value in the
    % current column the pivot element (for numerical stability)
    [~, pivotRow] = max(abs(A(i:n, i)));
    pivotRow = pivotRow + i - 1;
    if pivotRow ~= i
        A([i, pivotRow], :) = A([pivotRow, i], :);
    end
    
    % Make the diagonal element 1
    A(i, :) = A(i, :) / A(i, i);
    
    % Make the other elements in the current column 0
    for j = 1:n
        if i ~= j
            A(j, :) = A(j, :) - A(j, i) * A(i, :);
        end
    end
end

% Extract the solution
solution = A(:, end);

% Display the result
disp('The solution is:');
disp(solution);

The solution is:
    2.3333
   -1.6667
    3.3333

>>

打印页面