Computational Physics

A7 Differential equations 585

solution which is in general unknown). In practice, a theoretical approximation of
this optimal value is used, or it is determined empirically using the Jacobi method.
This method is calledsuccessive over-relaxation.
Summarising, we can say that the iterative methods discussed in this section are
subject to slow mode relaxation problems; some more than others. The reason for
the slow relaxation is that the update of the value of the solution at the grid points
takes only nearest neighbour points into account and acts hence on a short range.
It will take many iterations in order to update a long-wavelength mode. In the next
subsections we shall consider ways to get round this problem.

The conjugate gradient method for elliptic differential equations

For ease of notation and for generality we denote the PDE problem as

Ax=b (A.96)

wherexijis the solution;xandbare considered asvectorswhose indices are the
grid points(i,j)(in two dimensions), andAis amatrixwith similar indices. When
solving Poisson’s equation,Ais the matrix corresponding to the discretised Laplace
operator, as defined in (A.87). Using this notation, the Jacobi relaxation method
can be represented as

xnij+^1 −xnij=

∑

kl

Aij,klxnkl−bij. (A.97)

The result on the left hand side can be viewed as the extent to which the solution
xnfails to satisfy the equation(A.96). The left hand side is called theresidualrof
the trial solutionxn.
We would like to find the valuexfor which the residual vanishes. It turns out
that this is possible using the conjugate gradients method. To this end, we consider
the function
f(x)=^12 x·Ax−b·x. (A.98)

The gradient of this function is

∇f(x)=Ax−b. (A.99)

Therefore, when we can minimise this equation, we have solved the matrix equation
(A.96). But we have already encountered a method which is very efficient at doing
this: the conjugate gradient method (see Appendix A4). We now give this algorithm
in pseudocode for the present problem.

h=b;(bis the right hand side of (A.96)

g=b;

r 2 =g·g;

WHILEr 2 is not small enough DO