Computational Physics

586 Appendix A

k=Ah;

λ=r 2 /(h·k);

g=g−λk;

r 2 ′=g·g;

γ=r′ 2 /r 2 ;

x=x+λh;

h=g+γh;

r 2 =r′ 2 ;

END WHILE.

The algorithm can be directly related to the discussion of the conjugate gradient
method in Appendix A4.
The convergence of the algorithm is fast:O(L)steps are needed, each of which
takesLdfloating point operations (as a step involves a sparse matrix–vector mul-
tiplication). We see that we only need a matrix–vector multiplication routine. The
algorithm is so simple to implement and so efficient that it should always be
preferred over the Gauss–Seidel method.

Fast Fourier transform methods

The first method uses the fast Fourier transform (FFT), which enables us to calculate
the Fourier transform of a function defined on a one-dimensional grid ofNpoints
in a number of steps which scales asNlogN, in contrast with the direct method
of evaluating all the Fourier sums, which takesO(N^2 )steps. The FFT method will
be described in Appendix A9. We now explain how the FFT method can be used
in order to solve PDEs of the type discussed in the previous section. On a two-
dimensional square grid of sizeL×L, the FFT requiresO(L^2 logL)steps. In the
following we assume periodic boundary conditions.
The idea behind using Fourier transforms for these equations is that the Laplace
operator is diagonal in Fourier space. For our model problem of the previous section,
we have
∑

jm

(∇^2 Dψjm)ei(kxj+kym)=(2 coskx+2 cosky− 4 )ψˆkxky, (A.100)

whereψˆis the Fourier transform ofψ. The wave numberskx,kyassume the values
2 πn/L. We see that acting with the Laplace operator onψcorresponds to mul-
tiplying by the factor 2 coskx−2 cosky−4ink-space. If we are to determine the
solution to the PDE∇D^2 ψ=−ρ, we perform a Fourier transform on the left and
right sides of this equation, and we arrive at

( 4 −2 coskx−2 cosky)ψˆkxky=ˆρkxky. (A.101)