Solving Linear Systems: Augmented Matrices, Gaussian Elimination, and Gauss-Jordan

Sun Sep 13 2026

Part 4 of 5 in Matrix Algebra: Systems, Pivots, and Inverses

Solving a system of linear equations is one of the most fundamental tasks in applied mathematics. A 3×3 system like

2x + y − z = 8 −3x − y + 2z = −11 −2x + y + 2z = −3

can be solved by inspection or substitution for small systems, but when the number of variables grows, row reduction on an augmented matrix becomes the systematic approach taught in every linear algebra course.

In this post we write a Python implementation from scratch (no NumPy) that walks through each step: building the augmented matrix [A|b], demonstrating elementary row operations, then solving the same system twice — once via Gaussian elimination (forward elimination + back substitution) and once via Gauss-Jordan elimination (full reduction to reduced row echelon form). We also show how the method detects an inconsistent system.

The Three Elementary Row Operations

There are exactly three elementary row operations, each of which preserves the solution set:

1. Row scaling — multiply a row by a nonzero scalar
2. Row swapping — exchange two rows
3. Row replacement — add a multiple of one row to another

These operations are the only moves allowed during row reduction. The key insight is that every operation keeps the system’s solutions intact while reorganizing it into a form where the unknowns become easier to isolate.

Here’s how all three look applied to a 3×3 augmented matrix:

These operations are the atomic building blocks of every algorithm below.

The Code

The full script is in linear_systems.py. It’s structured as a single file with three demo functions — one for each example — plus helper routines that do the actual matrix manipulation. Here’s the complete source:

"""
Solve linear systems using augmented matrices, elementary row operations,
Gaussian elimination, and Gauss-Jordan elimination.
"""


def _print_matrix(M):
    """Pretty-print an augmented matrix with a vertical separator."""
    n_rows = len(M)
    n_cols = len(M[0]) if M else 0
    for i, row in enumerate(M):
        parts = []
        for j, val in enumerate(row):
            if isinstance(val, float) and val == int(val) and abs(val) < 1e10:
                formatted = str(int(val))
            elif isinstance(val, (int, float)) and abs(val) > 1e-15:
                formatted = f"{val:.2f}"
            else:
                formatted = "0.00"
            parts.append(formatted.rjust(8))
        sep_before_last = "|" if n_cols > 1 else ""
        line = f"  R{i+1}: {parts[0]}"
        for k in range(1, len(parts) - 1):
            line += f"  {parts[k]}"
        line += f"{sep_before_last} {parts[-1]}"
        print(line)


def elementary_row_operations_demo(equations):
    """Demonstrate the three types of elementary row operations."""
    M = [row[:] for row in equations]
    print(f"Original augmented matrix:")
    _print_matrix(M)

    # Operation 1: Row scaling
    print("Operation 1 — Row Scaling:")
    print("  R2 ← 3 * R2")
    M[1] = [x * 3 for x in M[1]]
    _print_matrix(M)

    # Operation 2: Row swapping
    print("Operation 2 — Row Swapping:")
    print("  Swap R1 <-> R3")
    M[0], M[2] = M[2][:], M[0][:]
    _print_matrix(M)

    # Operation 3: Row replacement
    print("Operation 3 — Row Replacement:")
    print("  R2 ← R2 + (-1) * R1")
    M[1] = [M[1][j] + (-1) * M[0][j] for j in range(len(M[1]))]
    _print_matrix(M)


def gaussian_elimination(A, b):
    """Solve Ax = b via forward elimination to row echelon form, then back substitution."""
    n = len(A)
    M = []
    for i in range(n):
        M.append([A[i][j] for j in range(n)])
        M[i].append(b[i])

    print(f"Augmented matrix [A|b]:")
    _print_matrix(M)

    # Forward elimination
    for col in range(n):
        max_row = col
        max_val = abs(M[col][col])
        for row in range(col + 1, n):
            if abs(M[row][col]) > max_val:
                max_val = abs(M[row][col])
                max_row = row

        if max_row != col:
            M[col], M[max_row] = M[max_row][:], M[col][:]
            print(f"  Pivot: Swap R{col+1} <-> R{max_row+1}")
            _print_matrix(M)

        pivot = M[col][col]
        if abs(pivot) < 1e-12:
            continue

        for row in range(col + 1, n):
            factor = M[row][col] / pivot
            M[row] = [M[row][j] - factor * M[col][j] for j in range(n + 1)]
            print(f"  Eliminate col {col+1}: R{row+1} <- ...")
            _print_matrix(M)

    # Back substitution
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        if abs(M[i][i]) < 1e-12:
            raise ValueError(f"Zero diagonal at row {i+1}.")
        x[i] = M[i][n]
        for j in range(i + 1, n):
            x[i] -= M[i][j] * x[j]
        x[i] /= M[i][i]
    return x


def gauss_jordan_elimination(A, b):
    """Solve Ax = b: reduce to reduced row echelon form [I|x]."""
    n = len(A)
    M = []
    for i in range(n):
        M.append([A[i][j] for j in range(n)])
        M[i].append(b[i])

    print(f"Augmented matrix [A|b]:")
    _print_matrix(M)

    for col in range(n):
        max_row = col
        max_val = abs(M[col][col])
        for row in range(col + 1, n):
            if abs(M[row][col]) > max_val:
                max_val = abs(M[row][col])
                max_row = row

        if max_row != col:
            M[col], M[max_row] = M[max_row][:], M[col][:]
            print(f"  Pivot: Swap R{col+1} <-> R{max_row+1}")
            _print_matrix(M)

        pivot = M[col][col]
        scale = 1.0 / pivot
        M[col] = [x * scale for x in M[col]]
        print(f"  Scale R{col+1} by {scale:.4f}:")
        _print_matrix(M)

        for row in range(n):
            if row == col:
                continue
            factor = M[row][col]
            if abs(factor) > 1e-12:
                M[row] = [M[row][j] - factor * M[col][j] for j in range(n + 1)]
                print(f"  Eliminate R{row+1}: ...")
                _print_matrix(M)

    x = [M[i][n] for i in range(n)]
    return x

Running It

Example 1: A 3×3 system

The first example solves the classic 3×3 system:

2x + y − z = 8 −3x − y + 2z = −11 −2x + y + 2z = −3

Both methods converge to the same answer: x = 2, y = 3, z = −1, and both pass verification. The Gaussian elimination path uses fewer steps because it only eliminates below the pivot, then back-substitutes. Gauss-Jordan does extra work by eliminating above each pivot too — which is why you can read the solution directly from the final column instead of doing a separate back substitution pass.

Example 2: A 4×4 system with pivoting

The second example is a 4×4 system where the top-left element is zero. Gaussian elimination would silently break without partial pivoting — swapping rows to bring the largest absolute value into the pivot position. The script handles this automatically:

1. At col 0: swaps R1 ↔ R3 (bringing 2.00 into the pivot spot)
2. At col 2: swaps R2 ↔ R4 (the row with the largest entry below the pivot)

Pivoting isn’t just a safety measure — it reduces numerical error by avoiding divisions by tiny numbers.

Example 3: Detecting an inconsistent system

The final example demonstrates what happens when no solution exists:

x + y = 2 x + y = 3

After eliminating col 0, row 2 reads [0.00, 0.00 | 1.00] — meaning 0 · x + 0 · y = 1, which is impossible. Any implementation that encounters a zero pivot with a nonzero right-hand side can immediately declare the system inconsistent.

Takeaway

Gaussian elimination and Gauss-Jordan reduction are the same algorithm viewed from two angles: Gaussian stops at row echelon form and back-substitutes (fewer operations), while Gauss-Jordan goes all the way to reduced row echelon form (the answer is right there in the last column). Both require three elementary row operations, both benefit from partial pivoting for numerical stability, and both naturally reveal when a system has no solution or infinitely many. The augmented matrix [A|b] is the data structure that makes this work — it carries both the coefficients and constants through every operation, so nothing gets lost.