FYTN03

HT08

Lecture Notes, FYTN03

Computational Physics

1 Introduction

Numerical methods in physics is simple. Just take

R your
P formula and replace every-
thing looking like δx with ∆x. This also implies → , ie.
Z b N
X
f (x)dx ≈ ∆xf (a + n∆x),
a n=0

with N ∆x = b − a.

When we do this we introduce an error. The smaller ∆x the smaller the (truncation)
error, but the required computational time also increases since we have to do a larger
number of function evaluations.

To understand the truncation error we typically simply consider the Taylor expansion:
∞
X ∆xn
f (x + ∆x) = f (n) (x).
n=0
n!

This is all you need to know.

Almost.

To reliably do numerical physics computations you need to have lots of experience

and learn lots of tricks. Especialy important is to understand tricks involving pseudo
random numbers. You also need to learn how to write programs. In this course you
will learn a lot of tricks. You will gain a beginning of experience. You will, however
not learn programming.

1
2 Errors. Interpolation and Extrapolation.

[NR: 1.3, 3.0, 3.1, 5.7 (3.2, 3.3)]

2.1 Errors

Numerical calculations are done in integer or floating-point arithmetic. Integer arith-

metic is exact, whereas floating-point arithmetic has roundoff errors.

The floating-point representation looks like


 s sign
s · M · 2p M mantissa

p integer exponent

Any real number can be represented this way, and the representation becomes unique
if we require eg. that 12 ≤ M < 1 for any non-zero number.

On the computer, each floating-point number becomes a string of bits, each bit being
0 or 1. Suppose the number of bits per floating-point number (the “wordlength”) is
32, with 1 sign bit, 8 bits for p and 23 bits for M . p and M may then (for example)
be represented as

p = x7 · 27 + x6 · 26 + ... + x0 · 20 − 27 xi ∈ {0, 1}
| {z }
−27 ,−27 +1,...,27 −1
−1
M = y1 · 2 + y2 · 2−2 + ... + y23 · 2−23 yi ∈ {0, 1}

Roundoff errors arise because the number of yi ’s is finite.

How large are the roundoff errors? The number 1 can be represented exactly (x7 =
x0 = y1 = 1, the other bits zero). The smallest number > 1 that can be represented
exactly is 1 + 2−22 (x7 = x0 = y1 = y23 = 1, the other bits zero). This gives
us an estimate of the typical relative precision ε of 32-bit floating-point numbers:
ε ∼ 2−22 ∼ 10−7 .

C has two data types for floating-point numbers, float and double. The size of a
float or a double is not specified by the language. To find out the size, one can use
the function sizeof which returns the number of bytes for a given type (1 byte = 8
bits).

2
In a C program you also have access to the precision, in you #include the standard
header file float.h: FLT_EPSILON gives the precision for float and DBL_EPSILON
gives the precision for double

Roundoff errors are especially troublesome when subtracting two numbers with small
relative difference. Consider eg. one of the solutions of the equation ax2 + bx + c = 0:
√
b − b2 − 4ac
x=−
2a
√ √
for ac ≪ b. Here the trick is simple. Just multiply with b + b2 − 4ac in both
denominator and nominator:
√ √
b− b2 − 4ac b + b2 − 4ac 2c
x=− × √ =− √ .
2a 2
b + b − 4ac b + b2 − 4ac

In addition to roundoff errors, most numerical calculations contain errors due to

approximations, such as discretization or truncation of an infinite series. Such errors
are called truncation errors. Truncation errors are errors that would persist even on
a hypothetical computer with infinite precision.

Example: Numerical derivation.

Suppose a derivative f ′ (x) is calculated by using

h2 ′′
f (x + h) = f (x) + hf ′ (x) + f (x) + ...
2
f (x + h) − f (x)
f ′ (x) = + O (h)
h

• The truncation error is εt ∼ h |f ′′ (x)|.

• The roundoff errors in f (x + h) and f (x) are & ε|f (x)|, where ε is the relative
floating-point precision. Therefore, a lower bound on the total roundoff error
εr is εr & ε|f (x)|/h.

Note that h large ⇒ εt large, whereas h small ⇒ εr large.

For the total relative error, we get the estimate

· ¸
εr + εt h |f ′′ | + ε |f | /h 1 ³p ′′
p ´2 p
′′
& = ′ h |f | − ε |f | /h + 2 ε |f f |
|f ′ | |f ′ | |f |

3
With an optimal choice of h (so that the square vanishes), this becomes
s
εr + εt √ |f ′′ f |
& ε
|f ′ | f ′2

So, the error is O(ε1/2 ) rather than O(ε), which does make a difference: ε = 10−7 ⇒
ε1/2 ≈ 3 · 10−4 .

How can we improve on this?

1. One way is to derive a better approximation of the derivate by direct manipu-

lation of the Taylor expansion.

h2 ′′
′ h3 ′′′
f (x + h) = f (x) + hf (x) + f (x) + f (3) + ...
2 3!
2 3
h h
f (x − h) = f (x) − hf ′ (x) + f ′′ (x) − f ′′′ (3) + ...
2 3!
2 4
f (x + h) − f (x − h) h h
= f ′ (x) + f ′′′ (x) + f (5) (x) + ...
2h 3! 5!
The truncation error is O(h2 ) for this “central difference”, which is one order
better than for the “forward difference” in the example above (as a result, the
total error scales as ε2/3 for the central difference, as can be easily verified).

2. Another possibility is to use Richardson extrapolation.

2.2 Richardson Extrapolation

Suppose we want to determine a quantity a0 that satisfies

a0 =lim a (h)
h→0

where h may be thought of as a step-size parameter. Suppose further that

1. the functional form of a(h) is known to be a (h) = a0 + a1 h2 + a2 h4 + ...

2. the value a (h) is known for h = h0 , 2h0 , 22 h0 , . . .

(a0 could be a derivate and a(h) the central-difference approximation).

4
We can then obtain an improved estimator â(h) of a0 in the following way:
½
a(h) = a0 + a1 h2 + a2 h4 + O(h6 ) ¡ ¢
2 4 6 ⇒ 4a (h) − a (2h) = 3a0 − 12a2 h4 + O h6
a(2h) = a0 + 4a1 h + 16a2 h + O (h )
4a (h) − a (2h) ¡ ¢
⇒ â (h) ≡ = a0 − 4a2 h4 + O h6
3

â(h) is an improved estimator of a0 because the truncation error is O(h4 ) instead of

O(h2 ).

We can now go on to eliminate the h4 term by forming

ˆ (h) ≡ 16â (h) − â (2h) = a0 + O h6

¡ ¢
â (verify this!)
15
ˆ
and so on. To calculate â(h) for one value h = h0 , we need to know a(h) for h = h0 ,
2h0 and 4h0 .

a(4h0 )
ց
a(2h0 ) → â(2h0 )
ց ց
ˆ 0)
a(h0 ) → â(h0 ) → â(h

Comments

The same procedure can be applied for more general versions of the assumptions 1
ˆ will then change.
and 2. The precise form of the expressions for â and â

Examples of methods that use this technique are Romberg’s integration method and
the Burlirsh-Stoer method for ordinary differential equations.

2.3 Interpolation and Extrapolation

Suppose we are given a set of data points (xi , yi ), i = 1, 2, . . . , N , and want to find
an approximating or interpolating function f (x; {cj }), where the cj ’s are parameters.
For a given functional form of f , the task then is to find optimal parameters cj .

Suppose f is taken to be a polynomial of degree M − 1,

f (x; {cj }) = c1 + c2 x + ... + cM xM −1

5
Whether or not it is possible to find cj ’s such that

M −1
 c1 + c2 x1 + ... + cM x1
 = y1
f (xi ; {cj }) = yi , i = 1, 2, ..., N ⇔ .
..

 c + c x + ... + c xM −1 = y
1 2 N M N N

depends on N and M — we have a linear system of N equations for the M unknown

parameters cj .

• M >N
This case is the least interesting one — there are too many parameters.

• M =N
In this case, there is a unique solution which is given by the Lagrange interpo-
lation formula
(x − x2 ) (x − x3 ) · · · (x − xN ) (x − x1 ) · · · (x − xN −1 )
f (x) = · y1 + ... + · yN
(x1 − x2 ) (x1 − x3 ) · · · (x1 − xN ) (xN − x1 ) · · · (xN − xN −1 )
N Q
X (x − xj )
= Q j6=i yi
i=1 j6=i (xi − xj )

provided that xi 6= xj whenever i 6= j (geometrically: through any two points

there is a unique straight line or first-degree polynomial; through any three
points there is unique second-degree polynomial; etc.).

• M <N
In general, there is no solution in this case, so one must approximate rather
than interpolate. This is often done by using the method of least squares; that
is, the cj ’s are determined by minimizing
N
X
[f (xi ; {cj }) − yi ]2
i=1

Taking f (x, {cj }) to be a polynomial is a convenient choice, for example, in numerical

integration. However, it is far from always the best choice when it comes to approx-
imating functions (consider, for example, the function (1 − x)−1 = 1 + x + x2 + . . .).
Two alternatives to polynomials are rational functions and cubic splines.

Rational functions

6
 • Interpolation.
Example: Assume that there are three data points (xi , yi ) and that f (x; {cj }) =
(x − c1 )/(c2 x + c3 ). Interpolation then amounts to solving a linear system of
three equations for the cj ’s:

f (xi ; {cj }) = yi ⇔ c1 + c2 xi yi + c3 yi = xi i = 1, 2, 3

In general we can use any rational polynomial

P
Pµ (x) 1 + µ1 pi xi
Rµν (x) = = Pν j
,
Qν (x) 0 qj x

where we choose µ + ν + 1 equal to the number of points N and arrive at the

equation system
µ ν
X X
1+ i
p i x k − yk qj xjk = 0
1 0

• Padé approximation. Suppose we want to approximate a function g(x) and

know g(x) as well as a number of derivatives g ′ (x), g ′′ (x), ... for x = 0. In
the Padé method, the approximating function f (x; {cj }) is rational, and the
parameters cj are determined so that

dk dk g (0)
f (x = 0; {cj }) = g (0) , f (x = 0; {c j }) =
dxk dxk
up to highest possible order k.

Cubic spline-interpolation

Let a = x1 < x2 < ... < xN = b. The cubic spline function is obtained by using one
cubic polynomial for each subinterval [xi , xi+1 ]. These polynomials are put together
in such a way that the resulting function as well as its first two derivatives become
continuous over the whole interval a 6 x 6 b.

This will give us N − 1 functions with four parameters each, ie. 4(N − 1) unknown.
We have 2(N − 1) equations forcing each function to go through the two connecting
points and 2(N − 2) equations for the continuity of the first and second derivatives.
This means we will have to specify two more conditions, eg. the third derivatives for
the last and first point.

Note that splines may severely misrepresent functions which contain kinks, or if the
points have some statistical fluctuations around function values. In the latter case
one should use function fitting rather than interpolation (see appendix).

7
3 Numerical integration (quadrature)

[NR: 4.0, 4.1, 4.2, 4.3, 4.5 (4.4, 4.6)]

We will discuss several methods for calculating integrals. For an integral

Zb
I= f (x) dx
a

all these methods can be written in the form

½
X Ai weights
I≈ Ai f (xi )
xi abscissas
i

As mentioned in the introduction, the more points we use, the better the precision
and the longer the computing time. So let’s try to choose the points, xi , as efficiently
as possible, and let’s see if we can find some tricks to allow us to reduce the number
of points without reducing the precision. Here we will investigate three different
possibilities:

1. Equidistant xi ’s (section 3.1).

2. Gaussian quadrature, where the xi ’s are zeros of polynomials (section 3.2).
3. Monte Carlo, where the xi ’s are randomly chosen (section 4.4).

Another possibility is to use methods for ordinary differential equations. This can be
done because I = y(b), where y(x) denotes the solution to the initial value problem
dy
= f (x) y(a) = 0
dx
This can be a good approach if the integral is dominated by a few sharp peaks.

3.1 Equidistant xi

3.1.1 Trapezoidal Rule, Simpson’s Rule,. . .

Consider
xN
Z
I= f (x) dx fi ≡ f (xi ) xi+1 = xi + h, i = 1, . . . , N − 1
x1

8
• The trapezoidal rule.
Use a linear approximation of f (x) on each subinterval [xi , xi+1 ]:
xZ
i+1
xi+1 − x x − xi h
f (x) ≈ fi + fi+1 ⇒ f (x) dx ≈ (fi + fi+1 )
h h 2
xi

From the expansion

¡ ¢ fi+1 − fi ¡ ¢
f (x) = fi + (x − xi ) f ′ (xi ) + O h2 = fi + (x − xi ) + O h2
h
it follows that the truncation error in f (x) is O(h2 ). This in turn suggests that
the integral has an error of O(h3 ) (the length of the interval is h). Adding up
the intervals to a fixed integration region will give us O(1/h) intervals and the
total error of the integral will be O(h2 )
• The Simpson rule.
This method is similar, but uses second-degree instead of first-degree polynomi-
als. This means that three points are needed for the interpolation, which makes
it necessary to look at intervals [xi , xi+2 ]. Lagrange’s interpolation formula and
integration give:
xZ
i+2
h
f (x) dx ≈ (fi + 4fi+1 + fi+2 )
3
xi

The same kind of error estimate as above indicates that the truncation error is
O(h4 ) in this case. It turns out, however, that the actual error is smaller than
that, O(h5 ), thanks to cancellation.

9
So far, we looked at subintervals. The full integral is obtained by adding the results
for the subintervals. Consider, for simplicity, the trapezoidal rule.

1. Closed formula. Apply the result above to each of the subintervals, which in
particular means that the end points x1 and xN are used. This gives
xN
Z
¡1 ¢
f (x) dx ≈ h f + f2 + ... + fN −1 + 12 fN
2 1
x1

The truncation error should not be worse than ∼ N h3 ∼ h2 (the total length
xN − x1 = (N − 1)h is fixed).
2. Open formula. Here, to avoid using the end points, the two end intervals are
Rtreated
1 √
differently. This is necessary if f has an end-point singularity (eg.
0
dx/ x), or has a limiting value at the endpoint (eg. sin(x)/x at 0). For the
end intervals, the following estimates can be used:
Zx2 xN
Z
¡ ¢ ¡ ¢
f (x) dx ≈ hf2 + O h2 f (x) dx ≈ hfN −1 + O h2
x1 xN −1

The fact that the estimates for these two subintervals are poorer does not affect
the scaling of the total truncation error, which remains O(h2 ) when adding two
O(h2 ) terms.

3.1.2 Romberg’s Method = Trapezoidal Rule + Richardson Extrapola-

tion

Romberg’s method is based on the Euler-Maclaurin summation formula, which we

state without proof. Let
T (h) = h( 21 f1 + f2 + ... + fN −1 + 12 fN )

10
be the closed trapezoidal rule. The Euler-Maclaurin summation formula says that
xN
Z
T (h) − f (x)dx = c1 h2 + c2 h4 + . . .
x1

where 
 B2k ¡ (2k−1) ¢

 ck = f (xN ) − f (2k−1) (x1 )


 (2k)!
X tn ∞
t
 Bn are Bernoulli numbers, defined by t = Bn


 e − 1 n=0 n!
 dn
¡ t ¢
 ie. Bn = dtn et −1 |t=0

In particular, this result shows that the truncation error of the trapezoidal rule con-
tains only even powers of h. If we use Richardson extrapolation, we will therefore
gain two powers of h at each step. The situation is exactly the same as for the
central-difference approximation of a derivative, so from the results in section 2.2, it
immediately follows that
xN
Z
4T (h) − T (2h) ¡ ¢
T̂ (h) = = f (x) dx + O h4
3
x1
xN
Z
ˆ 16T̂ (h) − T̂ (2h) ¡ ¢
T̂ (h) = = f (x) dx + O h6
15
x1

Comments

• The Euler-Maclaurin formula

– For fixed h, the sum on the RHS does not necessarily approach the LHS
as the number of terms tends to infinity (it is an asymptotic series).
– The formula is sometimes useful in the opposite direction — to calculate
a sum by replacing it with an integral.
– There is an analogous formula for the open trapezoidal rule.

• Trapezoidal rule + 1 Richardson step = Simpson’s rule

To see this, consider an arbitrary subinterval [xi , xi+2 ]:
½
T (2h) = h (fi + fi+2 )
T (h) = h ( 12 fi + fi+1 + 12 fi+2 )

11
 h
⇒ T̂ (h) = 34 T (h) − 13 T (2h) = . . . = (fi + 4fi+1 + fi+2 )
|3 {z }
Simpson’s formula

3.2 Gaussian Quadrature

Let f be a polynomial of degree < N . Lagrange’s interpolation formula tells us that

N
X Y x − xj
f (x) = Li (x)f (xi ) Li (x) =
i=1 j6=i
x i − xj

for all possible choices of x1 , . . . , xN such that xi 6= xj if i 6= j. This implies that, for
all such f , the integration formula

Zb N
X Zb
f (x) w (x) dx ≈ f (xi ) Li (x) w (x) dx
a i=1
|a {z }
≡Ai

is exact. Here, w(x) ≥ 0 is an arbitrary “weight” function and it can, of course,

be put to 1, but it can also be set to some other function. Eg. consider
√ a function
which has an integrable singularity
√ which is known to√behave like 1/ x as x → 0,
then we can use w(x) = 1/ x and replace f (x) by xf (x) which may be better
approximated by a polynomial.

In Gaussian quadrature, this formula becomes exact for all polynomials with degree
< 2N , by a careful choice of the xi ’s.

To see how this works, we need to define orthogonality. Using the same a, b and w(x)
as in the integration formula, we define the scalar product of two arbitrary functions
f and g by
Zb
hf |gi = f (x) g (x) w (x) dx ,
a

and we say that f and g are orthogonal if hf |gi = 0. For any given scalar product, it
is possible to construct a sequence of orthogonal polynomials ϕ0 , ϕ1 , ... with degrees
0, 1, . . . (ϕn has degree n) by using the Gram-Schmidt method:

• ϕ0 (x) = 1

12
 hϕ0 |xi
• ϕ1 (x) = x + aϕ0 (x) hϕ1 |ϕ0 i = 0 ⇒ a = −
hϕ0 |ϕ0 i

• ϕ2 (x) = x2 + bϕ1 (x) + cϕ0 (x)

hϕ1 |x2 i hϕ0 |x2 i
hϕ2 |ϕ1 i = hϕ2 |ϕ0 i = 0 ⇒ b = − and c = −
hϕ1 |ϕ1 i hϕ0 |ϕ0 i
• ...

• ϕi+1 (x) = (x − ai )ϕi (x) − bi ϕi−1 (x)

hxϕi |ϕi i hϕi |ϕi i

⇒ ai = and bi = ,
hϕi |ϕi i hϕi−1 |ϕi−1 i
• ...

By construction, ϕi+1 is orthogonal to ϕi and ϕi−1 . With induction we can then also
show that if ϕi is orthogonal to all ϕj with j < i, so is ϕi+1 . We have for n > 1

hϕi+1 |ϕi−n i = h(x − ai )ϕi − bi ϕi−1 |ϕi−n i = hxϕi |ϕi−n i = hϕi |xϕi−n i

which is zero since xϕi−n can be written as a linear combination of polynomials of

degree j < i.

In NR 4.5 you will find formulas for the orthogonal polynomials corresponding to
some useful weight functions.

Suppose that ϕN is the N th-degree polynomial obtained by the Gram-Schmidt pro-

cedure, for the particular scalar product defined above. The Gaussian quadrature-
theorem says that if x1 , ..., xN are taken to be the zeros of ϕN , then the integration
formula
Zb XN Zb
f (x) w (x) dx ≈ Ai f (xi ) Ai = Li (x) w (x) dx
a i=1 a

is exact for all polynomials with degree < 2N (it can be shown that a < xi < b for
all i).

“Proof”: Let f be an arbitrary polynomial with degree < 2N . We can then write
f = qϕN + r for some polynomials q and r with degree < N . It follows that (in
simplified notation) Z Z Z
fw = qϕN w + rw

13
 R
where qϕN w = 0 because q is a linear combination of ϕ0 , ..., ϕN −1 , all of which are
orthogonal to ϕN . The second term on the RHS can be written as
Z N
X
rw = Ai r (xi )
i=1

because we know that this formula is exact for polynomials with degree < N . But
N
X N
X N
X
Ai r (xi ) = Ai f (xi ) − Ai q (xi ) ϕN (xi ) ,
| {z }
i=1 i=1 i=1 0

where ϕN (xi ) = 0 because of the choice of xi ’s. This completes the proof.

Example: Determine x1 , x2 , A1 and A2 so that the formula

Z1
f (x) dx ≈ A1 f (x1 ) + A2 f (x2 )
−1

becomes exact for all polynomials with degree < 4.

Solution: Use Gaussian quadrature (with N = 2). The relevant scalar product is

Z1
hf |gi = f (x) g (x) dx (w (x) ≡ 1)
−1

The orthogonal polynomials are known in this case, so we do not have to carry out
the Gram-Schmidt procedure. They are called Legendre polynomials and the first
three are given by
¡ ¢
P0 (x) = 1; P1 (x) = x; P2 (x) = 12 3x2 − 1 ;

and the following are given by

(i + 1)Pi+1 (x) = (2i + 1)xPi (x) − iPi−1 (x)

The desired abscissas x1 and x2 are the zeros of P2 (x),

1 1
x1 = − √ and x2 = √ .
3 3

14
The weights A1 and A2 can be determined by simply performing the integrals in the
definition above. Another, slightly simpler, way is to use that the integration formula
must be exact for the functions f (x) = 1 and f (x) = x, which gives

 Z1




 A1 + A2 = dx = 2

−1
Z1 ⇒ A1 = A2 = 1


 1 1


 − √ A1 + √ A2 = xdx = 0
 3 3
−1

With these x1 , x2 , A1 and A2 , the integration formula becomes exact for all polyno-
mials with degree < 4.

4 Random Numbers and Monte Carlo

[NR: 7.0, 7.1, 7.2, 7.3, 7.6, 7.8 (only importance sampling in 7.8)]

Monte Carlo calculation is a widely used term that can mean different things. Com-
mon to such calculations is that random numbers are involved.

• Monte Carlo calculation can mean simulation of a process that indeed is stochas-
tic in nature (for example, scattering processes).

• But it can also be a calculation of an integral or a sum, in which the random

numbers serve just as a computational tool. This type of Monte Carlo calcu-
lation is common, for example, in statistical physics, where it can be used to
calculate ensemble averages.

The plan of this section is as follows. We begin with some basics on random variables
and probability theory. We then discuss methods for generating random numbers on
the computer. Finally, we discuss Monte Carlo integration and summation.

15
4.1 Random variables

4.1.1 Some Definitions

A random variable X is defined by its probability distribution (or frequency function)

p(x) ≥ 0 which has the following properties1
p(x)dx = P {x < X < x + dx} = the probability that x < X < x + dx
Zb
p (x) dx = P {a < X < b}
a
Z∞
p (x) dx = 1 (normalization)
−∞

The average hf (X)i of an arbitrary function f (X) is defined as

Z∞
hf (X)i = f (x) p (x) dx
−∞

Some standard quantities for characterizing the distribution of a random variable X

are:

• The mean
Z∞
hXi = xp (x) dx
−∞
which in general is not the same as the most probable value.

1
Throughout the text, P {· · ·} denotes the probability that the statement within curly brackets
is true.

16
 • The variance
2
σX = h(X − hXi)2 i = hX 2 i − 2hXi2 + hXi2 = hX 2 i − hXi2

which is a measure of the fluctuations in X. The square root of the variance,

σX , is the standard deviation.

• Higher-order moments hX n i (n = 1, 2, 3, . . .).

Example: The normal distribution with mean µ and variance σ 2 is given by

1 (x−µ)2
−
p (x) = √ e 2σ 2
2πσ 2
(verify that the mean and variance indeed are µ and σ 2 , respectively)

The simultaneous distribution p(x, y) of two random variables X and Y is called the
joint distribution of these variables, and is defined by

p (x, y) dxdy = P {x < X < x + dx and y < Y < y + dy} .

Some useful one-dimensional distributions that can be obtained from p(x, y) are:

• The marginal distributions in x and y:

Z∞ Z∞
p (x) = p (x, y) dy p (y) = p (x, y) dx
−∞ −∞

• The conditional probabilities p(x|y) and p(y|x).

p(x|y) is the probability of x given y. For each y, p(x|y) is a normalized distribution

in x. The function x 7−→ p(x, y) is proportional to p(x|y) but is not normalized. The
precise relation between these two functions of x is given by Bayes’ rule, which says
that
p (x, y) = p (x|y) p (y) ( = p(y|x)p(x) ) .

Two random variables X and Y are independent if

p (x, y) = p (x) p (y) .

If so, then p(x|y) = p(x, y)/p(y) = p(x) for all y.

17
Two random variables X and Y are uncorrelated if
h(X − hXi)(Y − hY i)i = hXY i − hXihY i = 0 .

Independence is stronger than uncorrelatedness; if X and Y are independent, it

immediately follows that
hXY i = hXihY i
so X and Y are then uncorrelated, too.

The converse is not true, as the following example shows.

Example of uncorrelated but dependent random variables.

Let X and Y be independent binary random numbers both with possible values 0
and 1, and assume that pX (0) = pX (1) = pY (0) = pY (1) = 1/2. Put U = X + Y and
V = |X − Y |. There are four equally probable states of this system:
X Y probability U V
0 0 1/4 0 0
0 1 1/4 1 1
1 0 1/4 1 1
1 1 1/4 2 0
The marginal distributions of U and V are:
   
pU (0) 1/4 µ ¶ µ ¶
 pU (1)  =  1/2  pV (0) 1/2
=
pV (1) 1/2
pU (2) 1/4
U and V are dependent, because
p (U = 0, V = 1) = 0 6= pU (0) pV (1) .
U and V are nevertheless uncorrelated, because


 hU V i = 41 · 0 + 41 · 1 + 14 · 1 + 41 · 0 = 1
2

hU i = 41 · 0 + 21 · 1 + 41 · 2 = 1 ⇒ hU V i = hU ihV i


 hV i = 1 · 0 + 1 · 1 = 1
2 2 2

4.1.2 Sums of Independent and Identically Distributed Random Vari-

ables

Suppose X1 , . . . , XN are independent and identically distributed (iid) random vari-

ables with mean µ and variance σ 2 . What is then the distribution of their sum? This
question has a surprisingly simple answer in the limit N → ∞.

18
This result, called the central limit theorem, says that the random variable
X1 + ... + XN − N µ
S̃N = √
σ N
becomes normally distributed with zero mean and unit variance as N → ∞.

“Proof”: Let us sketch how to prove this. For this purpose, it is convenient to
introduce the so-called characteristic function Φ(k) for the random variable
1
SN = (X1 + . . . + XN ) − µ ,
N
defined as
Φ(k) = heikSN i .
Φ(k) is the Fourier transform of the probability distribution pN (s) of SN , which means
that pN (s) can be obtained as the inverse Fourier transform of Φ(k),
Z∞ Z∞
1
Φ (k) = eiks pN (s)ds ⇒ pN (s) = e−iks Φ (k) dk
2π
−∞ −∞

For large N , ln Φ(k) can be estimated in the following way:

ik ik
ln Φ (k) = ln he N (X1 −µ) · · · e N (XN −µ) i
ik
= N ln he N (X−µ) i (the Xi ’s are iid)
µ ¶2
ik 1 ik ¡ ¢
= N lnh1 + (X − µ) + (X − µ)2 + O N −3 i
N 2 N
· ¸
k2 2 ¡ −3 ¢ k2σ2
= N ln 1 − σ + O N ∼ − N →∞
2N 2 2N
k2 σ 2
So, Φ(k) ∼ e− 2N as N → ∞, which gives
Z∞ ½ ¾
1 k2 σ 2 σ
pN (s) ∼ e−iks e− 2N dk = u = k√ =
2π 2N
−∞
√ Z∞
1 2N 2 +iu
√
= e−(u 2N s/σ)
du =
2π σ
−∞
√ Z∞ r
s N 2 s2 N
1 2N −(u+i )
= e σ 2 du ·e− 2σ2
2π σ
−∞
| {z
√
}
π

19
 s2
1 −
2σ 2 /N
pN (s) ∼ p e N →∞
2πσ 2 /N
This implies that the rescaled variable S̃N indeed is normally distributed with zero
mean and unit variance for large N .

It is very important to note that in this derivation we never looked at the precise form
of the probability distribution of the individual Xi ’s — the result holds irrespective of
the precise form of this distribution. This makes the central limit theorem extremely
useful.

4.1.3 Confidence intervals

Suppose we make N measurements of some quantity µ and that these can be viewed
as independent and identically distributed random variables Xi with
hXi i = µ
hXi2 i − hXi i2 = σ 2
for i = 1, . . . , N . The average
1
MN = (X1 + ... + XN )
N
provides an unbiased estimator of µ, since hMN i = µ. A not unreasonable error bar
on this estimate is MN ± σMN , where σMN is given by
2 2 1 X 2 2 σ2
σM = h(M N − hMN i) i = (hXi i − hXi i ) = .
N
N2 i N
However, there are two problems with this error estimate.

The first one is how to estimate the unknown quantity σMN . This is relatively easy.
If we put
" #2 Ã !2
1 X 1 X 1 X 2 1 X
YN = Xi − Xj = Xi − Xi
N i N j N i N i

then
1 X 2 1 X
hYN i = hXi i − 2 hXi Xj i =
N i N i,j
µ ¶
1 1 X 1 X N −1 2
= − 2 hXi2 i − 2 hXi Xj i= 2
σ = (N − 1)σM
N N i
| {z } N i6=j | {z } N N

σ 2 +µ2 µ2

20
This implies that
" #2
1 1 X 1 X
YN = Xi − Xj
N −1 N (N − 1) i
N j
2 2
provides an unbiased estimator of σM N
; that is, hYN i/(N − 1) = σM N
.

The second and more fundamental problem is how to assign a probabilistic meaning
to the error bar. This can be done by making use of the central limit theorem, which
says that S̃N = (MN − µ)/σMN becomes normally distributed with zero mean and
unit variance for large N . This means that
½ ¾ Z1
MN − µ 1 2
P {|MN − µ| < σMN } = P −1 < <1 ≈ √ e−t /2 dx ≈ 0.683
σMN 2π
−1

for large N . Hence, the probability is 68% that µ is in the interval MN ± σMN .
In other words, if we assign MN a statistical error of σMN , then we have chosen a
confidence level is 68%. An error of 2σMN corresponds to a confidence level of 95%.

4.1.4 Moments and Cumulants

Let Φ(k) = heikX i be the characteristic function for a random variable X. A Taylor
expansion of the exponential function yields
∞
X (ik)n n
Φ (k) = hX i
n=0
n!
where hX n i is the nth moment of X.

Closely related to the moments are the so-called cumulants, which are denoted by
hX n ic , n = 1, 2, . . .. The cumulants are defined by the Taylor expansion of ln Φ(k):
∞
X (ik)n
ln Φ (k) = hX n ic
n=1
n!
By expanding the LHS of this equation in terms of the moments, it is possible to see
the connection between cumulants are moments. It turns out that the nth cumulant
hX n ic can be expressed in terms of the first n moments. For the first two cumulants,
one finds (verify this)
hXic = hXi
hX 2 ic = hX 2 i − hXi2

21
It may seem unnecessary to introduce these quantities that anyhow can be expressed
in terms of the moments. However, the cumulants have two important properties
that make them quite useful:

• The cumulants provide a simple measure of how close a given distribution is to

the normal distribution, because all cumulants of order three and higher vanish
for the normal distribution.

• Cumulants of independent variables are additive; if X and Y are independent

and Z = X + Y , then
hZ n ic = hX n ic + hY n ic
for all n.

4.2 Transforming Random Numbers

On the computer, one typically has access to some random number generator that
delivers (pseudo) random numbers R uniformly distributed between 0 and 1,
½
1 0<r<1
p(r) =
0 otherwise

Such random numbers are sometimes called rectangularly distributed.

In this subsection, we look at how to obtain random numbers with other distributions,
assuming that we have uniformly distributed random numbers at our disposal. The
question of how to generate uniformly distributed random numbers is briefly discussed
in the next subsection.

4.2.1 The Transformation Method

Let X be an arbitrary random variable with distribution pX (x) and suppose Y =

f (X), where f is a monotonous function. What is then the distribution pY (y) of Y ?

Let PX (x) and PY (y) denote the cumulative distributions of X and Y , respectively;
that is,
Zx Zy
PX (x) = P {X ≤ x} = pX (t) dt and PY (y) = P {Y ≤ y} = pY (t) dt .
−∞ −∞

22
If y = f (x) is an increasing function of x, then
dPX dPY dy dy
PX (x) = PY (y) ⇒ pX (x) = = = pY (y)
dx dy dx dx
If, on the other hand, y = f (x) is a decreasing function of x, then
dy
PX (x) = 1 − PY (y) ⇒ pX (x) = −pY (y)
dx
This shows that the relation between pX (x) and pY (y) can be written as
¯ ¯
¯ dy ¯
pX (x) = pY (y) ¯¯ ¯¯
dx

for all monotonous f (and the expression remains valid in higher dimensions if |dy/dx|
is thought of as the Jacobian).

Assume now that X = R is uniformly distributed between 0 and 1. What should

then the transformation function f look like in order for Y = f (R) to have a certain
distribution pY (y)? Look for an increasing f . The cumulative distributions must
then satisfy
Zr
PY (y) = PR (r) = pR (t)dt = r ⇒ y = PY−1 (r) (0 ≤ r ≤ 1)
−∞

which means that f = PY−1 . Since 1 − R and R have the same distribution, we can
equally well take Y = PY−1 (1 − R). This makes Y a decreasing function of R.

This way of generating random numbers with different distributions will be called
the transformation method.

23
Example 1: The exponential distribution.

Suppose we have access to random numbers R uniformly distributed between 0 and

1 and want random numbers with the distribution
½ −λy
λe y≥0
p(y) = (λ > 0)
0 y<0

Use the transformation method. The first step is to calculate the cumulative distri-
bution
Zy
£ ¤y
P (y) = p (t) dt = −e−λt 0 = 1 − e−λy
−∞

The transformation r → y is then obtained by solving

1
P (y) = PR (r) = r ⇒ 1 − e−λy = r ⇒ y = − ln (1 − r)
λ
So, Y = −λ−1 ln(1 − R) has the desired distribution.

Example 2: The Box-Muller method for normally distributed random numbers.

The transformation method cannot be directly applied to the distribution

1 2
p(y) = √ e−y /2 (assume, for simplicity, zero mean and unit variance)
2π
because we cannot get a closed expression for P (y). It turns out, however, that this
problem can be circumvented by considering two independent variables Y1 and Y2
with the same distribution p(y),
1 −(y12 +y22 )/2
p (y1 , y2 ) = e .
2π
In polar coordinates (ρ, θ), this distribution becomes
¯ ¯
¯ ∂ (y1 , y2 ) ¯
p (ρ, θ) = p (y1 , y2 ) ¯
¯ ¯=ρ e−ρ2 /2 · 1 .
∂ (ρ, θ) ¯ | {z } |{z} 2π
| {z } p(ρ)
ρ p(θ)

The fact that this distribution factorizes, p(ρ, θ) = p(ρ)p(θ), implies that ρ and θ can
be generated independently.

24
 • The θ distribution is uniform. Take θ = 2πR1 , where R1 is uniformly distributed
between 0 and 1.

• Use the transformation method for ρ.

Zρ p
2 /2 2 /2
P (ρ) = te−t dt = 1 − e−ρ =r ⇒ ρ= −2 ln(1 − r)
0
p
So, ρ can be obtained as ρ = −2 ln(1 − R2 ), where R2 is another uniformly
distributed random number between 0 and 1.

A transformation back to Cartesian coordinates gives us two independent random

numbers with the desired distribution:
½ p
Y1 = p−2 ln(1 − R2 ) cos 2πR1
Y2 = −2 ln(1 − R2 ) sin 2πR1

4.2.2 The Accept/Reject Method

The transformation method is simple and convenient provided that a closed expres-
sion for P −1 can be obtained, but this is far from always the case. A more general
approach, which does not require knowledge of P −1 , is the accept/reject method. To
generate random numbers with a given distribution p(x), this method makes use of
an auxiliary function f0 (x), which must satisfy f0 (x) ≥ p(x). f0 should be chosen
so that it is easy to obtain random numbers with the distribution p0 (x) ∝ f0 (x) (f0
itself is not a normalized distribution). For a given f0 , the method is as follows (see
figure).

1. Draw a point (X, Y ) from the uniform distribution for the area under f0 . This
can be done in two steps:

(a) Draw X from the distribution p0 (x) ∝ f0 (x).

(b) Take Y = f0 (X)R, where R is uniformly distributed between 0 and 1.

2. Accept (X, Y ) if Y < p(X), and reject otherwise.

25
The accepted points (X, Y ) obtained this way will be uniformly distributed in the
region under p. This implies that the X component of the accepted points has the
distribution p(x).

Example: Kahn’s method for normally distributed random numbers.

Consider the distribution


 2 √1 e−x2 /2 x ≥ 0
p(x) = 2π
 0 x<0
(if a random number with this distribution is given a random sign, a normally dis-
tributed random number is obtained).

Use the accept/reject method with

 r
 2 −x+1/2
f0 (x) = e x≥0
 π
0 x<0
This choice of f0 is OK since
p(x) 2
= e−(x−1) /2 ≤ 1 ,
f0 (x)
and the corresponding probability distribution is
½ −x
f0 (x) e x≥0
p0 (x) = R ∞ =
−∞ 0
f (x)dx 0 x<0

We now follow the steps above:

26
 1. Draw (X, Y ) in the area under f0 . If R1 and R2 are uniformly distributed
between 0 and 1, we can put
(a) X = − ln R1 (exponential distribution; see previous example).
(b) Y = f0 (X)R2 .
2. Accept if
2 /2
Y < p(X) ⇔ R2 < e−(X−1)

This gives us a simple two-step “algorithm”:

1. Put X = − ln R1 .
2 /2
2. Accept X if R2 < e−(X−1) .

How efficient is this method? This depends crucially on the acceptance rate, which
is the area under p (which is 1) divided by the area under f0 . In this particular
example, we find an acceptance probability of
1
r Z∞ ≈ 0.76 .
2e
e−x dx
π
0

4.2.3 The Veto Algorithm

The ‘radioactive decay’ type of problems is very common. In this kind of problems
there is one variable t, which may be thought of as giving a kind of time axis along
which different events are ordered. The probability that ‘something will happen’ (a
nucleus decay) at time t is described by a function f (t), which is non-negative in the
range of t values to be studied. However, this naı̈ve probability is modified by the
additional requirement that something can only happen at time t if it did not happen
at earlier times t′ < t. (The original nucleus cannot decay once again if it already
did decay; possibly the decay products may decay in their turn, but that is another
question.)

The probability that nothing has happened by time t is expressed by the function
N (t) and the differential probability that something happens at time t by P(t). The
basic equation then is
dN
P(t) = − = f (t) N (t) .
dt

27
For simplicity, we shall assume that the process starts at time t = 0, with N (0) = 1.

The above equation can be solved easily if one notes that dN /N = d ln N :

½ Z t ¾ ½ Z t ¾
′ ′ ′ ′
N (t) = N (0) exp − f (t ) dt = exp − f (t ) dt ,
0 0

and thus ½ Z t ¾
′ ′
P(t) = f (t) exp − f (t ) dt .
0

With f (t) = c this is nothing but the textbook formula for radioactive decay. In
particular, at small times the correct decay probability, P(t), agrees well with the
input one, f (t), since the exponential factor is close to unity there. At larger t, the
exponential gives a dampening which ensures that the integral of P(t) never can
exceed unity, even if the integral of f (t) does. The exponential can be seen as the
probability that nothing happens between the original time 0 and the final time t. In
the parton-shower language of Quantum Chromo Dynamics, this corresponds to the
so-called Sudakov form factor.

If f (t) has a primitive function with a known inverse, it is easy to select t values
correctly:
Z t ½ Z t ¾
′ ′ ′ ′
P(t ) dt = N (0) − N (t) = 1 − exp − f (t ) dt = 1 − R ,
0 0

which has the solution

F (0) − F (t) = ln R =⇒ t = F −1 (F (0) − ln R) .

If f (t) is not sufficiently nice, one may again try to find a better function g(t), with
f (t) ≤ g(t) for all t ≥ 0. However to use the normal accept/reject method with
this g(t) would not work, since the method would not correctly take into account
the effects of the exponential term in P(t). Instead one may use the so-called veto
algorithm:

1. start with t0 = 0;

2. select ti = G−1 (G(ti−1 ) − ln R), i.e. according to g(t), but with the constraint
that ti > ti−1 ,

3. compare a (new) R with the ratio f (ti )/g(ti ); if f (ti )/g(ti ) ≤ R, then return to
point 2 for a new try;

28
 4. otherwise ti is retained as final answer.

It may not be apparent why this works. Consider, however, the various ways in which
one can select a specific time t. The probability that the first try works, t = t1 , i.e.
that no intermediate t values need be rejected, is given by
½ Z t ¾ ½ Z t ¾
′ ′ f (t) ′ ′
P0 (t) = exp − g(t ) dt g(t) = f (t) exp − g(t ) dt ,
0 g(t) 0

where the ratio f (t)/g(t) is the probability that t is accepted. Now consider the case
where one intermediate time t1 is rejected and t = t2 is only accepted in the second
step. This gives
Z t ½ Z t1 ¾ · ¸ ½ Z t ¾
′ ′ f (t1 ) ′ ′ f (t)
P1 (t) = dt1 exp − g(t ) dt g(t1 ) 1 − exp − g(t ) dt g(t) ,
0 0 g(t1 ) t1 g(t)

where the first exponential times g(t1 ) gives the probability that t1 is first selected,
the square brackets the probability that t1 is subsequently rejected, the following
piece the probability that t = t2 is selected when starting from t1 , and the final factor
that t is retained. The whole is to be integrated over all possible intermediate times
t1 . The exponentials together give an integral over the range from 0 to t, just as in
P0 , and the factor for the final step being accepted is also the same, so therefore one
finds that Z t
P1 (t) = P0 (t) dt1 [g(t1 ) − f (t1 )] .
0
This generalizes. In P2 one has to consider two intermediate times, 0 ≤ t1 ≤ t2 ≤
t3 = t, and so
Z t Z t
P2 (t) = P0 (t) dt1 [g(t1 ) − f (t1 )] dt2 [g(t2 ) − f (t2 )]
0 t1
µZ t ¶2
1 ′ ′ ′
= P0 (t) [g(t ) − f (t )] dt .
2 0

The last equality is most easily seen if one also considers the alternative region 0 ≤
t2 ≤ t1 ≤ t, where the rôles of t1 and t2 have just been interchanged, and the integral
therefore has the same value as in the region considered. Adding the two regions,
however, the integrals over t1 and t2 decouple, and become equal. In general, for
Pi , the i intermediate times can be ordered in i! different ways. Therefore the total
probability to accept t, in any step, is
∞ ∞ µZ t ¶i
X X 1 ′ ′ ′
P(t) = Pi (t) = P0 (t) [g(t ) − f (t )] dt
i=0 i=0
i! 0

29
 ½ Z t ¾ ½Z t ¾
′ ′ ′ ′ ′
= f (t) exp − g(t ) dt exp [g(t ) − f (t )] dt
0 0
½ Z t ¾
′ ′
= f (t) exp − f (t ) dt , (1)
0
which is the desired answer.

If the process is to be stopped at some scale tmax , i.e. if one would like to remain
with a fraction N (tmax ) of events where nothing happens at all, this is easy to include
in the veto algorithm: just iterate upwards in t at usual, but stop the process if no
allowed branching is found before tmax .

Usually f (t) is a function also of additional variables ~x. The methods of the preceding
section are easy to generalize if one can find a suitable function g(t, ~x) with f (t, ~x) ≤
g(t, ~x). The g(t) used in the veto algorithm is the integral of g(t, ~x) over ~x. Each time
a ti has been selected also an ~xi is picked, according to the conditional probability
g(~x|ti ), and the (t, ~x) point is accepted with probability f (ti , ~xi )/g(ti , ~xi ).

4.3 Random Number Generators

Random numbers are generally pseudo random numbers on the computer, obtained
by a deterministic algorithm. Whether or not these random numbers are “random
enough” depends on both the algorithm used and the problem at hand. The use
of deterministic algorithms may seem strange, but has two major advantages, speed
and reproducibility. Reproducibility can be a valuable property when debugging a
program.

Most programming languages provide some built-in random number generator, but
there are sometimes good reasons not to use these. One possible reason is portability.
Another possible reason is that built-in random number generators are not seldom
of quite poor quality.

4.3.1 The Linear Congruential Method

A widely used deterministic algorithm for generating random numbers is the linear
congruential method; many random number generators rely on this algorithm or
variants of it.

In its simplest version, this method uses a recursion formula of the form
ij+1 = aij + c mod m

30
to generate a sequence i1 , i2 , . . . of integers. Putting xj = ij /m, we obtain a nor-
malized sequence x1 , x2 , . . . of numbers between 0 and 1. The hope is that these
numbers, for a suitable choice of the integer parameters a, c and m (see table in
NR), will behave as approximately independent and uniformly distributed random
numbers.

Clearly, these xj ’s are not independent, and therefore it is important to test how
strong the correlations between the xj ’s are. One way to do this is as follows.

• x1 , x2 , . . . should be uniformly distributed on the unit interval from 0 to 1.

• (x1 , x2 ), (x3 , x4 ), . . . should be uniformly distributed on the unit square.

..
.

• (x1 , . . . , xn ), (xn+1 , . . . , x2n ), . . . should be uniformly distributed on the unit

cube in n dimensions.

For the linear congruential method with c = 0, it can be shown that all possible
n-tuples (xj+1 , . . . , xj+n ) fall onto one of at most (n! m)1/n different hyperplanes.
This number is not very large; if we, for example, take m = 232 and n = 10, then
(n! m)1/n ≈ 42. This shows that this method exhibits correlations that definitely may
cause problems in applications where many random numbers are needed.

It is also worth noting that the recursion formula above is periodic, with a period of
m or smaller, so m should be very large.

There are methods that are better than this one; for a discussion, see NR. Neverthe-
less, it should be kept in mind that random number generators are not perfect, and
they should be chosen with care in applications where lots of random numbers are
needed.

4.4 Monte Carlo Integration and Summation

Example: Monte Carlo calculation of π.

Consider the first quadrant of the unit circle. Its area (= π/4) can be written as
Z
I= dxdy
first quadrant

31
Let us see how this integral can be estimated by using random numbers.

Suppose we have N random points drawn from the uniform distribution on the unit
square. For each point i, we introduce a binary random variable χi such that χi = 1
if point i is in the first quadrant of the unit circle, and χi = 0 otherwise. The mean
of each χi is I, the area of the first quadrant (since the area of the square is 1). It
follows that
N
1 X 1
χi = {no. of points in the first quadrant} → I N →∞
N i=1 N

which gives us a method for estimating I (and thereby π).

Let us now generalize this by considering

Z
I = f (x)p(x)dx

where f (x) is an arbitrary function and p(x) is some probability distribution. Assume
that X1 , . . . , XN are independent random variables, all with the distribution p(x).
f (X1 ), . . . , f (XN ) are then independent and identically distributed random variables.
This means, according to the central limit theorem, that

N
1 X
IN = f (Xi )
N i=1

is approximately normally distributed for large N , with

Z
• mean hIN i = f (x)p(x)dx = I

2
• variance σN = (hf (X)2 i − hf (X)i2 )/N → 0 as N → ∞

32
Hence, IN may be used as an estimator of I for large N . The method can be
immediately generalized to higher dimensions.

Sums can be dealt with in a similar way. Consider

X
S= f (i)p(i)
i

where f (i) is an arbitrary function and p(i) is some discrete probability distribution.
If I1 , . . . , IN are independent random numbers with the distribution p(i), we can
estimate S by using
N
1 X
S ≈ SN = f (Ik )
N k=1
for large N .

4.4.1 Convergence Rate

What about the efficiency of Monte Carlo integration? Let us compare the efficency
of this method in one dimension, D = 1, with that of the Simpson rule.

Consider an integral over an interval of length L, and let Tε denote the amount of
computer time needed to achieve an accuracy of O(ε).

• Simpson’s rule
½
ε ∼ h4 (h step size)
⇒ Tε ∼ ε−1/4
Tε ∝ {no. of function values} ∼ L/h

• Monte Carlo
½
ε ∼ N −1/2 (N no. of points)
⇒ Tε ∼ ε−2
Tε ∝ N

This comparison shows that the convergence of the Monte Carlo method is typically
much slower than that of the Simpson rule for D = 1. The strength of the Monte
Carlo method is its generality. If, for example, D is large or the boundary of the
integration region is complex, there often are few alternatives to Monte Carlo.

33
4.4.2 Importance Sampling

A Monte Carlo calculation of a given integral

Z
I = f (x)dx

may use random numbers drawn from any probability distribution p(x) > 0. In fact,
with X1 , . . . , XN drawn from any given p(x) > 0, we may estimate I as
Z N
f (x) p 1 X f (Xi )
I= p(x)dx ≈ IN =
p(x) N i=1 p(Xi )

The mean of the estimator INp is, of course, independent of p, hINp i = I. The variance
h(INp )2 i − hINp i2 depends, by contrast, strongly on p. In order to have a reasonable
performance, it is therefore crucial to make a careful choice of p.

In principle, it is known what the optimal choice of p is, namely p(x) ∝ |f (x)| (see
NR 7.8). In practice, this is of little help because finding the proportionality constant
in this relation is as difficult as finding the integral we want to compute. However, it
is often possible to make an educated guess of p that is useful, although not perfect.

Example: The Ising model.

The Ising model is a simple model for ferromagnetism. The system consists of N
binary spin variables σi = ±1 that live on a lattice. In the absence of an external
magnetic field, the energy E of a configuration σ = (σ1 , . . . , σN ) is given by
X
E = −J σi σj
hiji

where the sum runs over all nearest-neighbor pairs on the lattice. The thermodynamic
behavior of the system is governed by the Boltzmann weight
p(σ) ∝ e−E(σ)/kT
where k is Boltzmann’s constant and T the temperature. The average of an observable
O, the total magnetization say, at temperature T is given by
X
O(σ)e−E(σ)/kT
X
hOi = O(σ)p(σ) = σ X
σ e−E(σ)/kT
σ

34
How can we calculate such an average?

• Exactly, by exhaustive enumeration of all possible states? No, because the

number of states, 2N , is too large; already for N = 100 (which is a very small
system), there are 2100 ∼ (103 )10 = 1030 possible states.

• “Naive” Monte Carlo? Here, we would draw configurations σ (1) , . . . , σ (J) from
the uniform distribution p0 (σ) = 1/2N = constant, and estimate
à J
! Ã J
!−1
1X (j) 1 X −E(σ(j) )/kT
hOi ≈ O(σ (j) )e−E(σ )/kT × e
J j=1 J j=1

But E is an extensive quantity, which means that there will be huge fluctuations
in e−E/kT . As a result, the variance is very large and the convergence very slow.

• Importance sampling. If instead we draw σ (1) , . . . , σ (J) from the Boltzmann

distribution p(σ), then we can use an estimate
J
1X
hOi ≈ O(σ (j) )
J j=1

that does not contain any Boltzmann factor e−E/kT . This is typically an enor-
mous improvement. The next question then is how to generate Boltzmann
distributed configurations. A widely used method for this is the Metropolis
algorithm.

4.5 The Metropolis Algorithm

Consider a system with state or configuration space T , and suppose we want to

sample some distribution p̃(σ), σ ∈ T . For simplicity, we assume T to be discrete (as
in the Ising model).

The Metropolis algorithm can be thought of as a guided random walk in the state
space T ,
σ1 → σ2 → σ3 → . . .
Here, σn denotes the state of the system at “time” n. The guidance is such that the
probability distribution pn of σn approaches p̃ for large n; that is,

lim pn (σ) = p̃(σ) . (2)

n→∞

35
This is meant to hold irrespective of what the initial distribution p1 is.

A fundamental property of the Metropolis algorithm is that pn+1 is entirely deter-

mined by pn ; in order to determine pn+1 , we do not have to know where the system
was at time n − 1, n − 2, . . .. A stochastic process with this property is called a
Markov chain. For a Markov chain, the time evolution can be described in terms of a
transition matrix W (σ, σ ′ ); W (σ, σ ′ ) being the conditional probability of finding the
system in state σ at time n + 1, given that it was in state σ ′ at time n. Another
important property of the Metropolis algorithm is that the transition matrix W (σ, σ ′ )
does not change with time n.

These two properties imply that the time evolution of pn is given by a simple vector-
matrix equation, X
pn+1 (σ) = W (σ, σ ′ )pn (σ ′ ) (3)
σ′

with a constant (n independent) matrix W .

The key question now is how to ensure that pn → p̃ as n → ∞, equation (1). Useful
information about this can be obtained from the theory for general Markov chains
with constant transition matrices (“stationary” Markov chains). For a general process
of this type, it can be shown that equation (1) does hold independent of the initial
distribution if the following two conditions are met:

1. The distribution p̃ is stationary. This means that pn = p̃ ⇒ pn+1 = p̃. An-

other way to say this is that p̃ should be an eigenvector of the matrix W with
eigenvalue 1; that is,
X
p̃(σ) = W (σ, σ ′ )p̃(σ ′ ) (for all σ)
σ′

2. The process is ergodic. Loosely speaking, this means that each state can be
reached from each other state. A somewhat more precise formulation of this
condition can be found in the hand-out for the first computer exercise.

We are not going to prove that these two requirements are sufficient to ensure that
equation (1) holds, but to give an idea of how it works, we will prove two weaker
statements. For this purpose, we need to define the distance between two arbitrary
distributions pa and pb , which can be taken as
X
kpa − pb k = |pa (σ) − pb (σ)|
σ

36
Statement 1: Suppose p̃ is stationary. Then the distance kpn − p̃k is a non-increasing
function of n.

Proof: Using equation (2) and that p̃ is stationary, we obtain

X
kpn+1 − p̃k = | pn+1 (σ) − p̃(σ)|
σ
X ¯¯ X ¯
′ ′ ′ ¯
= ¯ W (σ, σ )(pn (σ ) − p̃(σ ))¯
σ σ′
XX
6 W (σ, σ ′ )|pn (σ ′ ) − p̃(σ ′ )|
σ σ′
³X ´
= k pn − p̃k W (σ, σ ′ ) = 1
σ

Statement 2: Suppose, in addition to the stationarity of p̃, that W (σ, σ ′ ) > 0 for all
σ, σ ′ and that pn 6= p̃. Then the inequality above is strict,

k pn+1 − p̃k < k pn − p̃k

Proof: That pn 6= p̃ implies that pn (σ ′ ) − p̃(σ ′ ) takes on both positive and negative
values, since pn and p̃ both are normalized distributions. The same must then be
true for W (σ, σ ′ )(pn (σ ′ ) − p̃(σ ′ )), since W (σ, σ ′ ) > 0 for all σ, σ ′ . From this follows
that the inequality must be strict.

The Metropolis algorithm provides a simple and general way to ensure that condition
1 above is met. This is achieved by designing the basic update of the system in such
a way that detailed balance is fulfilled, a condition that is stronger than the condition
1 above. How this is done is discussed in the hand-out for the first computer exercise.

Comments

For simplicity, we have here considered discrete systems. The Metropolis algorithm
can be easily applied to continuous systems, too.

The fact that the states generated by the Metropolis algorithm are not independent
makes this method fundamentally different from methods such as the simple trans-
formation and accept/reject methods. Methods like these two are sometimes called
static, and methods like the Metropolis algorithm are then called dynamic.

37
5 Optimization (minimization/maximization)

[NR: 9.6, 10.0, 10.1, 10.4, 10.5, 10.6, 10.9]

Optimization is a wide field. Sometimes there exists a well-established method that

can be used in a black-box manner, but many optimization problems are true chal-
lenges.

Suppose we want to minimize some function f (x) in D dimensions, x ∈ RD . How to

proceed depends on a number of things, such as

• What is the dimensionality D?

• Are there constraints on x?

• Is the problem linear or non-linear?

• Is f such that minimization is a smooth downhill process, or are there traps in

the form of local minima?

• Do we want the global minimum of f , or is it sufficient to make a local mini-

mization?

• Do we have access to derivatives of f ?

In this section, we first look at a simple and general scheme called simulated annealing
that can be tried if the aim is to make a global optimization (section 5.1). We then
discuss a few different methods for local optimization; first, the downhill simplex
method without derivatives (section 5.2), and then the conjugate-gradient and quasi-
Newton methods that do use derivatives (sections 5.3.2 and 5.3.3). We end with a
few words on optimization in the presence of constraints (section 5.4).

But first of all, a general remark on precision. Suppose we want to minimize a function
f (x) (assume, for simplicity, D = 1). A Taylor expansion about the minimum xmin
gives
f (x) ≈ f (xmin ) + 12 (x − xmin )2 f ′′ (xmin ) (x near xmin )
since f ′ (xmin ) = 0. Near xmin , the roundoff error in f (x) is & ε|f (xmin )|, where ε is
the relative floating-point precision. If now |x − xmin | is so small that |f (x) − f (xmin )|
is comparable to the roundoff error, then we cannot expect to be able to come closer
to xmin , irrespective of what search method we use. This gives an estimate of the
smallest possible error, ∆x, in xmin :

38
 s
√ 2f (xmin )
1
2
(∆x)2 f ′′ (xmin ) & ε |f (xmin )| ⇒ ∆x & ε
f ′′ (xmin )

This shows that the error in xmin (& ε1/2 ) is typically much larger than that in f (xmin )
(& ε).

5.1 Global Optimization. Simulated annealing.

For functions with many local minima, global minimization is generally very hard,
because the system tends to get trapped in local minima. The simulated-annealing
method is an attempt to circumvent this problem by using Metropolis dynamics. The
function we want to minimize is then thought of as an energy E. With Metropolis
dynamics, steps upwards in E do occur (probability e−∆E/kT ), which is needed in
order for the system to be able to escape from local minima. Note also that the
Boltzmann factor ∝ e−E/kT gives a high statistical weight to low-E states at low
temperature.

In a simulated-annealing run, the “temperature” T serves as a control parameter.

The calculations are started at a high temperature where the mobility of the system
is high. The temperature is then gradually decreased till the system freezes. The
hope is that the final frozen state will be the global energy minimum.

Of course, this will not always be the case, so it is essential to repeat the experiment
for many different initial conditions. A key parameter in the simulations is the rate
of cooling. If the cooling is too rapid, it is likely that the system gets stuck in a local
minimum with non-minimal E.

The simulated-annealing method is young compared to the Metropolis method (pro-

posed in 1983 and 1953, respectively) but has been applied to a wide range of physical
and non-physical problems. Non-physical applications include various combinatorial

39
optimization problems, such as the traveling salesman problem. In this problem,
there are N cities to be visited, and each city is to be visited precisely once. The task
is to minimize the distance that the salesman has to travel. Solving this problem ex-
actly is possible only for small N , because the number of alternative routes, (N − 1)!,
grows rapidly with N . With simulated annealing it becomes possible to study larger
N . A simulated-annealing program for this problem can be found in NR 10.9.

5.2 Local Optimization without Derivatives.

The Downhill Simplex Method

When searching for the minimum of a function f (x) in D dimensions, it is of help to

know the gradient ∇f (x), because −∇f (x) is the direction in which the decrease of
f (x) is fastest (locally). In what direction should one search if the gradient ∇f (x)
is not available?

A simple and general scheme that does not rely on information on the gradient is the
downhill simplex method. A simplex is a geometrical object with one more vertex
than dimension: a line segment for D = 1, a triangle for D = 2, a tetrahedron for
D = 3, and so on. In the downhill simplex method, the simplex is a dynamic object
that can grow and shrink. When it reaches a minimum, it gives up and shrinks down
around it.

The method proceeds iteratively, starting from some simplex with vertices x1 , . . . , xD+1
which we assume ordered so that f (xD+1 ) ≤ f (xD ) ≤ . . . ≤ f (x1 ). The elementary
move can be seen as an attempt to improve the worst point, x1 , and is as follows.

Calculate the center of the face defined by the points x2 , . . . , xD+1 ,

D+1
1 X
xmean ≡ xi
D i=2

Since all these D points are better than x1 , it makes sense to try to move x1 in the
direction of xmean . Therefore, the next step is to reflect x1 across this face to

xa = xmean + (xmean − x1 )

40
Whether this point is accepted or not depends on the value of f (xa ):

• If f (xD+1 ) < f (xa ) < f (x2 ), replace x1 by xa .

• If f (xa ) < f (xD+1 ), try a larger step (the direction seems good) to
xb = xmean + 2(xmean − x1 )
x1 is then replaced by the best of xa and xb .
• If f (xa ) > f (x2 ), try instead a smaller step to
xc = xmean + 12 (xmean − x1 )
If f (xc ) < f (x2 ), replace x1 by xc . Otherwise, try an even smaller step to
xd = xmean − 12 (xmean − x1 )
If f (xd ) < f (x2 ), replace x1 by xd . If this attempt also fails, the simplex is too
large to give a good idea of what direction to choose. Therefore, we now give
up and shrink all the vertices towards the best one
xi → xi + 12 (xD+1 − xi ) i = 1, . . . , D

This update is iterated untill the values are no longer improving, according to some
stopping criterion.

41
5.3 Local Optimization with Derivatives

5.3.1 Successive Line Minimizations

Many methods for multidimensional minimization are based on successive line mini-
mizations. Usually, the gradient of the function is used to decide on what lines to be
considered. The actual minimization along each of these lines may or may not use
derivatives.

For one-dimensional minimization, there exist methods that are relatively robust and
fast. Assume that we have such a method at our disposal and want to minimize a
function f (x), x ∈ RD . We may then proceed as follows.

1. Pick a starting point x0 and a direction h0 .

2. Determine λ0 by minimization of f along the line λ 7→ x0 + λh0 .

3. Put x1 = x0 + λ0 h0 and select a new direction h1 .

4. Determine λ1 by minimization of f along the line λ 7→ x1 + λh1 .

5. . . .

But what directions hi should we use?

• The D basis vectors êi ?

No, we don’t want to be restricted to these directions (see figure in NR 10.5).

• hi = −∇f (xi )?
This is called steepest descent and may seem like a natural choice. However,
this method has an unwanted property that tends to make it inefficient (see
figure in NR 10.6): consecutive directions hi and hi+1 are, by construction,
orthogonal. To see this, note that xi+1 is obtained by minimizing f along the
line λ 7→ xi + λhi , which implies that

df (xi + λhi ) ¯¯ d(xi + λhi )

0= ¯ = ∇f (xi+1 ) · = (−hi+1 ) · hi
dλ λ=λi dλ

Conjugate directions

42
Near minima, quadratic approximations are generally good. Let us therefore consider
a quadratic function

 x and b ∈ RD
1
f (x) = 2 x · Ax − b · x + c A symmetric, positive definite D × D matrix

c number
(A is assumed positive definite so that f has a minimum). If we were to minimize
this f by successive line minimizations, how should we then choose the directions hi ?

When minimizing f , we are searching for an x at which all the components of the
vector ∇f (x) = Ax − b vanish. By construction, one component of ∇f (xi+1 ) has
to be zero, ∇f (xi+1 ) · hi = 0. Can we by a suitable choice of the hi ’s see to that
∇f (xi+1 ) · hi−1 = 0, too?
hi−1 · ∇f (xi+1 ) = hi−1 · [A(xi + λi hi ) − b] =
= hi−1 · [∇f (xi ) + λi Ahi ]
= λi hi−1 · Ahi
So,
hi−1 · ∇f (xi+1 ) = 0 if hi−1 · Ahi = 0
hi−1 and hi are called conjugate, or A orthogonal, if this condition is fulfilled.

5.3.2 The Conjugate Gradient Method

This is a method based on successive line minimizations in which the directions hi

are conjugate if the function is quadratic. This makes the method efficient for such
functions. For a general function, the method is expected to work well if we are
sufficiently close to the minimum.

The algorithm is simple to formulate. For convenience, put gi = −∇f (xi ), where f is
the function to be minimized. Let x0 be an arbitrary starting point and take h0 = g0 .
Pairs (x1 , h1 ), (x2 , h2 ),. . . are then generated by using the recursion formulas:
(
xi+1 = xi + λi hi where λi is determined by line minimization
gi+1 · gi+1
hi+1 = gi+1 + γi hi where γi =
gi · gi
(with γi = 0 this would be the steepest descent method).

The choice of the parameter γi is such that if

f (x) = 12 x · Ax − b · x + c

43
(as above), then hi · Ahi+1 = 0 and also hi+1 · Ahi = 0, since A is symmetric.

In fact we will show something more general, namely that for such f , the following
orthogonality relationships hold

hi · Ahj = 0 (i 6= j)
gi · gj = 0 (i =6 j)
gi · hj = 0 (i > j)

and that we can write

gi · hi gi · gi
λi = = (i = 0, 1, . . .)
hi · Ahi hi · Ahi

We will prove this using induction, i.e. assuming it holds for all i, j ≤ n, we will show
that it holds also for i, j ≤ n + 1.

First we look at some general relationships, noting that hi · gi+1 = 0 by construction.

So, we have

0 = hi ·gi+1 = hi ·(b−Axi+1 ) = hi ·(b−Axi −λi Ahi ) = hi ·(gi −λi Ahi ) = hi ·gi −λi hi ·Ahi

and
gi · hi
λi = .
hi · Ahi
Also,
gi · hi = gi · gi + γi−1 gi · hi−1 = gi · gi
so we can write
gi − gi+1
0 = hi+1 · Ahi = (gi+1 + γi hi ) ·
λi
where we know that gi+1 · gi = 0, hi · gi+1 = 0 and hi · gi = gi · gi , so
1 gi+1 · gi+1
0= (γi gi · gi − gi+1 · gi+1 ) ⇒ γi =
λi gi · gi
Now, assuming that gi · gj = hi · Ahj = 0 for all i 6= j with i, j ≤ n (this is certainly
true for n = 1) we need to show that it is true for all i 6= j with i, j ≤ n + 1:

gn+1 · gi = gn · gi − λn hn · Agi = −λn hn · Agi .

For i = 0 and i = n this is zero, otherwise we have

= −λn hn · A(hi − γi−1 hi−1 ) = 0.

44
Similarly we have for hn+1 · Ahi , which is zero for i = n, otherwise
gi − gi+1
hn+1 · Ahi = gn+1 · Ahi + γn hn · Ahi = gn+1 · = 0.
λi

We also need to show that everything holds for n = 1:

Clearly g1 · h0 = g1 · g0 = 0. Also
g0 − g1 1
h0 · Ah1 = h1 · Ah0 = (g1 + γ0 g0 ) · = (γ0 g0 · g0 − g1 · g1 ) = 0 .
λ0 λ0

In particular, this shows that the vectors gi are pairwise orthogonal. This can hold
for at most D different gi 6= 0 in D dimensions. Hence, gi = 0 for some i ≤ D. But if
gi = 0, then xi is the solution we want. This means that the method needs at most D
steps to find the minimum. This holds when the function is quadratic. For a general
f the orthogonality relations are at best approximate, and there is no guarantee of
convergence within a finite number of steps. In this case, the recursion formulas are
iterated untill some stopping criterion is fulfilled. Note that the algorithm is written
in such a way that it can be directly applied to general functions f .

Solution of linear equations by the conjugate gradient method

Suppose we want to solve a linear equation system

Ax = b

where A is a symmetric, positive definite matrix. One approach to this problem is

to make use of the fact that the solution must be the minimum xmin of the function

f (x) = 21 x · Ax − b · x ,

because xmin satisfies

∇f (xmin ) = Axmin − b = 0
So, in principle, we can solve our linear algebra problem by using a minimization
method such as the conjugate-gradient method. Is this a good approach? It can be
computationally convenient if the matrix A is large but sparse (that is, A has many
elements but relatively few of them are nonzero). The reason that the conjugate
gradient method is interesting in this case is that the matrix A enters the algorithm
only through expressions of the type A × (vector), which makes it relatively easy to
take advantage of the sparseness of the matrix.

45
5.3.3 The Quasi-Newton Method

Minimizing f (x) can be viewed as solving ∇f (x) = 0, which is a system of D generally

non-linear equations, 
 ∂1 f (x) = 0

..
 .
 ∂ f (x) = 0
D

Suppose we are at some point x. A Taylor expansion about this point gives

f (x′ ) = f (x) + (x′ − x) · ∇f (x) + 12 (x′ − x) · A(x)(x′ − x) + . . .

where A(x) is the so-called Hessian matrix with elements Aij (x) = ∂i ∂j f (x). By
taking the gradient with respect to the primed variables of both sides of this equation,
we obtain
∇f (x′ ) = 0 + ∇f (x) + A(x)(x′ − x) + . . .
From this we see that if the omitted higher-order terms can be neglected, and if x′
is given by
x′ − x = −A−1 (x) · ∇f (x) ,
then we have ∇f (x′ ) = 0. The Newton method would be to iterate this equation till
some stopping criterion is fulfilled.2

This method has, however, two disadvantages: first, the second derivatives ∂i ∂j f
are needed; and second, the equation system A(x)y = ∇f (x) must be solved for
y = A−1 (x) · ∇f (x) at each step. Usually, this makes the method impractical, but
there are exceptions, like the Levenberg-Marquardt method for χ2 minimization.

The so-called quasi-Newton method circumvents these two problems. This method
is based on successive line minimizations and can be written as

xi+1 = xi − λi Hi · ∇f (xi ) ,

where the parameter λi is determined by line minimization, and Hi is a D ×D matrix

constructed so that Hi → A−1 as i → ∞. For details, see NR 10.7.

The quasi-Newton method is often comparable in efficiency to the conjugate-gradient

method, but requires more memory if D is large; the memory requirement scales as
D2 for the quasi-Newton method and as D for the conjugate-gradient method.
2
The “standard” Newton method for a one-dimensional problem g(x) = 0 is given by the recur-
sion formula x′ − x = −g(x)/g ′ (x).

46
5.4 Constrained minimization

5.4.1 Linear optimization/linear programming

Constrained minimization with a linear cost function and linear constraints is referred
to as linear optimization or linear programming. For such problems there is a well-
established method called the Simplex method, which is described in NR 10.8.

For a simple example of such a problem, consider the cost function f (x1 , x2 ) = x1 +x2
with the constraints

2x1 + x2 ≤ 2 x1 + 2x2 ≤ 2 x1 , x 2 ≥ 0

Since the gradient ∇f = (1, 1) is constant, the desired minimum must be somewhere
on the boundary, and it is easy to see that the solution is (x1 , x2 ) = (0, 0).

5.4.2 Lagrange multiplier

Suppose we want to minimize a function f (x) subject to the constraint u(x) = 0,

x ∈ RD . One possible approach is to look for stationary points of the auxiliary
function f˜(x, λ) = f (x)+λu(x), where λ is called a Lagrange multiplier. A stationary
point of this function is a solution of the D + 1 equations

0 = ∇x f˜(x, λ) = ∇f (x) + λ∇u(x)

∂ f˜
0 = = u(x)
∂λ
The last equation is just the constraint, which says that x must fall onto the hy-
persurface defined by u(x) = 0. The gradient ∇u(x) is normal to this surface. The
first D equations imply that ∇f (x) is parallel (or anti-parallel) to ∇u(x). If t̂ is an

47
arbitrary unit vector in the tangential plane of the surface u(x) = 0, it follows that

∂f
= t̂ · ∇f = 0 .
∂t

This means that the derivative of f is zero in every “allowed” direction, which is
precisely what we want.

This method is widely used in analytical calculations. As a numerical method, it

has the disadvantage that a system of generally non-linear equations must be solved,
which can be tricky.

5.4.3 Soft constraints

Consider the same task once more; minimize f (x) subject to the constraint u(x) = 0.
A method that often works better numerically than the previous one is to introduce
the constraint in a soft manner, by forming the auxiliary function

f˜(x, Λ) = f (x) + Λu(x)2

This function is to be minimized. Numerically, this is generally more convenient than

solving a system of equations. The desired solution is obtained by extrapolating
results for large Λ, xΛ , to Λ = ∞, x = limΛ→∞ xΛ . Note that the penalty term
Λu(x)2 grows large as Λ → ∞ unless x satisfies u(x) = 0.

6 Ordinary Differential Equations (ODE)

[NR: 16.0, 16.1, 16.3, 16.6 (16.2, 16.4, 16.7)]

48
6.1 Introduction

This section deals with numerical integration of ODEs. We will discuss methods for
solving systems of first-order ODEs

dy1

 = f1 (x, y1 , . . . , yn )
 dx

..
 .

 dy n
 = fn (x, y1 , . . . , yn )
dx
for given initial values of all the components yi at some “time” x0 (for a discussion
of boundary value problems, see chapter 17 of NR). In vector notation, this initial
value problem can be written as
dy
= f (x, y) y(x0 ) = y0 .
dx

Looking only at first-order ODEs is not a strong limitation, because any nth order
ODE of the form µ ¶
dn y dy dn−1 y
= f x, y, , . . . , n−1
dxn dx dx
can be written as a system of first-order ODEs by the transformation

dy1
y1 = y 

 = y2


 dx
dy 
 dy2

y2 = = y3
dx ⇒ dx
.. ..
. 

 .
n−1 
d y 
 dyn
yn = n−1 

 = f (x, y1 , . . . , yn )
dx dx

For convenience, we will often assume that n = 1. This does not mean that these
methods work only for n = 1; on the contrary, generalizing to n > 1 is typically easy.

Example: Consider the one-dimensional Newton equation

d2 x
m = −V ′ (x)
dt2
Putting x1 = x and x2 = dx/dt, we obtain the first-order ODEs
½
dx1 /dt = x2
dx2 /dt = −V ′ (x1 )/m

49
which are the corresponding Hamilton equations.

(Hamilton’s equations of motion are given by dx/dt = ∂H/∂p and dp/dt = −∂H/∂x,
where p = m dx/dt and H = p2 /2m + V (x).)

6.2 Euler’s Method

Consider the first-order ODE

dy
= f (x, y) .
dx
To solve this equation numerically, we discretize x using a step size h:
xn = x0 + nh n = 0, 1, . . .
yn = y (xn )

Approximating the derivative with a simple forward difference,

dy ¯¯ yn+1 − yn
¯ = + O (h) ,
dx x=xn h
we obtain the recursion formula

yn+1 = yn + hf (xn , yn ) .

This can be used to calculate yn for arbitrary n, starting from a given y0 .

This is called the Euler method. Its local truncation error is ∼ h2 (one step), which
makes the global truncation error ∼ h (this is how the error scales for a fixed total
length in x).

A method whose global error scales as hn is called an nth -order method, so the Euler
method is a first-order method.

Low-order methods require small step sizes h, which makes the required number of
steps large. Therefore, the computational cost tends to be high.

6.3 Taylor Expansion

A straightforward way to reduce the truncation error is to make use of the Taylor
expansion ¯ ¯
dy ¯¯ h2 d2 y ¯¯
y (xn + h) = yn + h + + ...
dx ¯x=xn 2 dx2 ¯x=xn

50
where dy/dx = f and

d2 y df dy
2
= = ∂x f + ∂y f · = ∂x f + ∂ y f · f .
dx dx dx

Keeping the first three terms we obtain

h2 ¡ ¢
y (xn + h) = yn + hf (xn , yn ) + [∂x f (xn , yn ) + f (xn , yn ) ∂y f (xn , yn )] + O h3
2

which is a second-order method, or one order better than the Euler method. Clearly,
the same procedure can be used to obtain higher-order methods. This approach has,
however, the disadvantage that higher and higher derivatives of the function f will
be needed.

6.4 The Runge-Kutta Method

The Runge-Kutta method avoids this problem — no derivatives of f are needed.

Instead, the estimator yn+1 is constructed using values of f itself at carefully chosen
points between xn and xn+1 . By increasing the number of points, the order of the
method can be increased.

6.4.1 Second-Order Runge-Kutta (local error ∼ h3 )

Consider the ansatz

k1 = hf (xn , yn )
k2 = hf (xn + mh, yn + mk1 )
yn+1 = yn + ak1 + bk2

which combines two different f values and has three parameters: a, b and m.

The idea is to determine these three parameters so that the error becomes as small
as possible. For this purpose, we make a Taylor expansion of yn+1 about x = xn :

51
 ¡ ¢
k2 = hf (xn , yn ) + mh2 ∂x f (xn , yn ) + mh2 f (xn , yn ) ∂y f (xn , yn ) + O h3 ⇒
yn+1 = yn + (a + b) hf (xn , yn ) + bmh (∂x f (xn , yn ) + f (xn , yn ) ∂y f (xn , yn )) + O(h3 )
2

By comparing this expression with the Taylor expansion of the exact solution y(xn +h)
(see above), we see that the local error becomes O(h3 ) if

a+b=1 and bm = 1/2

This is the best that can be achieved with this ansatz — to eliminate the O(h3 ) error
term we would need more points.

The two equations for a, b and m have a one-parameter family of solutions. We can,
for example, take
a = 0, b = 1, m = 12
as in NR, or
a = b = 21 , m = 1
which is called Heun’s method.

6.4.2 Fourth-Order Runge-Kutta

Higher-order Runge-Kutta methods can be obtained by adding more points. This

gives improved accuracy at the cost of increased complexity. A very popular com-
promise is the fourth-order scheme given by Eq. 16.1.3 in NR, which requires four
f evaluations. One reason that this scheme is so popular is that in order to get a
fifth-order scheme, it turns out that six (and not five) points are needed.

52
 k1 = hf (xn , yn )
k2 = hf (xn + h/2, yn + k1 /2)
k3 = hf (xn + h/2, yn + k2 /2)
k4 = hf (xn + h, yn + k3 )
yn+1 = yn + k1 /6 + k2 /3 + k3 /3 + k4 /6 + O(h5 )

6.5 Adaptive Step Size

So far, we have assumed that the step size h is held constant throughout the integra-
tion. However, the function f (x, y) may be very different in different parts of space.
As a result, there may exist both “easy” regions where a large step can be used and
“difficult” ones where a much smaller step size is required.

Therefore, it is often very useful to vary the step size. A natural choice is to vary the
step size in such a way that the local error is kept at a constant level. For that, we
need to have an estimate of the local error.

A convenient way to estimate the local error is by step doubling. This means that
we repeat all calculations using a step size twice as large. By comparing the two
calculations, we get an estimate of the local error.

To see how this works, consider fourth-order Runge-Kutta, for which the local error
scales as h5 . Let y(x + h; h/2) and y(x + h; h) denote the estimates of the exact
solution y(x + h) that are obtained by using two steps of size h/2 and one step of
size h, respectively. For small h, these estimates should behave as
µ ¶5
h
y(x + h; h/2) − y(x + h) ∼ c(x) + O(h6 )
2
y(x + h; h) − y(x + h) ∼ h5 c(x) + O(h6 )

where c(x) is some unknown function. The quantity ∆ = |y(x + h; h) − y(x + h; h/2)|
provides a rough estimate of the error in y(x + h; h). To keep the local error under
control, we may require that ∆ < ∆tol , where ∆tol is a predetermined tolerance level.

If ∆ > ∆tol , we redo the calculation using a smaller step size h′ . How small the new
step size h′ should be can be estimated by using that ∆ ∼ h5 . This gives
µ ¶1/5
′ ∆tol
h ≈S h
∆

53
where h is the step size that gave us ∆. S is a “safety” factor that is supposed to be
less than 1.

If, on the other hand, ∆ ≪ ∆tol , the step size should be increased. How much it can
be increased can also be estimated by using that ∆ ∼ h5 .

The gain from using an adaptive step size can be huge compared to the cost of the
additional calculations that are needed to keep track of the local error.

6.6 The Modified Midpoint Method

Consider the same first-order ODE as before,

dy
= f (x, y) .
dx
A central-difference approximation of the derivative,
dy ¯¯ yn+1 − yn−1 ¡ ¢
¯ = + O h2 ,
dx x=xn 2h
gives us the so-called midpoint method,
yn+1 = yn−1 + 2hf (xn , yn ) ,
with a local error of O(h3 ). This is a “two-point” method, in which both yn and yn−1
are needed in order to obtain yn+1 .

(Compare this with the second order Runge-Kutta method with a = 0: yn+1 =
yn + hf (xn + h/2, yn + hf (xn , yn )/2) ”=” yn + hf (xn+1/2 , yn+1/2 ))

The modified midpoint method has three components:

1. An initial Euler step to get a second starting value:

y1 = y0 + hf (x0 , y0 )

2. N − 1 steps with the midpoint method:

yn+1 = yn−1 + 2hf (xn , yn ) n = 1, . . . , N − 1

3. A “correction” of the last value yN :

y(x0 + H; h) = 12 (yN + yN −1 + hf (xN , yN ))
where H = N h. y(x0 + H; h) is the final estimate, for step size h, of the exact
solution y(x0 + H).

54
For fixed x0 and H and even N , it has been shown (by Gragg) that

y (x0 + H; h) − y (x0 + H) = c1 h2 + c2 h4 + ...

In other words, the truncation error contains only even powers of h, a situation we
have encountered before (see Romberg’s integration method). This property makes
the method well suited for Richardson extrapolation.

The so-called Bulirch-Stoer method is based on the modified midpoint method and a
Richardson-like extrapolation (but with rational functions rather than polynomials).
This method can be a good choice if high accuracy is needed.

6.7 Predictor-Corrector Methods

All methods discussed so far have been explicit in that yn+1 could be calculated in
a direct manner, given one or more input values. There are also implicit methods,
in which an equation must be solved in order to obtain yn+1 . Implicit methods have
their advantages, as will be seen below in the discussion of stability. A simple example
of an implicit method is obtained by applying the trapezoidal rule for integration to
the equation dy/dx = f (x, y). This gives
xZ
n+1
h
yn+1 − yn = f (x, y) dx ≈ [f (xn , yn ) + f (xn+1 , yn+1 )] (i)
2
xn

This equation for yn+1 can be solved by using Newton’s method for root finding.
Another possibility is to use functional iteration. Consider the recursion formula
h
y (k+1) = g(y (k) ) g(y) ≡ yn + [f (xn , yn ) + f (xn+1 , y)] (ii)
2
The solution yn+1 to equation (i) is a fixed point of this formula; that is, yn+1 =
g(yn+1 ). If y (k) is close to yn+1 , we may linearize the recursion formula:
¡ ¢
y (k+1) = g yn+1 + y (k) − yn+1
¡ ¢h
≈ g (yn+1 ) + y (k) − yn+1 ∂y f (xn+1 , yn+1 )
2
¡ ¢ h
⇒ y (k+1) − yn+1 ≈ y (k) − yn+1 ∂y f (xn+1 , yn+1 )
2
This shows that if our initial guess y (0) is close enough to yn+1 for the linearization
to be OK, and if h is small enough that |h∂y f (xn+1 , yn+1 )/2| < 1, then |y (k) − yn+1 |
decreases with k. The desired value yn+1 can then be found by iteration.

55
Predictor-corrector methods use some recursion formula like (ii) to improve on an
initial guess of yn+1 . The recursion formula is applied a fixed number of times,
and not until some convergence criterion is fulfilled, which makes predictor-corrector
methods explicit rather than implicit.

A simple example of a predictor-corrector method, based on the trapezoidal rule, is:

( (0)
y = yn + hf (xn , yn ) “predictor” step
h
yn+1 = y (1) = yn + [f (xn , yn ) + f (xn+1 , y (0) )] “corrector” step
2
Note that this is equivalent to the second order Runge–Kutta with a = b = 1/2 and
m = 1.

A widely used, more advanced predictor-corrector method is the Adams-Bashforth-

Moulton method (see NR 16.7).

6.8 Stability

If n > 1, it may happen that different components yi evolve at very different speeds.
The problem is then called stiff. Stiff problems are generally hard to solve.

Example: Consider
dy
= −Ay y(0) = y0
dx
where y is a vector with two components and A is a 2 × 2 matrix. Assume that

Avi = λi vi (i = 1, 2) λ 2 ≫ λ1 > 0

and that y0 = c1 v1 + c2 v2 for some c1 , c2 . The exact solution is then given by

y(x) = c1 e−λ1 x v1 + c2 e−λ2 x v2

as can be easily verified. Here, the second term decays much faster than the first
one, because λ2 ≫ λ1 . As a result, the second term is not very relevant if we are
interested in “timescales” x & 1/λ1 .

Suppose we want to solve this problem numerically by using the Euler method. This
gives us the recursion formula

yn+1 = yn + h(−Ayn ) = (1 − hA)yn

56
(where 1 denotes the unit matrix) which can be easily solved:

yn = (1 − hA)n y0 =
= c1 (1 − hA)n v1 + c2 (1 − hA)n v2
= c1 (1 − hλ1 )n v1 + c2 (1 − hλ2 )n v2

If we take the limit h → 0 for fixed x = nh, we should recover the exact solution,
and this is indeed the case, because

(1 − hλi )n = (1 − λi x/n)n → e−λi x n→∞ (h → 0) .

Stability requires that |1−hλi | ≤ 1 for i = 1, 2 (otherwise two solutions corresponding

to slightly different initial values will diverge exponentially), and for this to hold it
is necessary to take h very small if λ2 is large,

1 − hλ2 ≥ −1 ⇒ h ≤ 2/λ2 .

Suppose now we want to study the behavior at large “timescales”, x & 1/λ1 . The
number of steps needed will then be very large, x/h & λ2 /2λ1 . Note that it is the
rapidly decaying, “uninteresting” component that dictates what step size we can use.

The Implicit Euler Method

If instead of a forward difference we use a backward difference to approximate the

derivative,
dy ¯¯ yn − yn−1
¯ ≈ ,
dx x=xn h
we obtain the so-called implicit Euler method. This may look like a minor modifi-
cation, but the stability properties change drastically. The recursion formula for the
implicit Euler method is

yn = yn−1 + h(−Ayn ) ⇒ yn = (1 + hA)−1 yn−1

which has the solution

yn = c1 (1 + hλ1 )−n v1 + c2 (1 + hλ2 )−n v2 .

Again, it is easy to verify that the exact solution is recovered as h → 0 for fixed
x = nh. The relative error in the second term will be large unless h is very small, but
this time this error is small in absolute numbers — |1/(1 + λi h)| ≤ 1 for all h > 0, so
the stability problems have disappeared. This is a major advantage compared to the
previous method. The disadvantage is that we have to solve an implicit equation for
yn for each n, which can be prohibitively time-consuming.

57
Example: Stiff problems are common in physics, especially for systems with many
degrees of freedom. As a not obvious example of a stiff problem, consider an ideal
ring polymer with potential energy
N
X
V = 1
2
(xn − xn+1 )2 (xN +1 = x1 )
i=1

and monomer units of mass 1.

The equations of motion for this system are (n = 1, . . . , N )

d2 xn
= −∇xn V =
dt2 £ ¤
= − 12 ∇xn (xn − xn−1 )2 + (xn − xn+1 )2 + {xn independent terms}
= xn−1 + xn+1 − 2xn

The general solution can be written as

N
X −1
£ ¡ ¢ ¡ ¢¤
xn (t) = A + Bt + Ck cos 2πkn
N
− ω(k)t + Dk sin 2πkn
N
− ω(k)t
k=1

where A, B, Ck and Dk are constant vectors (determined by the initial conditions)

and ω(k) = 2 sin πk/N . The fastest mode corresponds to k = N/2 (if N even) and
has a frequency of ωmax = 2. The k = 1 and k = N − 1 modes are the slowest ones,
with a frequency of ωmin = 2 sin π/N .

Now, if a method like the explicit Euler method is applied to this system, then the
step size must satisfy 1 & hωmax for the integration to be stable. This means that
a large number of steps will be needed in order to see a significant change in modes
with ω = ωmin , which are those with longest wavelength; the number of steps required
to integrate up to t = 1/ωmin is
1/ωmin ωmax 1 N
& = ∼ →∞ N →∞
h ωmin sin(π/N ) π

58
So, the number of steps needed to explore the long-wavelength or large-scale prop-
erties of the system increases with N . The computer time required for each step
increases, of course, too.

6.9 Molecular Dynamics (MD)

Classical MD (where quantum-mechanical and relativistic effects are ignored) amounts

to integrating Newton’s or Hamilton’s equations of motion for a system of particles,
N
X pi2
H= + V (x1 , . . . , xN )
i=1
2mi

dxi dpi
= pi /mi = −∇xi V (i = 1, . . . , N )
dt dt
MD simulations can be used to study both transport and equilibrium properties, and
are widely used not least in chemistry. The potential V may contain, for example,
Coulomb interactions and van der Waals terms.

Integrating the equations of motion above leaves, of course, the energy E = H

invariant (dH/dt = 0), as can be easily verified. Nevertheless, molecular dynamics
can be used for calculating equilibrium averages in the canonical ensemble, where the
energy is fluctuating. Such an average can be written as
Z
hOi = Oρcan d3N xd3N p

e−E/kT
ρcan = Z
e−E/kT d3N xd3N p

To be able to do this, it is clear that some mechanism must be added to enable

the system to make jumps in energy. One solution is to use a so-called Andersen
thermostat, which gives us the two-component algorithm:

1) Ordinary MD; that is, integration of the equation of motions.

2) “Stochastic collisions”. At randomly chosen or predetermined times t, the mo-

menta of the particles are refreshed, by drawing new values from the equilibrium
distribution:
¡ ¢
pi → p′i where P (p′i ) ∝ exp −(p′i )2 /2mi kT

59
Having added these collisions, canonical ensemble averages can be obtained as time
averages:
Z ZT
3N 3N 1
Oρcan d xd p = lim O (x (t) , p (t)) dt
T →∞ T
0

Assuming that the process is ergodic, this can be justified by showing that ρcan is
stationary. That ρcan is stationary under the stochastic collisions is obvious, since
the new p′i is drawn from this distribution. That ρcan is stationary under integration
of the equations of motion follows from Liouville’s theorem. This theorem says that
the flow in phase space (6N -dimensional) under the equations of motion is such that
the time evolution of a distribution ρ(x, p, t) must satisfy
N µ ¶
dρ ∂ρ X dxi dpi
= + ∇ xi ρ · + ∇pi ρ · = 0.
dt ∂t i=1
dt dt

So, ρ doesn’t change along the trajectories in phase space (the flow is incompress-
ible). If ρ = ρ(H(x, p), t), it can be easily shown that the sum in the expression for
dρ/dt vanishes, which implies that ∂ρ/∂t = 0. In particular, this shows that ρcan is
stationary, ∂ρcan /∂t = 0.

6.10 Numerical Integration of the Equations of Motion

Suppose we integrate the equations of motion above from time 0 to some time t. This
defines a mapping T (t) of phase space onto itself:

T (t) : R6N ∋ (x(0), p(0)) 7→ (x(t), p(t)) ∈ R6N

Two important properties of this mapping are:

60
 • By Liouville’s theorem, it conserves phase space volume — the Jacobian
¯ ¯
¯ ∂(x(t), p(t)) ¯
¯ ∂(x(0), p(0)) ¯ = 1
¯ ¯

• It is time reversible in the sense that

(xA , pA ) 7→ (xB , pB ) ⇒ (xB , −pB ) 7→ (xA , −pA )

When using a discretized version of the equation of motions, these properties are
generally lost — they hold approximately only. There are, however, discretization
methods such that these two properties remain exact.

6.10.1 The Leapfrog Method

This method is not very accurate (the local error is O(h3 )), but has the advantage
that it is time reversible and preserves phase space volume. It combines two different
elementary steps that represent Euler steps in x and p, respectively:
½
xi (t) 7→ xi (t + h) = xi (t) + hpi /mi
Ex (h) :
pi (t) 7→ pi (t + h) = pi (t)
½
xi (t) 7→ xi (t + h) = xi (t)
Ep (h) :
pi (t) 7→ pi (t + h) = pi (t) − h∇xi V (t)

for i = 1, ..., N .

Let S denote the mapping

½
xi → xi
(i = 1, . . . , N )
pi → −pi

An arbitrary mapping M is then time reversible (see above) if M −1 = SM S. It can

be easily verified that Ex (h) and Ep (h) both have this property, but

[Ex (h) Ep (h)]−1 = Ep (h)−1 Ex (h)−1

= SEp (h) Ex (h) S 6= SEx (h) Ep (h) S

The leapfrog method L(h) is obtained by taking a symmetric product,

L(h) = Ex (h/2)Ep (h)Ex (h/2)

61
(x and p can be interchanged). Thanks to the symmetry, L(h) becomes time re-
versible; that is,

L(h)−1 = SEx (h/2)Ep (h)Ex (h/2)S = SL(h)S

If we apply L(h) many times, the x and p variables will be updated following the
pattern in this figure:

L(h) is not only time reversible, but preserves phase space volume too; that is, the
Jacobian

¯ ¯
¯ ∂(x(h), p(h)) ¯
¯ ∂(x(0), p(0)) ¯ = 1 .
¯ ¯

This follows immediately from the time reversibility along with the fact that the
Jacobian for S 2 is 1 (because “the determinant of a product is the same as the
product of the determinants”). Alternatively, it can be shown that the Jacobians of
Ex (h) and Ep (h) both are 1. From this it follows that the Jacobian of L(h), which is
the product of three such factors, must be 1 as well.

62
7 Partial Differential Equations (PDE)

[NR 19.0, 19.1, 19.2, 19.5, 19.6 (19.3, 19.4)]

7.1 Introduction

Our discussion of PDEs will focus on three simple but important equations: the
diffusion or heat equation, the Poisson equation and the wave equation. The Poisson
equation differs from the other two in that there is no time dependence; this equation
is solved for given boundary conditions. The diffusion and wave equations involve
time and can be solved by evolving a given initial state in time, subject to given
boundary conditions.

The methods to be discussed are finite-difference methods, in which the system is

put on a space-time grid and the derivatives are approximated with finite differences.
This is a simple and useful approach, but there exist many other methods too.

An approach that is widely used, for example, in structural mechanics, is finite-

element methods. Here, the desired function u(x, t) is expanded in some basis func-
tions ϕi (x),
N
X
u(x, t) ≈ ai (t)ϕi (x)
i=1

The task then is to determine the expansion coefficients ai (t), by some suitable crite-
rion. Each basis function ϕi (x) is nonzero only in a local neighborhood (an element).

There are also variational methods. Many physical PDE problems can be reformu-
lated as variational problems; the desired function u(x, t) is an extremum of some
integral I,
Z
δI[u(x, t)] = δ dD xF [u(x, t)] = 0

(for an example, see hand-out for computer exercise 3). If we expand ϕi (x) as above,
this criterion gives a set of N equations for the expansion coefficients, ∂I/∂ai = 0,
which are called the Rayleigh-Ritz equations. If each basis function ϕi (x) is nonzero
only in a small element, this is a finite-element method, but this does not have to be
the case.

63
7.2 The Diffusion Equation

Consider a system of particles undergoing random collisions in d dimensions. Let

u(x, t) be the particle density at time t (x ∈ Rd ). The random motion of the particles
gives rise to a particle flow, or current, given by

j(x, t) = −D(x)∇u(x, t) ,

where D is the so-called diffusion coefficient. The dimensions of u and j are 1/lengthd
and 1/time·lengthd−1 , respectively, so D has dimension length2 /time. In equilibrium,
the net particle flow is zero, j = 0, and the distribution u is uniform.

From particle number conservation follows (via Gauss’ theorem) that the current
must obey a continuity equation,

∂u ∂u
+∇·j=0 ⇒ = ∇ · (D∇u) .
∂t ∂t

If there is also an external force F(x) = −∇U (x) acting on the particles, where U is
a potential, then the particles will have a drift velocity vdrift = µF, where µ is called
the mobility. The drift causes a particle flow, described by the current jdrift = uvdrift ,
and the total current now becomes j = −D∇u + jdrift . In equilibrium, the current
vanishes, which implies that
µu
∇u = F.
D
But we know that u is given by u ∝ e−U/kT in equilibrium. Combining these two
equilibrium equations, we obtain D = µkT , which is called the Einstein relation. It
follows that the time evolution of u in the presence of an external force F is given by
the equation · µ ¶¸
∂u 1
= ∇ · D ∇u + uF .
∂t kT
In our discussion, we will make the simplifying assumptions that there is no external
force, F = 0, and that the diffusion coefficient D is constant. The equation can then
be written as
à d
!
∂u X ∂ 2u
= D∆u ∆u = ∇ · ∇u =
∂t i=1
∂x2i

where ∆ is called the Laplace operator. We may furthermore assume that D = 1 (by
a suitable choice of units).

64
7.2.1 Example of an Exact Solution

It can be easily verified that

1 −x2 /4t
u(x, t) = √ e
4πt
3
solves the one-dimensional diffusion problem
∂u ∂2u
=
∂t ∂x2
u(x, t) → δ(x) tց0
This shows how a distribution initially localized at x = 0 evolves in time under
diffusive dynamics. The distribution is Gaussian with mean hxi = 0 for all t > 0,
and the variance increases linearly with t, hx2 i = 2t.

7.2.2 A Simple Finite-Difference Scheme

As a first example of a finite-difference scheme for the one-dimensional diffusion

equation, consider the following approximation:
u (xj , tn ) = unj
un+1
j − unj ∂u ∂ 2u unj+1 − 2unj + unj−1
≈ = ≈
h ∂t ∂x2 a2

3
R∞
The “delta function” δ(x) has the characteristic properties that −∞
δ(x)dx = 1 and that
δ(x) = 0 if x 6= 0.

65
This difference equation becomes particularly simple if the step sizes h and a are
chosen so that h/a2 = 1/2, which gives
¡ ¢
un+1
j = 1
2
unj+1 + unj−1 .

Schematically, this equation can be written as

Suppose we start from a state with one unj = 1 and all the others zero, to mimic the
exact calculation above. By iterating our discrete equation for the time evolution,
we obtain the solution

x→
0 0 0 0 1 0 0 0 0
0 0 0 12 0 21 0 0 0
t↓
0 0 41 0 24 0 41 0 0
0 18 0 38 0 83 0 18 0

which has similarities with the exact solution above.

7.2.3 A More Systematic Approach

Let us now solve the same one-dimensional problem in a more systematic manner.
To that end, we first discretize x only. Put

uj (t) = u(xj , t) (xj = ja)

2 ¯
∂ u¯ uj+1 − 2uj + uj−1
2
¯ ≈
∂x x=xj a2

The diffusion equation then becomes a system of first-order ODEs,

du
= Au ,
dt

66
where u(t) is a vector with components uj (t) and the matrix A is given by
 
 ... ... ... 
 0 

1  
1 −2 1 
A= 2 
a  1 −2 1 
 ... ... ... 
 0 

We have encountered this type of problem before. Three possible methods for solving
the system of equations are (un = u(tn ), tn = nh):

1. The explicit Euler method

un+1 − un
= Aun ⇒ un+1 = (1 + hA)un
h
2. The implicit Euler method
un+1 − un
= Aun+1 ⇒ un+1 = (1 − hA)−1 un
h
3. The trapezoidal rule
Z tn+1
n+1 n du h
u −u = dt ≈ (Aun + Aun+1 ) ⇒
tn dt 2
h h
un+1 = (1 − A)−1 (1 + A)un
2 2
This is called the Crank-Nicolson method.

The explicit and implicit Euler methods are first order in time, whereas the Crank-
Nicolson method is second order. We know, however, that the order of the method is
not the only issue in a problem like this; we must also look at the stability properties.

7.2.4 von Neumann Stability Analysis

For the continuous diffusion equation ∂u/∂t = ∂ 2 u/∂x2 , it is possible to find solutions
of the form u(x, t) = T (t)X(x). These are given by u = e−ω(k)t eikx , where ω(k) = k 2 .
For a discretized version of the equation, we can make a similar ansatz,
un = ξ(k)n ũ(k) ,

67
where ξ(k) is a number corresponding to e−ω(k)h , and ũ(k) is a vector with components
(k)
ũj = eikxj , corresponding to the function eikx . It is not unreasonable to expect ũ(k) to
be an eigenvector of the matrix A, because eikx is an eigenfunction of the “operator”
∂ 2 /∂x2 (with eigenvalue −k 2 ), and A is our discretized version of this operator. A
direct calculation shows that this is indeed the case:
1 ikxj−1
(Aũ(k) )j = (e − 2eikxj + eikxj+1 )
a2
1 ¡ ¢ (k)
= − 2 2 − e−ika − eika ũj
| a {z }
λk

So,
4 ka
Aũ(k) = λk ũ(k) λk = − 2
sin2
a 2
Note that λk ≈ −k 2 if ka ≪ 1.

Assuming that the initial state can be written as a superposition of solutions of the
form un = ξ(k)n ũ(k) , we can now draw conclusions about the stability of the three
methods mentioned above.

1. Explicit Euler

ξ(k)n+1 ũ(k) = ξ(k)n (1 + hA)ũ(k) ⇒ ξ(k) = 1 + hλk

Stability requires that |ξ(k)| 6 1 for all k, so we must have (note that λk 6 0)

2h 2 ka
1 + hλk > −1 ⇔ sin 61 (all k)
a2 2
This holds for all k only if h 6 a2 /2. This is a severe restriction.
√ For the exact
solution, we saw that the width of the distribution scaled as time. This means
that the number of steps needed before the distribution reaches a given width
L scales as
no. of steps = time/h ∼ L2 /h & (L/a)2

2. Implicit Euler
1
ξ(k)n+1 ũ(k) = ξ(k)n (1 − hA)−1 ũ(k) ⇒ ξ(k) =
1 − hλk

Since λk 6 0, it immediately follows that 0 6 ξ(k) 6 1 for all k and h > 0.

Hence, the method is stable for any h > 0.

68
 3. Crank-Nicolson

h h 1 + h2 λk
ξ(k)n+1 ũ(k) = ξ(k)n (1 − A)−1 (1 + A)ũ(k) ⇒ ξ(k) =
2 2 1 − h2 λk

This implies that |ξ(k)| 6 1 for all k and h > 0. This follows from the facts
that all λk 6 0 and that |(1 − x)/(1 + x)| 6 1 for all x > 0. So, this method is
also stable for all h > 0.

7.2.5 Implementing the Implicit Euler and Crank-Nicolson Methods

At each implicit Euler step, we need to solve (1 − hA)un+1 = un for un+1 . For
a general A this would be very time-consuming. However, A is tridiagonal, which
makes it doable.

Assume, for simplicity, that we are studying an interval 0 6 x 6 Ja and that

the boundary conditions are u(0, t) = C1 and u(Ja, t) = C2 , where C1 and C2 are
constants. Let u be the vector with the J − 1 components u1 , . . . , uJ−1 (u0 and uJ
are left out because they are constant). The problem can then be written
du
= Au + c
dt
where  
−2 1  
C1
 1 −2 1 0 
   0 

 1 −2 1 


 ..


A= ... ... ...  c= . 
   
   0 
 0 1 −2 1 
C2
1 −2

69
The (J − 1) × (J − 1)-matrix A is called tridiagonal because it has nonzero elements
only in three bands along its diagonal. The matrix 1−hA is tridiagonal too, and this
makes it possible to solve the equation (1 − hA)x = y for x in O(J) operations, by
Gauss elimination. This figure illustrates how this is done for an arbitrary tridiagonal
N × N -matrix:

 
b1 c 1
 a2 b 2 c 2 
  −a2 /b1

 a3 b 3 c 3 
 ←֓ →

 ... ... ... 


 
b1 c1
 0 b′ = b2 − a2 c1 /b1 c2 
 2 

 a3 b3 c 3 
 −a3 /b′2 → ... →

 . .. ... ... 
 ←֓

 
b1 c 1

 b′2 c2 


 b′3 c3 

 ... ...  → ...
 
 
 b′N −1 cN −1  "
b′N −cN −1 /b′N

There are ≈ N steps “downwards” and ≈ N steps “upwards” and each step requires
O(1) operations, which makes the total cost O(N ).

In the Crank-Nicolson method, there is also a tridiagonal system of equations to be

solved at each step. The same method can be used in this case.

7.2.6 Higher Dimensions

So far, we have restricted ourselves to one-dimensional diffusion (d = 1). Formally, it

is straightforward to extend these methods to higher d. To illustrate that, consider
the two-dimensional diffusion problem

∂u ∂ 2u ∂ 2u
= +
∂t ∂x2 ∂y 2

70
with fixed boundary conditions, u = 0. Assume, for simplicity, that the region studied
is rectangular. We can then discretize x and y as follows.

xj = ja (j = 0, . . . , J)
yl = la (l = 0, . . . , L)
ujl (t) = u(xj , yl , t)
dujl uj+1l + uj−1l + ujl+1 + ujl−1 − 4ujl
≈
dt a2
Put all “internal” (non-constant) ujl in a vector with a single index i, which may be
defined as i = (j − 1)(L − 1) + l − 1.

We can then write the discretized diffusion equation as

du
= Au (+ constant vector if there are nonzero boundary values)
dt
The Crank-Nicholson method, for example, can be directly applied to this equation,
which gives
un+1 = (1 − h2 A)−1 (1 + h2 A)un .
Formally, this looks precisely as before. However, the matrix A now has 5 bands
instead of 3, which makes the problem computationally harder. The cost of the
calculation of un+1 scales linearly with the number of rows or columns of the matrix
in this case too, but this number is larger, (J − 1) · (L − 1), and the prefactor in this
scaling law is higher when the number of bands is 5 instead of 3.

The Alternating Direction Implicit method (ADI)

ADI is a method that, like the Crank-Nicolson method, is second order in time and
stable for all h > 0. The advantage of ADI compared to Crank-Nicolson is that it

71
works with tridiagonal matrices. This can be achieved by splitting the matrix A into
two terms,
A = Ax + Ay ,
where Ax and Ay correspond to ∂ 2 /∂x2 and ∂ 2 /∂y 2 , respectively. Both Ax and Ay
are tridiagonal. To be able to take advantage of that, each ADI step is split into two
substeps:
½ n+1/2
u = un + h2 (Ax un+1/2 + Ay un ) “implicit in x, explicit in y ′′
un+1 = un+1/2 + h2 (Ax un+1/2 + Ay un+1 ) “explicit in x, implicit in y ′′
There are two equation systems to be solved, one for each substep, and they are both
tridiagonal,
(1 − h2 Ax )un+1/2 = (1 + h2 Ay )un
(1 − h2 Ay )un+1 = (1 + h2 Ax )un+1/2
The method can be summarized in one equation by eliminating un+1/2 , which gives
un+1 = (1 − h2 Ay )−1 (1 + h2 Ax )(1 − h2 Ax )−1 (1 + h2 Ay )un

7.3 The Poisson Equation

Consider the Poisson equation

∆u(x) = ∇ · ∇u(x) = −ρ(x) x∈Ω
with the boundary condition u(x) = f (x), x ∈ ∂Ω (∂Ω denotes the boundary of Ω).
Such a condition that specifies u itself on the boundary (f is a given function) is
called a Dirichlet boundary condition.

This boundary-value problem appears at different places in physics. An example

is electrostatics. u is then the electrostatic potential and ρ the charge density (in
suitable units).

Suppose we want to solve this problem by finite differencing. The problem is then
transformed into a linear algebra problem of the form Au = b, where A, as before,
corresponds to ∆ and b is a vector that is determined by ρ and the function f in the
boundary condition.

Solving this system of equations is, in principle, straightforward. The problem is the
size of the system; if, for example, we are in d = 2 and work with a 1000 × 1000 grid,
then there are 106 u components and 1012 elements in the matrix A. This makes it
necessary to make use of the fact that A is sparse (that is, only a tiny fraction of its
elements are nonzero).

72
7.3.1 Relaxation Methods

A widely used strategy for solving the boundary value problem above is to let the
system evolve in a fictitious time τ with diffusive dynamics. We then consider the
diffusion problem

∂u

 = ∆u + ρ x ∈ Ω, τ > 0
∂τ
 u=f x ∈ ∂Ω, τ > 0

some initial condition τ = 0

where ρ and f are the same as in the original problem. In general, u will approach
a stationary state as τ → ∞; that is, ∂u/∂τ → 0 as τ → ∞. This stationary state is
the desired solution.

Suppose we evolve u in τ by using a finite-difference method of the form

un+1 = Tun + constant vector .

The matrix T is called the iteration matrix. The convergence rate is governed by the
spectral radius ρT of this matrix, defined by

ρT = max |λ| .
eigenvalues λ of T

To see this, we first note that the desired solution must correspond to a fix point of
the recursion formula. Denote this fix point by uf . The deviation from the fix point
satisfies
un+1 − uf = T(un − uf ) .
which implies that
|un − uf | ∼ ρnT n→∞
(unless the projection of u0 − uf onto the eigenvector with maximum |λ| happens to
be zero). This first of all shows that the method will converge if ρT < 1. We also see
that for the convergence to be fast, ρT should be as small as possible.

7.3.2 The Jacobi method

A simple relaxation method is the Jacobi method, which is obtained by using an

Euler step in τ . The method is given by

un+1 = Tun + ρ ,

73
where the iteration matrix T = 1 + hA and the vector ρ represents the charge
density ρ(x) (for simplicity, we assume that f = 0). As usual, the matrix A is our
approximation of the Laplace operator ∆. It is the same as before.

Assume that we are in two dimensions and that the region Ω is a square,

Ω = {(x, y); 0 6 x, y 6 Ja}

It is then easy to find the eigenvalues and eigenvectors of A, which are given by

Av(kx ,ky ) = λkx ,ky v(kx ,ky ) kx , ky = 1, . . . , J − 1

µ ¶
4 2 kx π 2 ky π
λkx ,ky = − 2 sin + sin
a 2J 2J
(k ,k ) jkx π lky π
vjl x y = sin sin
J J
The eigenvalues of the iteration matrix T are given by

λTkx ,ky = 1 + hλkx ,ky .

It is clear that λTkx ,ky 6 1, since λkx ,ky 6 0. The lowest eigenvalue of T is

4h
λTJ−1,J−1 = 1 + hλJ−1,J−1 ≈ 1 − 2 .
a2
This shows that in order to have ρT 6 1, we must take h 6 a2 /4. In the Jacobi
method, one takes h = a2 /4. With this choice, the time evolution is given by
1 n a2
un+1
jl
n n n
= (uj+1l + uj−1l + ujl+1 + ujl−1 ) + ρjl .
4 4
To find ρT and thereby the convergence rate, we note that
π π2
max λTkx ,ky = 1 + hλ1,1 = 1 − 2 sin2 ≈1− 2
2J 2J
(J − 1)π (J − 1)π π2
min λTkx ,ky
= 1 + hλJ−1,J−1 = 1 − 2 sin 2
= cos ≈ −1 + 2
2J J 2J
2 2
So, ρT ≈ 1 − π /2J . The number of iterations, r, required to reduce the error by a
factor 10−p satisfies
−p ln 10 2p ln 10 2
ρrT ∼ 10−p ⇒ r≈ ≈ J .
ln ρT π2
Hence, the number of iterations scales quadratically with J. This makes the method
too slow to be of direct practical interest.

74
7.3.3 The Gauss-Seidel Method

This method is very similar to the Jacobi method. The only difference is that “new”
u values are used as soon as they become available. This gives a slightly better
convergence rate, but the scaling of the required number of iterations with system
size remains quadratic.

The Gauss-Seidel method can be written as

1 n a2
un+1 = (u + u n+1
+u n
+ u n+1
) + ρjl
jl
4 j+1l j−1l
↑
jl+1 jl−1
↑
4

The Jacobi and Gauss-Seidel method are slow but important as starting points for
more advanced methods.

Successive overrelaxation (SOR; see computer exercise 3) is based on the Gauss-Seidel

method. It is almost as simple as the Gauss-Seidel method to implement, but much
more efficient — the number of iterations needed scales as J rather than J 2 , provided
that the overrelaxation parameter ω is carefully chosen.

Multigrid methods (see NR) are more advanced. They are very efficient, but also
more complicated than SOR to implement.

7.3.4 The Fourier Transform Method

Consider the boundary value problem

 2
 ∂ u ∂2u

 + = −ρ 0 < x < X, 0 < y < Y
∂x2 ∂y 2

 u(0, y) = u(X, y), u(x, 0) = u(x, Y ) periodic boundary conditions

u(0, 0) = u0

75
A very efficient method for this problem can be obtained by using Fourier decompo-
sition:
∞ ∞ k x ky y
„ «
X X 2πi x +
X Y
u(x, y) = û(kx , ky )e
kx =−∞ ky =−∞
k x ky y
Z Z „ «
X Y
1 −2πi x +
X Y
û (kx , ky ) = dx dyu(x, y)e
XY 0 0

Do the same for ρ(x, y). Insertion into the Poisson equation gives
h¡ ¢ ³ ´2 i
kx 2
−4π 2
X
+ kYy û (kx , ky ) = −ρ̂ (kx , ky ) ,

which is a simple algebraic equation for û, in place of the original differential equation.

This gives us a three-step “algorithm” for solving the given problem:

1. Transform ρ → ρ̂.
2. Calculate û(kx , ky ) by using the algebraic equation above for all (kx , ky ) 6= (0, 0).
For kx = ky = 0, û is determined by the additional condition u(0, 0) = u0 .
3. Transform û → u.

Numerically, these calculations are done using discrete Fourier transforms,

N N
1 X ˆ 2πikn 1 X 2πikn
fn = √ fk e N fˆk = √ fn e− N
N k=1 N n=1

With the method of Fast Fourier Transform (FFT), the cost of a transform fn 7→ fˆk
(or vice versa) becomes ∼ N ln N rather than N 2 .

This makes this method very fast. However, it requires that the region considered is
rectangular, which is a strong limitation (although the boundary conditions do not
have to be periodic).

7.4 Waves

The one-dimensional wave equation

∂ 2u 2
2∂ u
− v =0
∂t2 ∂x2

76
describes, for example, the transverse motion of a free string. This can be seen by
considering a small string element dx. The string is assumed to be homogeneous
with a linear density ρ, so the mass of the small element is ρdx. Furthermore, it is
assumed that the tension T is constant, and that the motion is “small” in the sense
that the angle θ in the figure below stays small. The equation for the transverse
motion of the small element then is
∂ 2u
ρdx = T sin θ|x+dx − T sin θ|x ≈ T tan θ|x+dx − T tan θ|x ⇒
∂t2
∂ 2u ∂u ¯¯ ∂u ¯¯ ∂ 2u ∂2u T ∂2u p
ρdx 2 ≈ T ¯ − T ¯ ≈ T dx 2 ⇒ ≈ (v = T /ρ)
∂t ∂x x+dx ∂x x ∂x ∂t2 ρ ∂x2

The one-dimensional wave equation can be solved by changing variables to ξ = x − vt

and η = x + vt. The equation then becomes

∂ 2u
=0 ⇒ u = f (ξ) + g(η) = f (x − vt) + g(x + vt) ,
∂ξ∂η

where the functions f and g are determined by the initial values of u itself and the
velocity ∂u/∂t (the equation is second order in time). f (x−vt) and g(x+vt) represent
right- and left-moving waves, respectively (v > 0).

77
Instead of the wave equation, we will, for simplicity, consider
∂u ∂u
+v = 0,
∂t ∂x
which has right-moving solutions only; the same change of variables as before gives
∂u
=0 ⇒ u = f (ξ) = f (x − vt) .
∂η

It is instructive to solve this problem by separation of variables, which means that

we write the solution as a superposition of solutions of the form u(x, y) = X(x)T (t).
This ansatz gives
X ′ (x) 1 T ′ (t)
=− ≡ ik = constant ⇒ X(x) ∝ eikx , T (t) ∝ e−ikvt ,
X(x) v T (t)
where the constant k must be real for X to remain bounded as x → ±∞. This shows
that uk (x, t) = ei(kx−ω(k)t) , where ω(k) = kv, is a solution to the given equation for
all real k. The general solution can be written as a continuous superposition of such
solutions, Z ∞
1
u(x, t) = √ c(k)ei(kx−ω(k)t) dk ,
2π −∞
where c(k) is determined by initial data. Suppose a (right-moving) “wave packet” is
expressed in this way. The fact that ω(k) = kv implies that all plane-wave compo-
nents uk of the wave packet have the same phase velocity (v). This in turn implies
that the wave packet propagates without changing shape.

If the phase velocity vf = ω(k)/k were k dependent, then the shape of a wave packet
would change with time. This is called dispersion.

This may be further generalized by considering a complex ω(k) = Ω(k)−iγ(k), where

Ω and γ both are real. In this case, it can be shown (by using Parseval’s formula)
that
Z ∞ Z ∞ ¯Z ∞ ¯2
2 1 ¯ i(kx−ω(k)t)
¯
|u(x, t)| dx = ¯ c(k)e dk ¯ dx
−∞ 2π ¯ ¯
Z−∞
∞
−∞
1 ′ ′
= dx dk dk ′ c(k)c̄(k ′ )ei(kx−ω(k)t) e−i(k x−ω̄(k )t)
2π −∞
Z ∞ Z ∞
1 ′ ′ i(ω̄(k′ )−ω(k))t) ′
= dk dk c(k)c̄(k )e dxeix(k−k )
2π
Z ∞ −∞ −∞

= |c(k)|2 e−2γ(k)t dk .
−∞

78
 R
From this result we see that the normalization |u|2 dx is preserved as long as ω is
real (γ = 0). If, on the other hand, γ > 0, the normalization integral decreases with
time. This is called damping.

7.4.1 Discrete waves

Let us now look at how to solve the equation ∂u/∂t = −v∂u/∂x by finite differencing.
As usual, we put unj = u(xj , tn ), where xj = ja and tn = nh. For the space derivative,
we use a central difference,
∂u ¯¯ unj+1 − unj−1
¯ = + O(a2 )
∂x t=tn 2a

The Euler Method

With an Euler step in time, we get

un+1
j − unj unj+1 − unj−1 hv
= −v ⇒ un+1
j = unj − (unj+1 − unj−1 )
h 2a 2a
To check the stability of the method, we make a von Neumann ansatz, unj = ξ(k)n eikxj .
This gives
hv hv
ξ(k) = 1 − (eika − e−ika ) = 1 − i sin ka ⇒
2a a
µ ¶2
2 hv
|ξ(k)| = 1 + sin2 ka > 1
a
showing that |unj | = |ξ(k)|n grows exponentially with n for all h > 0. So, this method,
which may not seem unreasonable, is a disaster — it is always unstable.

The Lax Method

The method above can be improved by a simple modification, which gives the so-
called Lax method. The modification is that unj is replaced by the average of unj+1
and unj−1 . This gives
1 hv
un+1
j = (unj+1 + unj−1 ) − (unj+1 − unj−1 ) .
2 2a
n n ikxj
This time the ansatz uj = ξ(k) e gives
1 hv hv
ξ(k) = (eika + e−ika ) − (eika − e−ika ) = cos ka − i sin ka ,
2 2a a
and from the figure below it can be seen that |ξ(k)| 6 1 for all k if hv/a 6 1. So, the
method is stable if h 6 a/v. How does the discrete solution behave?

79
The k dependence of ξ(k) corresponds to an ω(k) given by

hv
ξ(k) = e−iω(k)h = cos ka − i sin ka .
a
If hv/a = 1, we see that hω = ka. So, ω = ka/h = kv, which happens to be the
exact result.

If, on the other hand, hv/a < 1, then

µ ¶2
2 −2hγ 2 hv
|ξ| = e = cos ka + sin2 ka < 1 .
a

This means that γ > 0, so the discrete solution is, in contrast to the exact one,
damped. The damping is, however, small if ka is small. The real part Ω of ω satisfies
µ ¶
(hv/a) sin ka
hΩ = arctan ≈ arctan hvk ≈ hvk
cos ka

if ka is small. This shows that ω(k) for the discrete solution is approximately correct
if the wavelength is much larger than a, so that ka ≪ 1. Short-wavelength modes
have large relative errors but are damped, so these errors do not destroy the long-
wavelength properties of the solution.

The Leapfrog Method

The Lax method is first order in time and second order in space. A method that is
second order in time too can be obtained by using a central-difference approximation
for the time derivative, which gives

vh ¡ n ¢
un+1
j = un−1
j − uj+1 − unj−1 .
a

80
For this method, the ansatz unj = ξ(k)n eikxj gives
s µ ¶2
vh vh
ξ(k) = −i sin ka ± 1 − sin2 ka .
a a

From this it can be easily verified that this method, like the Lax method, is stable if
hv
6 1.
a
This is called the leapfrog method because the space-time grid separates into two
sub-lattices that do not influence each other. This may cause numerical problems;
the sublattices may “diverge” due to rounding errors. An artificial coupling term is
then needed.

7.5 Schrödinger’s equation

Consider the one-dimensional Schrödinger equation

∂ψ ∂2ψ
i = − 2 + V (x) ψ = Hψ
∂t ∂x
(in units such that ~ = 1, m = 1/2). For V = 0, this equation describes diffusion in
imaginary time; in fact, with τ = it and V = 0, we obtain ∂ψ/∂τ = ∂ 2 ψ/∂x2 .

After discretizing x, the Schrödinger equation becomes

dψ
i = HD ψ
dt
where ψ is a vector with components ψj (t) = ψ(xj , t) and HD is a discrete version
of H that is assumed to be Hermitian, H†D = HD .

A fundamental property of the exact solution is that the normalization is preserved,

Z ∞
|ψ(x, t)|2 dx = 1 .
−∞

81
Does this remains true after discretizing x? The discrete analog is
µ ¶†
d † 1 1
(ψ ψ) = HD ψ ψ + ψ † HD ψ = 0 .
dt i i
The eigenvalues EI of HD are real, since HD is Hermitian. Denote the eigenvectors
by ψ I ,
HD ψ I = EI ψ I .
The eigenvectors may be assumed to be orthonormal.

Let us now consider three possible ways of discretizing t.

1. The explicit Euler method gives

ψ n+1 = ψ n − ihHD ψ n = Tψ n ,
where the iteration matrix T = 1 − ihHD . The eigenvalues of T are λI =
1 − ihEI and satisfies |λI | > 1 for all I. Suppose now
X X
ψ0 = cI ψ I ⇒ ψn = cI λnI ψ I .
I I

From the fact that the eigenvectors are orthonormal, it follows that the nor-
malization X
|ψ n |2 = ψ n† ψ n = |cI |2 λ2n
I .
I
n 2
This implies that |ψ | is not constant, but increases exponentially with n (for
large n).
2. With an implicit Euler step instead, we obtain
ψ n+1 = ψ n − ihHD ψ n+1 ⇒ ψ n+1 = Tψ n
with T = (1 + ihHD )−1 . The eigenvalues of T are this time given by λI =
1/(1 + ihEI ), so |λI | 6 1 for all I. It follows that |ψ n |2 decreases steadily with
n.
3. Finally, we consider the Crank-Nicholson method (or trapezoidal rule),
h
ψ n+1 = ψ n + (−iHD ψ n − iHD ψ n+1 ) ,
2
which can be written ψ n+1 = Tψ n with
µ ¶−1 µ ¶
h h
T = 1 + i HD 1 − i HD .
2 2

82
 The eigenvalues of this T are

h
1 − i EI
λI = 2 .
h
1 + i EI
2
Since λI is a ratio between two numbers that are the complex conjugates of each
other, it follows that |λI | = 1 (for all I), and therefore that |ψ n |2 is constant.
So, the Crank-Nicolson method has the advantage that the normalization is
preserved.

8 Appendix: χ2 Minimization

[NR 15.5]

When discussing minimization, we mentioned the Newton method for root finding.
This method can be used to solve the equation system ∇f (x) = 0, and thereby find
the minimum of a given cost function f . Usually, this is an impractical method, but
it is useful for function fitting.

Suppose we want to fit a set of data points (xi , yi ) with errors σi , i = 1, . . . , N , to a

function y(x; a), where a = (a1 , . . . , aM ) are parameters (M < N ). This is usually
done by minimizing
N µ ¶2
2
X yi − y(xi , a)
χ (a) =
i=1
σi
with respect to a.

If the function y depends linearly on the parameters ai , this amounts to solving

a relatively simple linear algebra problem. To see this, assume that y is a linear
combination of some basis functions ϕk ,
M
X
y(x; a) = ak ϕk (x) .
k=1

χ2 can then be expressed as

N
à PM !2 N
à M
!2
X yi − k=1 ak ϕk (xi )
X X
χ2 (a) = = bi − Aik ak = |b − Aa|2 ,
i=1
σi i=1 k=1

83
where A is the matrix with elements Aik = ϕk (xi )/σi and b is the vector with
components bi = yi /σi . Putting ∇χ2 (a) = 0, we get a linear equation system for the
parameters ai ,
AT Aa = AT b .
These equations are called the normal equations and can be solved, for example, by
singular value decomposition (SVD; see NR 2.6 and p. 676).

If y is a non-linear function of the parameters ai , it is necessary to take a different

approach. A popular choice is to use the Levenberg-Marquardt method, which can
be thought of as a combination of steepest descent and the Newton method. In the
Levenberg-Marquardt method, we need the gradient of χ2 ,
N
∂χ2 X yi − y(xi , a) ∂y(xi , a)
= −2
∂ak i=1
σi2 ∂ak

as well as the Hessian,

N · ¸
∂ 2 χ2 X 1 ∂y(xi , a) ∂y(xi , a) ∂ 2 y(xi , a)
=2 − (yi − y(xi , a)) .
∂ak ∂al σ2
i=1 i
∂ak ∂al ∂ak ∂al

Here the last term should be small if we are near the minimum and the fit is good
(it vanishes if y depends linearly on the parameters). For convenience, this term is
neglected, which does not affect the final ak values and which gives
N
∂ 2 χ2 X 1 ∂y(xi , a) ∂y(xi , a)
≈ Akl = 2 .
∂ak ∂al σ 2 ∂ak
i=1 i
∂al

Put
1 ∂χ2 1
βk = − αkl = Akl δak = a′k − ak .
2 ∂ak 2
In this notation, a steepest descent step has the form δak = constant · βk . The
Levenberg-Marquardt method does not use a proper steepest descent step but rather
1
δak = βk ,
λαkk

where λ is a parameter (see below) and αkk is introduced because the different com-
ponents ak may behave very differently; dividing by αkk makes sense dimensionally.
Note that αkk is guaranteed to be positive because we are using the simplified version
of the Hessian.

84
A Newton step can be written as
X
a′ = a − A(a)−1 ∇χ2 (a) ⇒ αkl δal = βk .
l

The Levenberg-Marquardt method can be seen as an elegant way to interpolate be-

tween these two types of step, by changing the parameter λ. This is achieved by
considering
½
X αkk (1 + λ) k = l
α̃kl δal = βk where α̃kl =
αkl k 6= l
kl

For λ → 0, this becomes a Newton step, because the matrices αkl and α̃kl are the
same in this limit. For large λ, on the other hand, we get a modified steepest descent
step, because the diagonal elements dominate the matrix α̃kl .

The idea is to start with a large λ which gives a steepest descent-like step, and then
gradually decrease λ as we get closer to the minimum, so that the step becomes more
and more Newton-like.

Schematically, the algorithm can be written as:

1. Pick some initial a and a large λ.

2. Solve the equation system for δa and calculate χ2 (a + δa).

3. If χ2 (a + δa) < χ2 (a), update a to anew = a + δa, decrease λ, and continue at

2. If χ2 (a + δa) ≥ χ2 (a), increase λ (decrease the step size) and go back to 2
without changing a.

Iterate till some stopping criterion is fulfilled.

85