PHY400W — cp 2003 1
A On Computational Physics
When a problem is clearly stated, it is of no further interest to
— P. Debye
Analysing a general potential system with two degrees of freedom
is beyond the capabilities of modern science
— V.I. Arnol’d
Physics is the study of how the world works.
From experiment we derive some sort of intuition about how the world op-
erates, and try to distill this into some kind of physical law or model. From
then on, it’s just a matter of expressing this in mathematical terms, and
exploring the quantitative outcome of our intuition.
In a typical physics course, we learn how to express the ideas of physics
as, say, diﬀerential equations, which we can then solve. Indeed, much of
the physics course is devoted to all those lovely equations and their neat
closed-form (“analytical”) solutions: harmonic oscillators, point masses in
gravitational ﬁelds, hydrogen atoms, and so on. Such is the power of the
diﬀerential and integral calculus. It’s very easy to gain the impression that
all of the equations of physics is of this nature, and that it is just a measure
of our lack of mathematical sophistication that we can’t ﬁnd an analytical
solution to some problem.
Unfortunately, this is not true. Most problems do not have nice analytical
solutions, and there is no hope of ﬁnding one. So we have to turn to numerical
solutions in order to solve our problems.
Arnol’d quotes the set of solvable problems in classical mechanics:
one-dimensional problems, motion of a point in a central ﬁeld,
Eulerian and Lagrangian motions of a rigid body, the problem of
two ﬁxed centers, and motion along geodesics on the ellipsoid.
PHY400W — cp 2003 2
To this list we can add the instructive case of the Toda lattice. This was only
found to be solvable (integrable) after numerical solutions had hinted at it.
But let’s consider one of these integrable cases, the pendulum. This has the
equation of motion
= − sin θ
This is solvable only in that the solution can be expressed in terms of known
functions, the elliptic integrals. To evaluate these functions is still a problem.
But let’s take a step back: to evaluate sin θ is a problem. We don’t consider
it so these days, we just punch a few keys on a pocket calculator. (A few
decades ago we would have looked up the value in a set of tables, calculated
by a human). But our calculator (whether electronic or human) only does
arithmetic. In order to calculate the sine it must go through some arithmetic
procedure, a step by step recipe, to calculate the value required using arith-
metical operations. Developing such a step by step recipe, an algorithm, is
the basis of ﬁnding our solutions in computational physics. We might think
this is rather simple. After all, we can write a series expansion for the value
of our sine. But we want a solution in a ﬁnite time, so we would need to
truncate the series. This will incur a truncation error, so we will need to
make sure that it is small enough. Then we ﬁnd that there are numerical
techniques to enhance the accuracy with fewer series terms, by slight mod-
iﬁcation of the coeﬃcients. In fact, the algorithm used in our calculator is
quite far from the series that you might expect.
A.2 Computers and physics
The development of computers has been driven by the needs of physics.
By this I mean not just the collection and analysis of data, but the compu-
tation of mathematical results that cannot be obtained by analytical means:
in other words, number crunching. The most powerful and advanced systems
of the day have always been designed for this.
This means that programming for physics is a much more complicated task
than, say programming an operating system. We do have to worry about
eﬃciency and program size; and not only must the coding be correct, but
also the mathematical representation of the problem.
PHY400W — cp 2003 3
A.3 What is computational physics?
There are two major aspects to computation in physics:
• Obtaining solutions to known equations which cannot be solved ana-
• Exploring model systems for which no solution is known.
These both exploit algorithms for numerical methods of solution on comput-
To apply these techniques we have to understand their limitations, e.g the
accuracy of approximations, the stability of results against numerical per-
turbations, and the limitations of the representation of real numbers on the
The algorithms then have to be translated into a list of instructions for
the computer. This usually involves a computer language in which these
instructions are given.
PHY400W — cp 2003 4
I don’t know what the language of the year 2000 will look like,
but I know it will be called Fortran.
– C.A.R. Hoare
B.1 Choice of language
In order to represent a problem on the computer we require a programming
language which we can use to encode the processes we want to carry out.
For instance, we need to be able to tell the computer “add this number to
that variable” or “print out this value”.
Such languages are of course formal languages that have a precise grammar
and syntax. We can construct a program which is a list of sentences that
accord with these rules of grammar and syntax. The elements of this language
have some semantic meaning in terms of the numerical tasks we require to
The program has to be run or executed on the computer. There are two basic
ways of achieving this. The program can be interpreted by another program,
which carries out the tasks speciﬁed in the program. This process can be
rather slow. A more eﬃcient scheme is to compile the program into machine
code (the code understood by the CPU). This machine code then directly
controls the hardware. (Typically, a program is compiled into object code,
which is then linked to other object ﬁles to form an executable ﬁle. Some of
these other object ﬁles may be in the form of libraries which hold commonly
The language of computational physics has traditionally been Fortran, but
these days C and C++ are being used to increasing extents, partly due
to the availability of compilers rather than any real beneﬁt. Indeed, array
operations are rather problematic in C-based languages due to the internal
representations, and C++ is only recently approaching the ease of use and
speed of Fortran, with the use of exotic template techniques that can defeat
In this course we will take a slightly easier route and give examples in Python.
Why use a barely known and obscure language like Python?
PHY400W — cp 2003 5
• Python is free and freely available for many platforms.
• Python is a modern language, with object oriented features. It is much
easier to learn than C++.
• Python is interpreted; this makes it easy to ‘see what happens if we. . . ’
• Python is easily extensible and has a large variety of available exten-
sions and libraries.
• With the Numerical Python extension, Python has many of the array-
handling features of Fortran 90/95.
• The VPython extension gives us an easy way to produce graphs and
Python is an interpreted language, which makes it somewhat simpler to code
and run a program (there are no intermediate compile/link steps).
Interpreted languages like Python do tend to produce slow programs, how-
ever. it might be proﬁtable to write a program that take a long time to run
in some compiled language. Here are some comments on available languages.
Fortran The latest incarnation of Fortran is Fortran95 (standardised in
1995). Unfortunately Fortran compilers tend to be high-cost items. A sub-
set compiler, for a language called F, is available for free download from
www.meltingsand.com. This is suitable for writing new Fortran95 programs,
but not for compiling old FORTRAN77 code.
Fortran95 oﬀers many beneﬁts over languages like C which were not designed
for numeric computing. (Of course, these other languages do have some
features it would be nice to have in Fortran.) Some beneﬁts of Fortran are:
• Parameterized ﬂoating point precision.
• Flexible array handling.
• Arbitrary ranges in array subscripts.
• Array slices.
• Complex numbers.
PHY400W — cp 2003 6
• Intrinsic functions are part of the language, so greater optimisation is
• HPF extensions permit writing programs for parallel execution on sev-
• Arbitrary operators can be deﬁned.
Fortran is not object-oriented.
C It’s better to start with C++ these days.
C++ C++ is a powerful and complex language which permits object-
oriented programming. (This is carefully phrased. C++ also permits non-OO
The gcc compiler is freely available for many platforms. Some beneﬁts of
• Object oriented features: classes, inheritance.
• The STL (Standard Template Library) is part of standard C++.
Eﬃcient handling of matrices and arrays in numerical calculations seems to
require exotic template techniques which may not be handled well by many
Java Java is compiled to run on a “Java virtual machine” which may in
turn recompile the code into native machine code (‘just in time” compilation).
Some beneﬁts of Java are:
• Object oriented features: classes, inheritance.
• Eliminates some areas of C++ that can cause problems, e.g. pointers.
• Large set of supporting libraries, including platform-independent GUI
(graphical user interface) construction.
PHY400W — cp 2003 7
Some disadvantages of Java are:
• No overloading of operators. In particular, complex numbers become,
• Multidimensional arrays are not eﬃciently handled.
• There are problems associated with the Java representation of ﬂoating
Java is not recommended for numerical work.
B.2 Basics of programming
How do we transform a physics problem into something that can be solved
on a computer?
There are a number of steps in the process:
• Analyse the problem to ﬁnd a suitable numerical method of solution.
• Derive an algorithm to perform the solution.
• Write a computer program to embody the algorithm.
• Get the program debugged and running.
• Solve the problem!
An algorithm is essentially a set of steps that has to be performed in order
to accomplish a task — a recipe, if you like, that has to be followed. The
computer language that is used to translate this algorithm into a program
must oﬀer useful facilities that match the algorithm. These are usually in the
form of data elements, control structures and input and output facilities. In
addition, the language should enable the proper structuring of the program.
A useful approach to designing an algorithm is that of structured program-
ming and step-wise reﬁnement. An algorithm is broken down into a number
PHY400W — cp 2003 8
of broad steps, each of which can be represented by a procedure in a program-
ming language. Each of these can then be reﬁned (by a recursive process)
until the algorithm is done. (This is also known as a top-down approach).
At a higher level, one can take a procedural or an object-oriented approach.
The latter attempts to view the problem as a set of objects which can be
manipulated to solve the problem. This approach tends to be rather vo-
luminous, but does limit the potential interactions between parts of large
programs. Most simple problems, such as we will be dealing with, are prob-
ably best approached from a straight-forward procedural approach.
B.3 Data elements
With numerical programming, the basic data elements are of course numbers.
In physics, we deal with various kinds of numbers: integers, real numbers and
complex numbers. Computers, and computer languages, attempt to supply
some realisations of these numbers for us to use, but there are of course major
diﬀerences from the idealised mathematical quantities.
Integers are available in some limited range, with a maximum value of 231 − 1
for “32-bit” integers.
The “real numbers” available on a computer have little resemblance to the
real numbers of mathematics. The range is limited, they are not complete and
continuous, and the spacing between available numbers varied from number
to number. The system will be described below.
Complex numbers are found in Fortran, but not by default in most other lan-
guages. The recently standardised C++ has a complex library, but complex
numbers are not part of the language.
B.4 Control structures
Computer languages provide statements to implement these. There are a
number of control structures that are useful:
Sequence Usually statements are executed in sequence. Some languages
(e.g. HPF) can parallelize code so that the same code operates in
PHY400W — cp 2003 9
parallel on several processors.
Repetition These structures are generically known as loops. Fortran pro-
vides do loops, and C for loops.
Decision Logical choice is implemented by if/then/else structures. Multiple
choice is implemented as a case statement.
B.5 Higher-level structures
Programs are usually assembled out of higher level structures called proce-
dures. C provides a function for this purpose, while in C++ the basic unit of
structure is the class. A C program starts running in the main program. (A
C++ program is less obvious). Fortran has a speciﬁc program construct to
deﬁne where the program will start: it also provides two types of procedure,
and a module construct for providing cross-ﬁle structuring.
PHY400W — cp 2003 10
The purpose of computing is insight, not numbers.
— R.W. Hamming
Physics is a quantitative science: numbers obviously play a large role.
But many of the mathematical ideas that are used in physics are there to
get away from numbers and arithmetic. Much of what we do in calculations
involves the integral and diﬀerential calculus; in fact, solving many problems
in physics means ﬁnding the solution to a diﬀerential equation.
Of course, the “physics” in the problem is in setting up the equations, but
to see if the physics is correct we still have to solve them.
In solving problems on the computer, we deal with a simpliﬁed number sys-
tem, with none of the topological niceties of the real numbers like continuity
and completeness. So notions of diﬀerentiability do no enter: the number sys-
tem is discrete. We have to fall back on arithmetic in order to do anything,
and everything has to be done by tedious arithmetic manipulations.
(What of symbolic packages like Mathematica and Maple: don’t these allow
us to use the familiar tools of the calculus? Yes, but they are just a handy
way of remembering what the integral, say, of a function is. Think of them
as a handy table of integrals.)
C.1 The computer number system
There are two types of number represented on the computer: integers and
ﬂoating point numbers, the counterparts of real numbers. These are usually
represented in a binary manner in the hardware.
The representations depend also on the hardware. The PC will be used as
the example; other types of hardware might have diﬀerences in the ranges of
The integers we will use are usually stored in four bytes (assuming we are
dealing with “32-bit” programs!), each byte consisting of eight bits. 31 bits
are used to represent an integer, and one bit for the sign. Thus integers in the
PHY400W — cp 2003 11
range −231 to 231 − 1 (-2147483648 to 2147483647) are represented. (There
are usually other integers available, e.g 16-bit numbers, and some languages
might oﬀer signed and unsigned variants).
Real numbers are represented in the form of ﬂoating point numbers. These
are written in the form
sM B e−eo
where s is a sign bit, M is a mantissa, B is the base of the number system
(usually 2), e is an exponent and eo is an oﬀset. The mantissa consists of a
certain number of signiﬁcant digits (bits in a binary representation). There
is an implied decimal point at the start of the mantissa, and the exponent
is adjusted so that the leading digit (bit) is non-zero. (This adjustment
process is known as normalisation of the number; all ﬂoating point numbers
Thus, if we take a mantissa of 3 decimal digits and a signed exponent of 2
decimal digits as an example, a ﬂoating point number might be
−0.123 × 1001
Clearly, there is only a ﬁnite number of such numbers. These are only ap-
proximations of real numbers, and properties of the reals like continuity are
The ﬂoating point number will be stored in 32 or 64 bits (on the PC). One
bit is used for the sign, 23 or 52 for the mantissa and 8 or 11 for the exponent.
The leading bit is always set, and the exponent is adjusted while the bits of
the mantissa are shifted to achieve this. This process is known as normalisa-
tion. After an arithmetic operation, the ﬂoating point result is normalised.
If this cannot be done, an error condition (denormalised number) occurs.
Precision sign mantissa exponent
single 1 23 8
double 1 52 11
extended 1 64 15
The range of representable numbers on the PC is approximately
PHY400W — cp 2003 12
Precision minimum maximum sig. decimals
single 1.2E-38 3.4E+38 6-9
double 2.2E-308 1.8E+308 15-17
extended 3.4E-4932 1.2E+4932 18-21
Even if there is no error condition, there can be problems with the result.
Most important is the disappearance of low-signiﬁcance bits in the process
of subtracting two similar numbers.
We can get a feeling for this by considering base 10 with four digits in the
mantissa, i.e. take B=10 and, say, e0 = 50. Then the number a = 10.3213
can be represented as
+ 5 2 1 0 3 2
Note that there has already been some truncation of the number to ﬁt into
the ﬂoating point scheme.
The subtraction a = 10.3213 − 10.3122
+ 5 2 1 0 3 2 - + 5 2 1 0 3 1
results in the unnormalised number
+ 5 2 0 0 0 1
The normalised result is then:
+ 4 9 1 0 0 0
Note that signiﬁcant digits have been lost from the right hand side. This can
be a major problem in many calculations.