SC-07

Course 7

Iterative Methods for Linear Systems

For large sparse systems, iterative methods are often more efficient than direct methods. Learn Jacobi, Gauss-Seidel, and SOR methods, and understand when they converge.

Learning Objectives

Understand matrix splitting and iterative formulations
Implement Jacobi, Gauss-Seidel, and SOR methods
Analyze convergence using spectral radius
Apply sufficient conditions for convergence
Choose optimal relaxation parameters
Understand block iterative methods

1. Jacobi Method

The Jacobi method solves $Ax = b$ by isolating each variable and iterating.

Definition 7.1: Matrix Splitting

Write $A = D - L - U$ where:

$D$ = diagonal of $A$
$-L$ = strict lower triangular part of $A$
$-U$ = strict upper triangular part of $A$

Definition 7.2: Jacobi Iteration

From $Dx = (L + U)x + b$ :

x^{(k+1)} = D^{-1}(L + U)x^{(k)} + D^{-1}b

Component-wise:

x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right)

Algorithm Jacobi Method

Input: $A, b, x^{(0)}$ , tolerance $\epsilon$

For $k = 0, 1, 2, \ldots$ until convergence:
For $i = 1, \ldots, n$ : compute $x_i^{(k+1)}$ using only $x^{(k)}$
If $\|x^{(k+1)} - x^{(k)}\| < \epsilon$ , stop

Remark:

Parallelizable: All components of $x^{(k+1)}$ can be computed independently since they only depend on $x^{(k)}$ .

2. Gauss-Seidel Method

Definition 7.3: Gauss-Seidel Iteration

Use updated values immediately as they become available:

x^{(k+1)} = (D - L)^{-1}Ux^{(k)} + (D - L)^{-1}b

Component-wise:

x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \right)

Example: Comparison

For system $4x_1 - x_2 = 1, -x_1 + 4x_2 = 1$ with $x^{(0)} = (0, 0)$ :

$k$	Jacobi $x^{(k)}$	Gauss-Seidel $x^{(k)}$
0	(0, 0)	(0, 0)
1	(0.25, 0.25)	(0.25, 0.3125)
2	(0.3125, 0.3125)	(0.328, 0.332)

Exact solution: $(1/3, 1/3)$ . Gauss-Seidel converges faster.

Remark:

Gauss-Seidel typically converges faster than Jacobi (roughly twice as fast for many problems), but it cannot be parallelized as easily due to data dependencies.

3. Successive Over-Relaxation (SOR)

Definition 7.4: SOR Iteration

Introduce relaxation parameter $\omega$ :

x_i^{(k+1)} = (1 - \omega)x_i^{(k)} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j^{(k+1)} - \sum_{j > i} a_{ij} x_j^{(k)} \right)

In matrix form: $x^{(k+1)} = (D - \omega L)^{-1}((1-\omega)D + \omega U)x^{(k)} + \omega(D - \omega L)^{-1}b$

Theorem 7.1: SOR Convergence

SOR can only converge if $0 < \omega < 2$ .

Note:

Special cases:

$\omega = 1$ : Gauss-Seidel
$\omega < 1$ : Under-relaxation (more stable, slower)
$\omega > 1$ : Over-relaxation (can accelerate convergence)

Theorem 7.2: Optimal SOR for Tridiagonal SPD

For symmetric positive definite tridiagonal matrices with Jacobi iteration matrix $B_J$ :

\omega_{\text{opt}} = \frac{2}{1 + \sqrt{1 - \rho(B_J)^2}}

This minimizes $\rho(B_{\text{SOR}})$ .

4. Convergence Analysis

Definition 7.5: Spectral Radius

The spectral radius of matrix $B$ is:

\rho(B) = \max_{i} |\lambda_i|

where $\lambda_i$ are the eigenvalues of $B$ .

Theorem 7.3: Convergence Criterion

The iteration $x^{(k+1)} = Bx^{(k)} + c$ converges to the unique solution for any initial guess $x^{(0)}$ if and only if:

\rho(B) < 1

Definition 7.6: Convergence Rate

The asymptotic convergence rate is:

R = -\log_{10}(\rho(B))

This measures decimal digits of accuracy gained per iteration.

Theorem 7.4: Sufficient Conditions

The following ensure convergence of both Jacobi and Gauss-Seidel:

Strict diagonal dominance: $|a_{ii}| > \sum_{j \neq i} |a_{ij}|$ for all $i$
Symmetric positive definite: $A = A^T$ and $x^T A x > 0$ for all $x \neq 0$

Example: Convergence Check

For matrix $A = \begin{pmatrix} 4 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 4 \end{pmatrix}$ :

Row 1: $|4| = 4 > |-1| + |0| = 1$ ✓
Row 2: $|4| = 4 > |-1| + |-1| = 2$ ✓
Row 3: $|4| = 4 > |0| + |-1| = 1$ ✓

Strictly diagonally dominant → Both methods converge.

5. Practical Considerations

Definition 7.7: Stopping Criteria

Common stopping criteria:

Residual: $\|b - Ax^{(k)}\| < \epsilon \|b\|$
Relative change: $\|x^{(k+1)} - x^{(k)}\| < \epsilon \|x^{(k+1)}\|$
Maximum iterations: $k < k_{\max}$

When to Use Iterative Methods

Large sparse matrices (thousands or millions of unknowns)
When matrix structure can be exploited (e.g., matrix-free methods)
When only approximate solutions are needed
When storage is limited

Block Methods

Partition unknowns into blocks and solve block systems at each iteration. Often converges faster than point methods.

Remark:

Modern iterative methods like Conjugate Gradient (for SPD) andGMRES (for general matrices) are often preferred over classical stationary methods. They are covered in advanced courses.

Practice Quiz

Iterative Linear Systems Quiz

Questions

Correct

Accuracy

In the Jacobi method, the update for

x_i^{(k+1)}

uses:

Easy

Not attempted

Gauss-Seidel differs from Jacobi by:

Easy

Not attempted

The SOR method introduces a parameter

\omega

that should satisfy:

Easy

Not attempted

An iterative method

x^{(k+1)} = Bx^{(k)} + c

converges for all initial guesses iff:

Medium

Not attempted

A matrix is strictly diagonally dominant if:

Medium

Not attempted

For a strictly diagonally dominant matrix, which methods converge?

Medium

Not attempted

The convergence rate of an iterative method is determined by:

Hard

Not attempted

For SPD matrices with SOR, the optimal

\omega

satisfies:

Hard

Not attempted

Block Jacobi differs from point Jacobi by:

Hard

Not attempted

A common stopping criterion for iterative methods is:

Medium

Not attempted

Frequently Asked Questions

When should I use iterative vs. direct methods?

Direct methods: For small to medium dense systems, or when you need to solve with many right-hand sides.

Iterative methods: For large sparse systems where direct methods would be too slow or use too much memory. Also useful when only approximate solutions are needed.

How do I choose between Jacobi, Gauss-Seidel, and SOR?

Jacobi: Simple to implement and parallelize. Often slowest to converge.

Gauss-Seidel: Usually faster than Jacobi (about 2x). Use when parallelization isn't critical.

SOR: Fastest when optimal $\omega$ is known. Worth the effort for repeated solves of similar systems.

My iteration isn't converging. What should I check?

Check these in order:

Is $A$ non-singular?
Is $A$ diagonally dominant or SPD?
For SOR, is $0 < \omega < 2$ ?
Try scaling the system (row equilibration)
Consider preconditioning or Krylov methods

What is preconditioning?

Preconditioning transforms $Ax = b$ into an equivalent system with better convergence properties. Instead of solving $Ax = b$ , solve $M^{-1}Ax = M^{-1}b$ where $M \approx A$ but $M^{-1}v$ is easy to compute. Good preconditioners dramatically speed up convergence.

How do I estimate the spectral radius without computing eigenvalues?

Run a few iterations and observe the ratio $\|x^{(k+1)} - x^{(k)}\| / \|x^{(k)} - x^{(k-1)}\|$ . This ratio approaches $\rho(B)$ asymptotically. Alternatively, use the power method on the iteration matrix (covered in the next chapter).