Numerical Methods · Learn
How to solve a big linear system without ever inverting the matrix: guess, sweep through the equations, repeat. Why it works, when it's guaranteed to work, and why swapping two rows can turn divergence into convergence.
Solve $A\mathbf{x}=\mathbf{b}$ for $n$ unknowns. Gaussian elimination does it directly in about $n^3/3$ operations. That's fine for $n=10$ and hopeless for $n=40{,}000$, the size of a modest 200 × 200 finite-difference grid. Elimination also fills in zeros as it goes, so memory explodes too.
But engineering matrices are sparse. A node on a heat-conduction grid talks only to its four neighbours, so each row has five nonzeros out of 40,000. Iterative methods exploit that: each sweep touches only the nonzeros, costing about as much as one matrix–vector product, and you stop as soon as the answer is accurate enough.
Solve equation $i$ for unknown $i$, using the previous sweep's values for everything else:
Split the matrix into its diagonal, strictly lower and strictly upper parts, $A = D + L + U$. Then a whole Jacobi sweep is
Every unknown updates from the same old vector, so the updates are independent of one another. That makes Jacobi trivially parallel, which is why it's still used on GPUs.
Same formula, one change: as soon as $x_1^{(k+1)}$ is computed, use it when updating $x_2$, and so on down the list.
It needs only one copy of $\mathbf{x}$ in memory, and it usually converges faster because each update already reflects part of the current sweep.
Gauss–Seidel's corrections tend to point the right way but fall short. Successive over-relaxation stretches each one by a factor $\omega$:
$\omega = 1$ is plain Gauss–Seidel, $\omega>1$ over-relaxes, and $\omega<1$ under-relaxes, which helps tame a stubborn iteration. SOR can only converge for $0<\omega<2$. For grid problems the best $\omega$ sits close to 2, and getting it right can cut the number of sweeps from thousands to tens.
Each method has the form $\mathbf{x}^{(k+1)} = G\,\mathbf{x}^{(k)} + \mathbf{c}$, so the error obeys $\mathbf{e}^{(k+1)} = G\,\mathbf{e}^{(k)}$. The iteration matrices are
The iteration converges from every starting guess exactly when the spectral radius, the largest eigenvalue magnitude, is below one:
That's the true criterion, but eigenvalues are hard to see by eye. These conditions are easy to check, and each one is sufficient:
| If A is… | then… | typical source |
|---|---|---|
| strictly diagonally dominant, $|a_{ii}| > \sum_{j\ne i}|a_{ij}|$ in every row | Jacobi and Gauss–Seidel both converge | heat balances, circuits |
| weakly dominant ($\ge$) in every row, strict in at least one, and irreducible (no part of the system is cut off from the rest) | Jacobi and Gauss–Seidel both converge | finite-difference grids with a fixed-temperature boundary |
| symmetric positive definite | Gauss–Seidel converges, and SOR too for any $0<\omega<2$ | energy methods, structures |
The middle row matters more than it looks. An interior node of the 5-point Laplacian has $|a_{ii}| = 4$ and four neighbours of $-1$. That's dominance with equality, not strict. Weak dominance alone isn't enough: a fully insulated plate is singular and has no unique answer. The fixed-temperature edges supply the strict rows that make the guarantee work.
Dominance depends on which equation sits in which row. Shuffle the rows and a dominant system stops being dominant, and Gauss–Seidel can blow up. Swapping them back fixes it without changing the solution. Play Make It Converge →
Each sweep multiplies the error by roughly $\rho(G)$, so reaching a tolerance $\varepsilon$ takes about
When $\rho$ is close to 1, as it is on fine grids, $\ln\rho$ is tiny and the sweep count balloons. For a large class of matrices, including tridiagonal systems and the standard finite-difference Laplacian, $\rho(G_{\text{GS}}) = \rho(G_{\text{J}})^2$. Gauss–Seidel then needs half as many sweeps as Jacobi, and optimal SOR does dramatically better than both.
| Method | Uses | Memory | Speed | Watch out for |
|---|---|---|---|---|
| Jacobi | last sweep only | two vectors | slowest | slow when $\rho$ is near 1 |
| Gauss–Seidel | newest values | one vector | often ~2× Jacobi | sequential, depends on equation order |
| SOR | GS plus overshoot $\omega$ | one vector | much faster at a good $\omega$ | a bad $\omega$ is slower, or diverges outside $(0,2)$ |
Where these systems come from: every sweep of the 2D conduction lab is SOR on a finite-difference grid. The irregular-conduit lab builds its matrix automatically from the mesh, and the Ritz–Galerkin method produces small dense systems that are simplest to solve directly.
EngineeringCandy · Learn · the theory behind the playground