Gauss-Seidel Iteration

Speed up Jacobi by using each updated value the moment it is computed within the same sweep.

The one-line change

Gauss-Seidel is Jacobi with immediate reuse: as you sweep through the unknowns, each x[i] is computed using the already-updated x[0..i-1] from this sweep and the not-yet-updated x[i+1..] from the last. New information propagates within a sweep instead of waiting for the next one.

import numpy as np
A=np.array([[4.0,1.0,0.0],[1.0,5.0,1.0],[0.0,1.0,3.0]])
b=np.array([5.0,7.0,4.0]); x=np.zeros(3); n=3
for it in range(100):
    xold=x.copy()
    for i in range(n):
        s=sum(A[i,j]*x[j] for j in range(n) if j!=i)
        x[i]=(b[i]-s)/A[i,i]           # uses updated x in place
    if np.max(np.abs(x-xold))<1e-10: break
print(np.round(x,5),'in',it,'iters')  # fewer iters than Jacobi

Convergence

Gauss-Seidel converges for any symmetric positive-definite matrix and for strictly diagonally dominant matrices, and it usually takes about half as many iterations as Jacobi. The gain comes from feeding fresh values forward, effectively a better approximation to the inverse at each sweep.

The trade-off

Because each update depends on the ones just computed, the basic method is inherently sequential - harder to parallelize than Jacobi. Red-black ordering recovers parallelism by splitting the grid into two independent color sets that each update in a Jacobi-like way, a common trick in PDE solvers.

Toward multigrid

Gauss-Seidel is an excellent smoother: it kills high-frequency error quickly but low-frequency error slowly. Multigrid methods exploit exactly this by moving the slow, smooth error to coarser grids where it looks high-frequency, achieving near-optimal convergence for elliptic problems such as Grad-Shafranov.