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

>>
列印頁面
上一頁
下一頁
廣告