22  Jacobi and Gauss-Seidel Methods

Standard stationary iterations based on diagonal and triangular splittings.

22.1 Jacobi Iteration

Decomposes \(\bA\) into diagonal \(\bD\) and remainder \(\bR\): \(\bA = \bD + \bR\).

NoteDefinition: Jacobi Scheme

\[ \begin{align} \bx^{(k+1)} = \bD^{-1}(\bb - \bR\bx^{(k)}) \end{align} \]

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

  • Concurrency: All components update independently using values from the previous step.

Tip

Performance: Jacobi iteration is inherently parallel, as all components can be updated simultaneously. However, it often exhibits slow convergence and requires double buffering of the solution vector.

22.2 Gauss-Seidel Iteration

Splits \(\bA = \bL_* + \bU_*\), where \(\bL_*\) includes the diagonal.

NoteDefinition: Gauss-Seidel Scheme

\[ \begin{align} \bL_* \bx^{(k+1)} = \bb - \bU_* \bx^{(k)} \end{align} \]

  • Component form: \(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)\).

  • Key Idea: Use updated values immediately within the same sweep.

Tip

Performance: Sequential dependencies restrict efficient parallelization, but Gauss-Seidel typically converges twice as fast as Jacobi and allows for in-place updates.

22.3 Convergence Conditions

NoteTheorem: Diagonal Dominance

If \(\bA\) is strictly diagonally dominant (\(|a_{ii}| > \sum_{j \neq i} |a_{ij}|\)), both Jacobi and Gauss-Seidel converge for any \(\bx^{(0)}\).

NoteTheorem: SPD Convergence

If \(\bA\) is SPD, Gauss-Seidel is guaranteed to converge. Jacobi may still diverge.

WarningExercise

  1. For \(\bA = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\), compute one step of Jacobi and Gauss-Seidel from \(\bzero\).

  2. Compare \(\rho(\bT_J)\) and \(\rho(\bT_{GS})\). Does \(\rho_{GS} \approx \rho_J^2\) hold?

22.4 Successive Over-Relaxation (SOR)

Accelerates Gauss-Seidel by overshooting the update.

NoteDefinition: SOR Update

\[ \begin{align} x_i^{(k+1)} = (1-\omega)x_i^{(k)} + \omega \left[\text{GS update}\right] \end{align} \]

  • Convergence Condition: \(0 < \omega < 2\) (Kahan’s Theorem).
Tip

Convergence Acceleration: For the 1D Poisson problem on a grid with spacing \(h\):

  • Jacobi and Gauss-Seidel require \(O(n^2)\) iterations.

  • Optimal SOR requires only \(O(n)\) iterations.

SOR effectively improves the convergence rate from \(O(h^2)\) to \(O(h)\).

WarningExercise

  1. Optimal Tuning: For \(\rho_J = 0.95\), calculate the optimal \(\omega^* = 2/(1 + \sqrt{1-\rho_J^2})\).

  2. Convergence Sweep: Implement SOR and plot iterations to convergence vs. \(\omega \in [1, 2)\). Identify the ``cliff’’ at \(\omega^*\).

22.5 Stopping Criteria

  1. Absolute Residual: \(\|\br^{(k)}\| < \varepsilon\).

  2. Relative Residual: \(\|\br^{(k)}\|/\|\bb\| < \varepsilon\). (Scale invariant).

  3. Increment: \(\|\bx^{(k+1)} - \bx^{(k)}\| < \varepsilon\). (Cheap, but may stop prematurely if convergence is slow).

WarningExercise

  1. Why is a small residual not enough to guarantee a small error \(\be^{(k)}\)? (Condition number).

  2. Implement a robust solver that combines relative residual with a max\_iter budget.