Numerical Linear Systems

Solving Ax = b efficiently and accurately. Learn Gaussian elimination cost and pivoting, LU decomposition for multiple right-hand sides, iterative methods (Jacobi, Gauss–Seidel, SOR), convergence conditions (diagonal dominance, spectral radius), conjugate gradient for SPD systems, sparsity and the condition number.

FEM and CFD models produce systems with millions of unknowns but mostly-zero (sparse) matrices. Direct elimination becomes impractical; iterative methods and preconditioning make them solvable.

An analogy

Iterative solvers are like tuning a guitar string by string: each adjustment uses the latest tuning of the others, and the whole instrument gradually converges.

Direct methods

Gaussian elimination with partial pivoting (swap to the largest pivot) is robust. LU factorisation reuses the factorisation for many right-hand sides (e.g. many load cases).

Iterative methods

Jacobi updates all unknowns from the previous sweep; Gauss–Seidel uses new values immediately (usually twice as fast); SOR over-relaxes with 1 < ω < 2 to accelerate further.

Convergence

Guaranteed for strictly diagonally dominant or SPD matrices (GS). Rate depends on the spectral radius of the iteration matrix; convergence stalls for poorly conditioned problems.

Krylov methods

Conjugate gradient (SPD), GMRES and BiCGSTAB (general) with preconditioners (ILU, multigrid) are the workhorses of large-scale simulation. Sparse storage formats (CSR) and reordering reduce memory and fill-in.

Formulas

Worked example: Diagonal dominance check

Given: A = [[4, −1, 0], [−1, 4, −1], [0, −1, 4]]

Find: Will Jacobi/GS converge?

Answer: Yes — strictly diagonally dominant

Worked example: Cost comparison

Given: n = 10⁵ unknowns, ≈ 5 non-zeros per row

Find: Dense elimination vs one iterative sweep.

Answer: Iterative sweeps are ~10⁸× cheaper each

Gauss–Seidel iteration (Python/NumPy)

import numpy as np

# Diagonally dominant system (e.g. 1D heat / resistor ladder)
A = np.array([[4., -1, 0], [-1, 4, -1], [0, -1, 4]])
b = np.array([15., 10, 10])

x = np.zeros(3)
for k in range(1, 100):
    x_old = x.copy()
    for i in range(3):                       # Gauss-Seidel sweep
        s = b[i] - A[i, :i] @ x[:i] - A[i, i+1:] @ x_old[i+1:]
        x[i] = s / A[i, i]
    if np.max(np.abs(x - x_old)) < 1e-10:
        break
print("Gauss-Seidel:", np.round(x, 8), f"after {k} sweeps")
print("Direct solve:", np.round(np.linalg.solve(A, b), 8))
# Gauss-Seidel: [4.91071429 4.64285714 3.66071429] after 14 sweeps
# Direct solve: [4.91071429 4.64285714 3.66071429]

Diagonal dominance guarantees convergence. Each sweep costs O(nonzeros), so for sparse FEM matrices many cheap sweeps beat one expensive factorisation.

Common mistakes

Applications

Key takeaways

Going further