Burden & Faires, 9th Ed. - §7.3–7.5
Chapter 14: Iterative Methods for Linear Systems
Jacobi's Method · Gauss-Seidel · Successive Over-Relaxation (SOR) · Convergence criteria
This closing chapter solves $A\mathbf x=\mathbf b$ a completely different way than a linear-algebra course would - not by row-reducing $A$ exactly (Gaussian elimination, $O(n^3)$ and exact up to round-off), but by iterating toward the solution, the same way every root-finder in Chapters 2–4 iterated toward a root. For the very large, very sparse systems that come from discretizing differential equations, iterative methods can reach an acceptable answer in far fewer operations than a direct solve - provided they converge at all, which is exactly what this chapter's criteria pin down.
Section 7.3
Jacobi's Method
Derivation - Solve Each Equation for Its Own Diagonal Variable
Write $A=D+L+U$ (diagonal, strictly lower, strictly upper parts). Solving $A\mathbf x=\mathbf b$ formally for the diagonal part:
$$\mathbf x = D^{-1}(\mathbf b - (L+U)\mathbf x)$$
gives an obvious fixed-point iteration (exactly Chapter 3's $x=g(x)$ idea, now for vectors):
$$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j\neq i} a_{ij}x_j^{(k)}\right), \qquad i=1,\ldots,n$$
Every component of $\mathbf x^{(k+1)}$ uses only the previous iterate $\mathbf x^{(k)}$ - none of the newly computed components are used until the next full sweep. This makes Jacobi trivially parallelizable, at the cost of being the slower of the two classical methods.
📘 Example 14.1 - Jacobi's Method on a $3\times3$ System
Solve $\begin{pmatrix}10&1&1\\2&10&1\\2&1&10\end{pmatrix}\mathbf x = \begin{pmatrix}12\\13\\13\end{pmatrix}$ (true solution $\mathbf x=(1,1,1)$) by Jacobi's Method from $\mathbf x^{(0)}=\mathbf 0$.
Solution
$$x_1^{(k+1)}=\frac{12-x_2^{(k)}-x_3^{(k)}}{10}, \quad x_2^{(k+1)}=\frac{13-2x_1^{(k)}-x_3^{(k)}}{10}, \quad x_3^{(k+1)}=\frac{13-2x_1^{(k)}-x_2^{(k)}}{10}$$
Steady linear convergence toward $(1,1,1)$ - the spectral radius of Jacobi's iteration matrix here is $\rho(T_J)\approx0.2562$, so each sweep shrinks the error by roughly that factor.
| $k$ | $x_1$ | $x_2$ | $x_3$ | $\|\mathbf x^{(k)}-\mathbf x\|_\infty$ |
|---|---|---|---|---|
| 0 | 0 | 0 | 0 | 1.000 |
| 1 | 1.200 | 1.300 | 1.300 | 0.300 |
| 2 | 0.940 | 0.930 | 0.930 | 0.070 |
| 3 | 1.014 | 1.019 | 1.019 | 0.019 |
| 6 | 0.99975 | 0.99969 | 0.99969 | 0.00031 |
Section 7.3 continued
The Gauss-Seidel Method
Derivation - Use Updated Values Immediately
The obvious improvement: as soon as $x_1^{(k+1)}$ is computed, use it (instead of the stale $x_1^{(k)}$) when computing $x_2^{(k+1)}$, and so on down the sweep:
$$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j\lt i} a_{ij}x_j^{(k+1)} - \sum_{j\gt i} a_{ij}x_j^{(k)}\right)$$
In matrix form, $(D+L)\mathbf x^{(k+1)} = \mathbf b - U\mathbf x^{(k)}$. This single change - using the freshest information available rather than waiting a full sweep - typically (not always) roughly doubles the convergence rate over Jacobi, at the cost of being inherently sequential rather than parallelizable.
📘 Example 14.2 - Gauss-Seidel, Same System
Repeat Example 14.1 with Gauss-Seidel.
Solution
Gauss-Seidel reaches $10^{-7}$ accuracy by $k=5$ where Jacobi (Example 14.1) was still at $\sim10^{-3}$ by $k=6$ - its iteration-matrix spectral radius here is $\rho(T_G)\approx0.0447$, roughly the square of Jacobi's $0.2562$, which is a common (though not universal) pattern for this class of matrices.
| $k$ | $x_1$ | $x_2$ | $x_3$ | $\|\mathbf x^{(k)}-\mathbf x\|_\infty$ |
|---|---|---|---|---|
| 0 | 0 | 0 | 0 | 1.000 |
| 1 | 1.200 | 1.060 | 0.954 | 0.200 |
| 2 | 0.9986 | 1.0049 | 0.9998 | 0.0049 |
| 3 | 0.99953 | 1.00011 | 1.00008 | 0.00047 |
| 5 | 1.000000 | 0.999999 | 1.000000 | $4\times10^{-7}$ |
Section 7.3 continued
Successive Over-Relaxation (SOR)
Extrapolating Past the Gauss-Seidel Step
Let $x_i^{\text{GS}}$ be the ordinary Gauss-Seidel update. SOR takes a weighted step beyond (or short of) it:
$$x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \omega\, x_i^{\text{GS}}$$
$\omega=1$ recovers plain Gauss-Seidel exactly; $\omega>1$ ("over-relaxation") extrapolates further in the direction Gauss-Seidel was already moving, which can accelerate convergence substantially for the right matrix - or slow it down, or cause divergence, for the wrong one. A necessary condition for convergence at all is $0<\omega<2$.
Iterations needed to reach $\|\mathbf x^{(k)}-\mathbf x\|_\infty<10^{-6}$ as a function of $\omega$, for the Example 14.1/14.2 system. For this particular matrix, the minimum sits almost exactly at $\omega=1$ (plain Gauss-Seidel) - a reminder that over-relaxation is not automatically better; the optimal $\omega$ genuinely depends on the matrix's spectral properties, and for some systems (especially large ones from discretized PDEs) $\omega^*$ well above 1 gives dramatic speedups.
📘 Example 14.3 - SOR Sensitivity to $\omega$
Using the Example 14.1 system, compare iterations-to-convergence ($\|\mathbf x^{(k)}-\mathbf x\|_\infty<10^{-6}$) for several $\omega$.
Solution
The curve is sharply U-shaped around the optimum, and choosing $\omega$ poorly (here, too large) can cost far more than it saves - SOR is a genuinely useful tool, but only once $\omega^*$ has been estimated (analytically for special matrix classes, or empirically as here) rather than guessed.
| $\omega$ | Iterations to $10^{-6}$ |
|---|---|
| 0.5 (under-relaxed) | 21 |
| 0.9 | 8 |
| 1.0 (plain Gauss-Seidel) | 5 |
| 1.1 | 7 |
| 1.5 | 22 |
| 1.9 | did not converge in 30 iterations |
Section 7.3–7.5
Convergence Criteria
Theorem 14.1 - The General Convergence Criterion
Writing any of these methods as $\mathbf x^{(k+1)} = T\mathbf x^{(k)} + \mathbf c$ for the appropriate iteration matrix $T$ (Jacobi's $T_J=-D^{-1}(L+U)$, Gauss-Seidel's $T_G=-(D+L)^{-1}U$), the iteration converges to the unique solution for any starting vector $\mathbf x^{(0)}$ if and only if
$$\rho(T) < 1$$
- the exact vector/matrix generalization of Chapter 3's scalar condition $|g'(p)|<1$, using Chapter 13's spectral radius.
A Practical Sufficient Condition - Diagonal Dominance
Checking $\rho(T)<1$ directly requires computing eigenvalues - exactly the problem this chapter is trying to avoid needing a full solve for. A much cheaper sufficient (not necessary) condition: if $A$ is strictly diagonally dominant,
$$|a_{ii}| > \sum_{j\neq i}|a_{ij}| \quad \text{for every row } i$$
then both Jacobi's and Gauss-Seidel's methods are guaranteed to converge, for any starting vector - no eigenvalue computation needed. The system in Examples 14.1–14.3 is strictly diagonally dominant in every row ($10>1+1$, $10>2+1$, $10>2+1$), which is exactly why both methods converged so cleanly.
Diagonal Dominance Is Sufficient, Not Necessary
A system can fail strict diagonal dominance and still have $\rho(T)<1$ for one or both methods - the condition is a convenient, cheap-to-check guarantee, not a requirement. Conversely, a badly non-diagonally-dominant system offers no guarantee either way; convergence must be checked directly (or simply attempted, with a maximum-iteration safeguard as in every iterative method throughout this course).
Practice
Practice Problem
Try this problem yourself before expanding the solution below.
✏️ Practice Problem 14.1 - Gauss-Seidel on a New $3\times3$ System
Solve
$$7x-2y+z=3,\qquad -2x+8y-z=5,\qquad x-y+5z=4$$
using Gauss-Seidel starting from $x^{(0)}=(0,0,0)$, for three iterations. First confirm the coefficient matrix is strictly diagonally dominant.
Solution
Diagonal dominance: row 1: $|7|>|-2|+|1|=3$; row 2: $|8|>|-2|+|-1|=3$; row 3: $|5|>|1|+|-1|=2$. All three hold, so convergence is guaranteed from any starting vector.
Gauss-Seidel update, using the newest values as soon as they're available:
The exact solution is $(0.554217,\,0.871486,\,0.863454)$.
After only 3 iterations the estimate is within $0.004$ of the exact solution in every component, and a 4th iteration tightens that to within $0.0002$ - fast convergence, as expected from a diagonally dominant system.
Gauss-Seidel update, using the newest values as soon as they're available:
| Iter | $x$ | $y$ | $z$ |
|---|---|---|---|
| 1 | 0.428571 | 0.732143 | 0.860714 |
| 2 | 0.514796 | 0.861288 | 0.869298 |
| 3 | 0.550468 | 0.871279 | 0.864162 |
Connections Across the Course
- Every iterative method in this chapter is fixed-point iteration (Chapter 3) lifted from scalars to vectors - same convergence logic, same error-bound structure, now governed by spectral radius instead of $|g'(p)|$.
- $\rho(T)<1$ and the diagonal-dominance shortcut both draw directly on Chapter 13's matrix norms and eigenvalue theory.
- The tridiagonal systems from cubic splines (Chapter 7) are exactly the kind of large, sparse, diagonally dominant systems where iterative methods (or the direct $O(n)$ tridiagonal solve used there) genuinely outperform general-purpose Gaussian elimination.