Solving Linear Systems with Gaussian Elimination - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import matplotlib.pyplot as plt
from sympy import Matrix, symbols, pprint

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

The Intuition - solution by elimination.¶

Gaussian elimination is a formal algorithm for the most familiar and intuitive algebraic solution technique: substitution and elimination.

Let’s start with a classic problem:

You are organizing a fundraising event and need to buy chairs and tables. Chairs cost $20 each and tables cost $50 each. You have a budget of $700 and need a total of 20 pieces of furniture. How many chairs and tables should you buy?

Symbolic Manipulation¶

Let $c$ and $t$ be the number of chairs and tables respectively. We can define our two equations:

(1) $20c + 50t = 700$
(2) $c + t = 20$

Solve (2) for $c$ :
$c = 20 - t$

Substitute into (1):
$20(20 - t) + 50t = 700 \implies t = 10$

Substitute $t=10$ back into (2) to find $c=10$ . This works easily because our first step carries the unknown symbol $t$ .

The Numerical Approach¶

Computers don’t like carrying symbols. Let’s solve it numerically using only the coefficients.

(1) $20c + 50t = 700$
(2) $c + t = 20$

Multiply (2) by 20:
$20c + 20t = 400$

Subtract from (1):

\begin{aligned} 20c + 50t &= 700 \\ -(20c + 20t) &= -400 \\ 30t &= 300 \end{aligned}

(1)

We eliminated $c$ from the second equation! It is now trivial to solve for $t$ , and then use back-substitution to find $c$ .

Transformation to Matrix Form¶

Notice how the linear system was fundamentally changed (without altering the final answer). In matrix form, we transformed it into:

\begin{pmatrix} 20 & 50 \\ 0 & 30 \end{pmatrix} \begin{pmatrix} c \\ t \end{pmatrix} = \begin{pmatrix} 700 \\ 300 \end{pmatrix}

(2)

The coefficient matrix $\mathbf{A}$ is now upper triangular! This is the core goal of Gaussian Elimination.

Upper Triangular Matrices¶

An upper triangular matrix $\mathbf{U}$ is a matrix whose elements below the main diagonal are exactly zero:

\begin{bmatrix} u_{1,1} & u_{1,2} & u_{1,3} & u_{1,4}\\ 0 & u_{2,2}' & u_{2,3}' & u_{2,4}'\\ 0 & 0 & u_{3,3}' & u_{3,4}' \\ 0 & 0 & 0 & u_{4,4}' \end{bmatrix}

(3)

This is incredibly useful because the system $\mathbf{U}\mathbf{x} = \mathbf{b}'$ can be solved from the bottom-up in $\mathcal{O}(n^2)$ time via back-substitution.

The Gaussian Elimination Algorithm¶

Gaussian Elimination transforms $\mathbf{A}\mathbf{x} = \mathbf{b}$ into $\mathbf{U}\mathbf{x} = \mathbf{b}'$ .

We build an augmented matrix $\mathbf{Ab} = [ \mathbf{A} | \mathbf{b} ]$ . To eliminate variables, we make a sequence of passes row-by-row, choosing a pivot row to eliminate elements in the equations below it.

We are allowed the following fundamental row operations:

Operation	Effect on solution	Effect on $\lvert\mathbf{A}\rvert$
Exchange 2 rows	None	Flips sign
Multiply a row by scalar $a$	None	Multiplied by $a$
Subtract 2 rows	None	Unchanged

(Note: While the solution $\mathbf{x}$ remains identical, the determinant, norm, and condition number of the matrix are actively changing during elimination!)

Example: Guassian elimination¶

Let’s use the algorithm to solve:

\begin{aligned} 4x_1 + 3x_2 - 5x_3 &= 2 \\ -2x_1 - 4x_2 + 5x_3 &= 5 \\ 8x_1 + 8x_2 &= -3 \end{aligned}

(4)

In the following, note the definition of coefficients $m_{i,j}$ - they will be reference when we talk about LU decomposition.

Step 1: Get the augmented matrix $[\mathbf{A} | \mathbf{b}]$

[\mathbf{A} | \mathbf{b}] = \begin{bmatrix} 4 & 3 & -5 & 2\\ -2 & -4 & 5 & 5\\ 8 & 8 & 0 & -3\\ \end{bmatrix}

(5)

First Elimination Pass¶

Step 2: Eliminate $x_1$ from Row 2 Choose Row 1 as the pivot. Multiply Row 1 by $m_{2,1} = -0.5$ and subtract it from Row 2:

\begin{bmatrix} 4 & 3 & -5 & 2\\ 0 & -2.5 & 2.5 & 6\\ 8 & 8 & 0 & -3\\ \end{bmatrix}

(6)

Step 3: Eliminate $x_1$ from Row 3 Multiply Row 1 by $m_{3,1} = 2$ and subtract it from Row 3:

\begin{bmatrix} 4 & 3 & -5 & 2\\ 0 & -2.5 & 2.5 & 6\\ 0 & 2 & 10 & -7\\ \end{bmatrix}

(7)

Second Pass & Substitution¶

Step 4: Eliminate $x_2$ from Row 3 Now Row 2 is the pivot. Multiply Row 2 by $m_{3,2} = -0.8$ and subtract it from Row 3: $$

\begin{bmatrix} 4 & 3 & -5 & 2\\ 0 & -2.5 & 2.5 & 6\\ 0 & 0 & 12 & -2.2\\ \end{bmatrix}

(8)

$$

$\mathbf{A}$ is now upper triangular!

Step 5: Back-substitution

Bottom row: $12x_3 = -2.2 \implies x_3 = -0.183$
Middle row: $-2.5x_2 + 2.5(-0.183) = 6 \implies x_2 = -2.583$
Top row: $4x_1 + 3(-2.583) - 5(-0.183) = 2 \implies x_1 = 2.208$

Computational Complexity¶

If we analyze the two phases of standard Gaussian Elimination separately, we find:

Elimination Phase (Forward): Requires $\approx n^3/3$ operations, meaning it scales as $\mathcal{O}(n^3)$ .
Back Substitution: Requires $\approx n^2/2$ operations, meaning it scales as $\mathcal{O}(n^2)$ .

For massive systems (e.g., $n = 10,000$ ), the elimination phase completely dominates the computational time.

Gauss-Jordan Elimination¶

An obvious extension is to conduct the elimination passes both downwards and upwards, normalizing the pivot rows to 1. This creates the reduced row echelon form (rref):

\begin{bmatrix} 1 & 0 & 0 & 2.208\\ 0 & 1 & 0 & -2.583\\ 0 & 0 & 1 & -0.183 \end{bmatrix}

(9)

Since the coefficient matrix is now the Identity matrix $\mathbf{I}$ , the answer is simply the right-hand vector!

Is Gauss-Jordan Faster? No.¶

Because Gauss-Jordan eliminates the need for back-substitution, one might think it is more efficient. Unfortunately, it is not:

The full up-and-down elimination phase requires $\approx n^3/2$ operations.
Standard Gaussian Elimination requires $\approx n^3/3 + n^2/2$ operations.

Therefore, standard Gaussian Elimination is mathematically preferred for solving large systems.

Package Implementations¶

Because Gaussian Elimination has largely been surpassed by matrix decompositions (like LU decomposition), you will rarely use raw GE in numerical packages like numpy.

However, it is still incredibly useful for symbolic manipulation, and you will find it natively in sympy under rref():

# Define the coefficient matrix and the right-hand side vector
A = Matrix([[1, 2, 3],
            [4, 5, 6],
            [7, 8, 9]])
b = Matrix([2, 5, -3])

# Combine into an augmented matrix
Ab = A.row_join(b)

# Calculate the reduced row echelon form (rref)
Ab_rref = Ab.rref()

print("Augmented matrix in rref:")
pprint(Ab_rref[0])

# Calculate the echelon form of the coefficient matrix
Ab_echelon = Ab.echelon_form()

print("\nCoefficient matrix in echelon form:")
pprint(Ab_echelon)