Direct vs. iterative linear solvers
Gaussian elimination and LU/Cholesky's exact answer in a fixed number of steps vs. Jacobi/Gauss-Seidel/SOR/CG's cheap-per-sweep approximation.
For a small-to-medium DENSE system, or whenever you need an exact answer (up to round-off) with a predictable operation count, use a direct method -- Gaussian elimination, or its LU/Cholesky factorization if you'll solve against several right-hand sides or the matrix is symmetric positive definite. For a LARGE, SPARSE system (finite-difference and finite-element matrices are the classic case), an iterative method's per-sweep cost of touching only the nonzero entries wins decisively -- start with Jacobi or Gauss-Seidel if the matrix is diagonally dominant (both are then guaranteed to converge), SOR to accelerate Gauss-Seidel with a relaxation factor, or conjugate gradient specifically when the matrix is also symmetric positive definite, where it's unbeatable.
Direct methods -- exact (up to round-off), O(n^3)
Reorders rows so elimination's multipliers never exceed 1 in magnitude -- the round-off-control step plain elimination skips.
Elimination's own multipliers, saved as A=LU instead of thrown away -- reused across several right-hand sides for free.
Half the arithmetic of LU, for symmetric positive definite systems specifically.
Section
Why elimination costs O(n^3)
Iterative methods -- cheap per sweep, only approximate
The sufficient condition that guarantees Jacobi AND Gauss-Seidel both converge, from any starting vector.
If A is strictly diagonally dominant -- \lvert a_{ii}\rvert> \sum_{j\ne i}\lvert a_{ij}\rvert for every row i -- then both the Jacobi and Gauss-Seidel iterations converge to the unique solution of A\mathbf{x}=\mathbf{b}, from every starting vector \mathbf{x}^{(0)}.
Show proof
Direct, via a contraction argument in the infinity norm (Chapter 7, Section 6). Write the Jacobi iteration as \mathbf{x}^{(k+1)}= T_J\mathbf{x}^{(k)}+\mathbf{c}, where T_J=-D^{-1}(L+U) has entries (T_J)_{ij}=-a_{ij}/a_{ii} for j\ne i and (T_J)_{ii}=0. The exact solution \mathbf{x}^* satisfies the SAME fixed-point equation, \mathbf{x}^*=T_J\mathbf{x}^*+\mathbf{c}, since \mathbf{x}^* solving A\mathbf{x}=\mathbf{b} is exactly what makes every row's own rearrangement (Section 2) come out true. Subtracting the two equations gives \mathbf{x}^{(k+1)}-\mathbf{x}^*=T_J\bigl(\mathbf{x}^{(k)}- \mathbf{x}^*\bigr).
The infinity norm of T_J is its largest absolute row sum: \lVert T_J\rVert_\infty=\max_i\frac{\sum_{j\ne i}\lvert a_{ij}\rvert} {\lvert a_{ii}\rvert}. Strict diagonal dominance means every one of these ratios is strictly less than 1, so \lVert T_J\rVert_\infty=\delta for some \delta<1. Taking infinity norms of the subtracted equation and using that an induced matrix norm is submultiplicative (\lVert T_J\mathbf{v}\rVert_ \infty\le\lVert T_J\rVert_\infty\lVert\mathbf{v}\rVert_\infty, Chapter 7, Section 6): \bigl\lVert\mathbf{x}^{(k+1)}-\mathbf{x}^*\bigr\rVert_\infty\le\delta \bigl\lVert\mathbf{x}^{(k)}-\mathbf{x}^*\bigr\rVert_\infty. Applying this bound repeatedly, by induction on k, \bigl\lVert\mathbf{x}^{(k)}-\mathbf{x}^*\bigr\rVert_\infty\le\delta^k \bigl\lVert\mathbf{x}^{(0)}-\mathbf{x}^*\bigr\rVert_\infty. Since 0\le\delta<1, \delta^k\to0 as k\to\infty, regardless of how far \mathbf{x}^{(0)} started from \mathbf{x}^* -- so \mathbf{x}^{(k)}\to\mathbf{x}^* from ANY starting guess. This is exactly the "geometric-decay to zero" pattern Chapter 2's convergence-order vocabulary already named linear convergence, with rate at most \delta.
Gauss-Seidel's iteration matrix T_{GS}=-(D+L)^{-1}U satisfies the same bound, \lVert T_{GS}\rVert_\infty\le\delta, by an analogous but more delicate row-by-row argument: each row's bound on the already-updated entries (using \delta<1 inductively on earlier rows) combines with the same diagonal-dominance margin on the not-yet-updated entries to keep the overall ratio under 1 as well. The convergence conclusion for Gauss-Seidel then follows by the identical contraction argument above, applied to T_{GS} in place of T_J.
The general convergence test when diagonal dominance doesn't hold -- converges exactly when this is under 1.
Choosing
Let A be a nonsingular n\times n matrix, let \mathbf{x} solve A\mathbf{x}=\mathbf{b} exactly, and let \hat{\mathbf{x}} be any approximation to \mathbf{x} with residual \mathbf{r}=\mathbf{b}-A\hat{\mathbf{x}}. Then
\frac{\lVert\mathbf{x}-\hat{\mathbf{x}}\rVert}{\lVert\mathbf{x}\rVert}\le \kappa(A)\,\frac{\lVert\mathbf{r}\rVert}{\lVert\mathbf{b}\rVert},where \kappa(A)=\lVert A\rVert\,\lVert A^{-1}\rVert is A's condition number in the same matrix norm as \lVert\cdot\rVert.
A poorly conditioned system undermines a direct method's own exactness -- worth checking before trusting either approach's answer.
Show proof
Since A\mathbf{x}=\mathbf{b},
\mathbf{x}-\hat{\mathbf{x}}=A^{-1}\mathbf{b}-\hat{\mathbf{x}}=A^{-1}(\mathbf{b}-A\hat{\mathbf{x}})=A^{-1}\mathbf{r}.Taking norms and using the induced-norm inequality \lVert M\mathbf{v}\rVert\le\lVert M\rVert\,\lVert\mathbf{v}\rVert (true of every matrix norm induced by a vector norm, by the very definition of \lVert M\rVert as the largest such ratio):
\lVert\mathbf{x}-\hat{\mathbf{x}}\rVert=\lVert A^{-1}\mathbf{r}\rVert\le\lVert A^{-1}\rVert\,\lVert\mathbf{r}\rVert.The same inequality applied to \mathbf{b}=A\mathbf{x} gives \lVert\mathbf{b}\rVert\le\lVert A\rVert\,\lVert\mathbf{x}\rVert, so
\lVert\mathbf{x}\rVert\ge\frac{\lVert\mathbf{b}\rVert}{\lVert A\rVert}.Dividing the first bound by this lower bound on \lVert\mathbf{x}\rVert:
\frac{\lVert\mathbf{x}-\hat{\mathbf{x}}\rVert}{\lVert\mathbf{x}\rVert} \le\frac{\lVert A^{-1}\rVert\,\lVert\mathbf{r}\rVert}{\lVert\mathbf{b}\rVert/\lVert A\rVert} =\lVert A\rVert\,\lVert A^{-1}\rVert\,\frac{\lVert\mathbf{r}\rVert}{\lVert\mathbf{b}\rVert} =\kappa(A)\,\frac{\lVert\mathbf{r}\rVert}{\lVert\mathbf{b}\rVert}.Sanity check: if \hat{\mathbf{x}}=\mathbf{x} exactly, then \mathbf{r}=\mathbf{0} and both sides are 0 — the bound is tight in the one case where there's nothing left to bound.