MATLAB - Gauss-Jordan Elimination

MATLAB - Gauss-Jordan Elimination

Gauss-Jordan elimination is a fundamental algorithm in linear algebra for solving systems of linear equations and finding matrix inverses. The method systematically applies elementary row operations to transform a matrix into its reduced row echelon form (RREF). While understanding this algorithm is crucial for grasping linear algebra concepts, MATLAB provides highly optimized, numerically stable, and efficient tools for these tasks that should be preferred in practice.

This tutorial will cover two approaches:

• The MATLAB Way: Using the built-in rref function and the backslash operator (\) for fast and reliable solutions. This is the best practice for real-world applications.
• The Algorithmic Deep Dive: A step-by-step implementation of the Gauss-Jordan algorithm to understand the underlying mechanics.

The Problem: A System of Linear Equations

Let’s consider the following system of linear equations:

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

We can represent this system in matrix form, Ax = b, where A is the coefficient matrix, x is the vector of variables, and b is the vector of constants.

% Coefficient matrix A
A = [1 1 2; 2 4 -3; 3 6 -5];

% Constants vector b
b = [9; 1; 0];

% The augmented matrix combines A and b
augmentedMatrix = [A, b];
% which is [1 1 2 9; 2 4 -3 1; 3 6 -5 0]

Part 1: The MATLAB Way (Best Practice)

For practical problem-solving, always use MATLAB’s built-in, highly optimized functions.

Method 1: Using the rref Function

The rref function directly computes the reduced row echelon form of a matrix. This is the most direct way to perform Gauss-Jordan elimination in MATLAB.

% Define the augmented matrix
augmentedMatrix = [1 1 2 9; 2 4 -3 1; 3 6 -5 0];

% Calculate the reduced row echelon form
R = rref(augmentedMatrix);

% The last column of R is the solution
solution = R(:, end);

disp('The reduced row echelon form is:');
disp(R);

disp('The solution vector [x; y; z] is:');
disp(solution);

The output will show a 3x3 identity matrix on the left, with the solution vector [1; 2; 3] in the last column, meaning x=1, y=2, and z=3.

Method 2: Using the Backslash Operator ()

When you only need the solution to Ax = b and not the RREF, the backslash operator (\), also known as mldivide, is the most efficient and numerically robust method in MATLAB. It uses a different, more advanced set of algorithms (like LU decomposition) optimized for speed and accuracy.

% Define the coefficient matrix and the constant vector
A = [1 1 2; 2 4 -3; 3 6 -5];
b = [9; 1; 0];

% Solve for x using the backslash operator
solution = A \ b;

disp('The solution vector [x; y; z] is:');
disp(solution);

Part 2: The Algorithmic Deep Dive

Now, let’s implement the Gauss-Jordan algorithm from scratch to understand how it works. This is for educational purposes and should not be used in production code where rref or \ are available.

The process involves two main phases: Forward Elimination (creating zeros below the diagonal) and Backward Elimination (creating zeros above the diagonal), along with Pivoting for numerical stability.

function solution = solveWithGaussJordan(A, b)
    % SOLVEWITHGAUSSJORDAN Solves a system of linear equations Ax=b
    % using a manual implementation of Gauss-Jordan elimination.
    % It includes partial pivoting for numerical stability.
    % Returns the solution vector x.

    % Form the augmented matrix
    [numRows, numCols] = size(A);
    if numRows ~= numCols
        error('Coefficient matrix A must be square.');
    end
    augmentedMatrix = [A, b];
    n = numRows;

    % --- Main Gauss-Jordan Elimination Loop ---
    for i = 1:n
        % 1. Partial Pivoting
        % Find the row with the largest pivot element in the current column
        [~, maxRowIndex] = max(abs(augmentedMatrix(i:n, i)));
        maxRowIndex = maxRowIndex + i - 1; % Adjust index to be absolute

        % Swap the current row with the row containing the max pivot
        if maxRowIndex ~= i
            augmentedMatrix([i, maxRowIndex], :) = augmentedMatrix([maxRowIndex, i], :);
        end

        % Check for singularity
        if abs(augmentedMatrix(i, i)) < 1e-10 % Using tolerance for floating point comparison
            error('Matrix is singular or nearly singular. Cannot solve the system.');
        end

        % 2. Normalize the pivot row
        % Divide the pivot row by the pivot element to make the pivot element 1
        augmentedMatrix(i, :) = augmentedMatrix(i, :) / augmentedMatrix(i, i);

        % 3. Eliminate other entries in the current column
        for j = 1:n
            if i ~= j
                % Subtract a multiple of the pivot row from all other rows
                factor = augmentedMatrix(j, i);
                augmentedMatrix(j, :) = augmentedMatrix(j, :) - factor * augmentedMatrix(i, :);
            end
        end
    end

    % Extract the solution from the last column
    solution = augmentedMatrix(:, end);

    % --- Verification (Best Practice) ---
    residual = A * solution - b;
    if norm(residual) > 1e-6
        warning('Solution might be inaccurate. Residual norm: %e', norm(residual));
    end
end

% --- How to use the function ---
% A = [1 1 2; 2 4 -3; 3 6 -5];
% b = [9; 1; 0];
% mySolution = solveWithGaussJordan(A, b);

Code Explanation

• Function Definition: The code is wrapped in a function solveWithGaussJordan for reusability and clarity.
• Input Validation: It first checks if the matrix A is square, a requirement for this simple solver.
• Partial Pivoting: For each column i, it finds the row below it with the largest absolute value. This row is swapped with the current row i. This technique, called partial pivoting, is crucial for minimizing floating-point errors and avoiding division by zero.
• Singularity Check: After pivoting, it checks if the pivot element is close to zero. If so, the matrix is singular (or close to it), and a unique solution doesn’t exist. The program throws an error.
• Normalization: The pivot row is divided by the pivot element A(i,i) to make the new pivot element equal to 1.
• Elimination: For all other rows j, it subtracts a multiple of the (newly normalized) pivot row to make the element A(j,i) zero.
• Solution Extraction: After the loop completes, the coefficient part of the matrix is an identity matrix, and the last column contains the solution vector.
• Verification: A good practice is to check the solution. The residual (A*x - b) should be a vector of zeros. We check its norm (magnitude) to see if it’s acceptably small.

Conclusion & Best Practices

While implementing algorithms like Gauss-Jordan elimination is a great learning exercise, it’s vital to use the right tools for the job in a professional or academic setting.

• For Solving Ax = b: Always prefer the backslash operator: x = A \ b.
• For Finding RREF: Use the built-in function: R = rref(augmentedMatrix).
• Numerical Stability: MATLAB’s built-in methods are professionally developed and tested to handle a wide range of numerical challenges far better than a simple script.
• Performance: Built-in functions are implemented in lower-level languages (like C++ or Fortran) and are significantly faster than interpreted MATLAB code in a for loop.