Jacobi, SOR and conjugate gradients on the Poisson equation

The problem is $-{\nabla }^{2}u=f$ on the unit square with $u=0$ on the edge, on an $n$ by $n$ grid of interior points with spacing $h=1/(n+1)$. The five-point stencil turns it into $n^2$ linear equations,

one per grid point, and $Au=b$ for short. Direct elimination is possible at this size, but the iterative methods here need only the stencil, so they are what gets used when the grid is large or three-dimensional.

The forcing $f=2[x(1-x)+y(1-y)]$ has the exact solution $u=x(1-x)y(1-y)$. The second difference of a quadratic is exact, so the grid solution equals $u$ at every grid point, and whatever error remains is the solver's.

Every method starts from $u=0$ and stops when $\Vert b-Au{\Vert }_2/\Vert b{\Vert }_2<10^{-8}$. The Python twin is posts/refs/jacobi-sor-and-cg-on-the-poisson-equation.py (output beside it).

Jacobi

A Jacobi sweep solves each equation for its own unknown, using the neighbors' values from the previous sweep: $u_{ij}\leftarrow \frac{1}{4}(u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}+b_{ij})$. All points update at once, so it is a few whole-array operations.

Its error obeys $e\leftarrow (I-A/4)e$. The eigenvectors of $A$ are the grid sine modes $\sin (j\pi x)\sin (k\pi y)$, with eigenvalues $4-2\cos j\pi h-2\cos k\pi h$, so a sweep multiplies each mode by $\frac{1}{2}(\cos j\pi h+\cos k\pi h)$. The smoothest mode, $j=k=1$, is multiplied by $\cos \pi h$, the largest factor, and after enough sweeps it is the only one left: the residual then shrinks by exactly $\cos \pi h$ per sweep, $0.98079$ at $h=1/16$. Getting from 1 to $10^{-8}$ at that rate takes $\ln 10^{-8}/\ln \cos \pi h=949$ sweeps.

Gauss-Seidel and SOR

Gauss-Seidel does the same update point by point, in order, and uses each new value as soon as it has it. For this matrix its error factor per sweep is exactly the square of Jacobi's, ${\cos }^2\pi h$, so it takes half the sweeps.

SOR (successive over-relaxation) takes the Gauss-Seidel update and goes further in the same direction: $u_{ij} \leftarrow u_{ij} + \omega\,(u^{GS}_{ij}

• u_{ij})$, with $1<\omega <2$. David Young worked out the best factor for

this matrix in 1950: $\omega =2/(1+\sin \pi h)$, which makes the error factor $\omega -1=(1-\sin \pi h)/(1+\sin \pi h)\approx 1-2\pi h$. Jacobi's is $\cos \pi h\approx 1-{\pi }^2{h}^2/2$. The distance from 1 is what sets the number of sweeps, so Jacobi's count grows like $1/h^2$ and SOR's like $1/h$. At $h=1/16$, $\omega =1.6735$.

At the optimal $\omega$ the eigenvalues of the SOR iteration matrix pair up and it cannot be diagonalized, so the residual is not a constant times $(\omega -1{)}^k$: the multiple grows, to 111 at sweep 30 and 147 at sweep 60. That is why SOR takes 60 sweeps here, not the $\ln 10^{-8}/\ln (\omega -1)=47$ the factor alone would give.

Conjugate gradients

CG treats $Au=b$ as minimizing $\frac{1}{2}{u}^{T}Au-b^{T}u$ and searches along directions that are conjugate, $p_i^{T}A{p}_j=0$, so no step undoes an earlier one. It needs only products with $A$. Its error bound depends on the condition number $\kappa$, the ratio of the largest eigenvalue of $A$ to the smallest: after $k$ steps the error, measured in the norm $\Vert e{\Vert }_A=\sqrt{e^{T}Ae}$, has shrunk by at least $2{\left(\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\right)}^k$. For the five-point matrix $\kappa ={\cot }^2(\pi h/2)$, which is 103 at $h=1/16$. The bound promises $10^{-8}$ within $\frac{1}{2}\sqrt{\kappa }\ln (2/10^{-8})=97$ iterations. It is a worst case over right-hand sides, and it bounds that error rather than the residual the stopping test measures, so it is a guide to the count, not a prediction of it. This right-hand side needs 25.

Like SOR, CG's count grows like $\sqrt{\kappa }\propto 1/h$. Unlike SOR, it needs no parameter tuned to the matrix.

The four residuals

The figure plots the residual of each method against the iteration count. The dashed lines are the theory's rates: $(\cos \pi h{)}^k$ for Jacobi and $({\cos }^2\pi h{)}^k$ for Gauss-Seidel. After the first few sweeps each measured line runs parallel to its dashed one.

Halving h

The last cell solves the same problem on 23 by 23 and 31 by 31 grids ($h=1/24$ and $1/32$) with Jacobi, SOR and CG, and plots the counts against $1/h$ on log axes. Doubling $1/h$ from 16 to 32 multiplies Jacobi's count by 4.01 and SOR's and CG's by 2.00 and 2.08. The dashed lines are the predictions, scaled to pass through the $h=1/16$ point: $\ln 10^{-8}/\ln \cos \pi h$ for Jacobi, $\ln 10^{-8}/\ln (\omega -1)$ for SOR and $\sqrt{\kappa }$ for CG.