4 Iterative methods

An implementation of the five-point stencil

$-4u_0+u_E+u_W+u_S+u_N=h^2{f}_0$

leads to a system of simultaneous equations in the unknowns. This system of equations can be dealt with using methods seen in HELM booklet  30, but here we show ways in which systematic iterative methods can be derived directly from the numerical stencil .

The general approach is as follows:

1. Start with an initial guess for the unknowns. Call this initial guess $u_{i,j}^0$ .
2. Use some means to improve the guess. Call the improvement $u_{i,j}^1$ .
3. And so on. In general we derive a new set of approximations $u_{i,j}^{n+1}$ in terms of the previous approximations $u_{i,j}^n$ .

4.1 Jacobi iteration

The approach we adopt here is to update the approximation at the centre of the stencil using the four old values around the edge of the stencil. That is

$-4u_0^{n+1}+u_E^n+u_W^n+u_S^n+u_N^n=h^2{f}_0$

rearranging this gives

$u_0^{n+1}=\frac{1}{4}\left(u_E^n+u_W^n+u_S^n+u_N^n-h^2{f}_0\right)$

The following Example uses the same data (rounded to four decimal places here) as in Example 6.

Example 7

Suppose that $u=u\left(x,y\right)$ satisfies Laplace’s equation

$u_{xx}+u_{yy}=0$

in the square region $0<x,y<1$ with $u=x^{2}y$ on the boundary. Assuming a mesh size of $h=\frac{1}{3}$ use the Jacobi iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Solution

Putting in the initial guesses for the four unknowns $u_{11},u_{12},u_{21},u_{22}$ we obtain the situation depicted below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

The first iteration involves using

$4u_0^1=u_E^0+u_W^0+u_N^0+u_S^0-h^2{f}_0$

where, in this case, $h^2{f}_0=0$ . So the first iteration gives us

$u_{11}^1=0.0000u_{21}^1=0.0833u_{12}^1=0.0278u_{22}^1=0.2778$

The second iteration begins by putting these new approximations to the interior values into the grid. This gives

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0.0278\hfill & \hfill 0.2778\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0833\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

We now apply $4u_0^2=u_E^1+u_W^1+u_N^1+u_S^1$ to obtain

$u_{11}^2=0.0278u_{21}^2=0.1528u_{12}^2=0.0972u_{22}^2=0.3056$

In practice, using a computer to carry out the arithmetic, we would continue iterating until the results settle down to a converged value. Using a computer spreadsheet, for example, we can see that a total of 15 iterations is enough to achieve results converged to four decimal places. We noted earlier that, to four decimal places, $u_{11}=0.0833$ , $u_{21}=0.1944$ , $u_{12}=0.1389$ and $u_{22}=0.3611$ .

The following Task uses the same data as the preceding Task (pages 23-24), except that we have rounded the boundary data to four decimal places instead of using the exact fractions.

Task!

Suppose that $u=u\left(x,y\right)$ satisfies Poisson’s equation

$u_{xx}+u_{yy}=-2$

in the square region $0<x,y<1$ with $u=xy$ on the boundary. Assuming a mesh size of $h=\frac{1}{3}$ use the Jacobi iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.3333\hfill & \hfill 0.6667\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Answer

The second iteration begins by putting these new approximations to the interior values into the grid. This gives

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.3333\hfill & \hfill 0.6667\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0.1389\hfill & \hfill 0.3889\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0.0556\hfill & \hfill 0.1389\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

We now perform the second iteration $4u_0^2=u_E^1+u_W^1+u_N^1+u_S^1-h^2{f}_0$ again, but with the new values. We obtain

In the case above 17 iterations are required to achieve results that have converged to 4 decimal places. We find that $u_{11}=0.2222$ , $u_{12}=0.3333$ , $u_{21}=0.3333$ and $u_{22}=0.5556$ .

4.2 Gauss-Seidel iteration

In the implementation of the Jacobi method we used old values for the southerly and westerly points when new values had already been calculated.

The Gauss-Seidel method uses the new values as soon as they are available. Stating this formally we have

$\begin{array}{ccc}\hfill & \multicolumn{2}{c}{\text{new values here}}\\ \hfill & \downarrow \hfill & \downarrow \hfill \\ \hfill u_0^{n+1}=\frac{1}{4}\left(u_E^n+\right.& u_W^{n+1}+\hfill & \left.u_S^{n+1}+u_N^n-h^2{f}_0\right)\hfill \end{array}$

Example 8 below uses the same data as Examples 6 and 7.

Example 8

Suppose that $u=u\left(x,y\right)$ satisfies Laplace’s equation

$u_{xx}+u_{yy}=0$

in the square region $0<x,y<1$ with $u=x^{2}y$ on the boundary. Assuming a mesh size of $h=\frac{1}{3}$ , use the Gauss-Seidel iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Solution

Putting in the initial guesses for the four unknowns we obtain the situation depicted below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

The first iteration involves using

$4u_0^1=u_E^0+u_W^1+u_N^0+u_S^1-h^2{f}_0$

where in this case $h^2{f}_0=0$ . So the first iteration gives us

The second iteration begins by putting these new approximations to the interior values into the grid. This gives

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.1111\hfill & \hfill 0.4444\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0.0278\hfill & \hfill 0.3056\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0833\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

We now apply $4u_0^2=u_E^1+u_W^2+u_N^1+u_S^2-h^2{f}_0$ to obtain

$u_{11}^2=0.0278u_{21}^2=0.1667u_{12}^2=0.1111u_{22}^2=0.3472$

(And, using a computer spreadsheet, for example, we can see that a total of 7 iterations is enough to achieve results converged to four decimal places. This compares well with the 15 iterations required by Jacobi in Example 7.)

Task!

Suppose that $u=u\left(x,y\right)$ satisfies Poisson’s equation

$u_{xx}+u_{yy}=-2$

in the square region $0<x,y<1$ with $u=xy$ on the boundary. Assuming a mesh size of $h=\frac{1}{3}$ use the Gauss-Seidel iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below.

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.3333\hfill & \hfill 0.6667\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Answer

The second iteration begins by putting these new approximations to the interior values into the grid. This gives

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.3333\hfill & \hfill 0.6667\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0.1528\hfill & \hfill 0.4653\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0.0556\hfill & \hfill 0.1528\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

We now apply $4u_0^2=u_E^1+u_W^2+u_N^1+u_S^2-h^2{f}_0$ again, but with the new values. We obtain

and we can write this information in the form

$\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.3333\hfill & \hfill 0.6667\hfill & \hfill 1.0000\hfill \\ \hfill 0.0000\hfill & \hfill 0.2882\hfill & \hfill 0.5330\hfill & \hfill 0.6667\hfill \\ \hfill 0.0000\hfill & \hfill 0.1319\hfill & \hfill 0.2882\hfill & \hfill 0.3333\hfill \\ \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Again, a computer can be used to continue iterating until convergence. This method applied to this Task needs 8 iterations to achieve 4 decimal place convergence, a fact which compares very well with the 17 required by the Jacobi method.

4.3 Convergence

We now summarise some important points

1. For the problems discussed in these pages, the Jacobi and Gauss-Seidel methods will always converge for any initial guesses $u_{ij}^0$ . (Of course, very poor initial guesses will result in more iterations being required.)
2. For a given problem and given starting guesses $u_{ij}^0$ , the Gauss-Seidel method will, in general, converge in fewer iterations than Jacobi . (That is, using the new, improved values as soon as they are available speeds up the process.)
3. One possible advantage with the Jacobi approach is that it can be parallelised , that is, it is in theory possible to do all the calculations for a given iteration simultaneously. In other words, everything we will need to know to carry out an iteration is known before the iteration begins. This is not the case with Gauss-Seidel in which during an iteration, most calculations use a result from within the current iteration. This advantage with Jacobi only manifests itself when using computers with a parallelisation option and for large problems.

Exercises

1. Suppose that $u=u\left(x,y\right)$ satisfies Laplace’s equation

   $u_{xx}+u_{yy}=0$

   in the square region $0<x,y<1$ . Assuming a mesh size of $h=\frac{1}{3}$ use the Jacobi iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below:

   $\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.2500\hfill & \hfill 0.7500\hfill & \hfill 1.0000\hfill \\ \hfill 0.4000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.8000\hfill \\ \hfill 0.8000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.4000\hfill \\ \hfill 0.0000\hfill & \hfill 0.7500\hfill & \hfill 0.2500\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

2. Suppose that $u=u\left(x,y\right)$ satisfies Laplace’s equation

   $u_{xx}+u_{yy}=0$

   in the square region $0<x,y<1$ . Assuming a mesh size of $h=\frac{1}{3}$ use the Gauss-Seidel iteration, with starting values $u_{ij}^0=0$ , to perform two iterations. The boundary data are as given in the schematic below.

   $\begin{array}{c}\hfill y\uparrow \hfill \\ \hfill 0.0000\hfill & \hfill 0.2500\hfill & \hfill 0.7500\hfill & \hfill 1.0000\hfill \\ \hfill 0.4000\hfill & \hfill u_{12}\hfill & \hfill u_{22}\hfill & \hfill 0.8000\hfill \\ \hfill 0.8000\hfill & \hfill u_{11}\hfill & \hfill u_{21}\hfill & \hfill 0.4000\hfill \\ \hfill 0.0000\hfill & \hfill 0.7500\hfill & \hfill 0.2500\hfill & \hfill 0.0000\hfill & \hfill \to x\hfill \end{array}$

Answer