COMPUTATIONAL PHYSICS
UNIT 3/PH/SS (part 2)
Dr PA MULHERAN
AUTUMN TERM 2003
[Version date: 9 OCTOBER 2003]
5. INTRODUCTION TO THE UNIT
Introduction
This unit continues the theme of Computational Physics (3/PH/SS; Part 1) and students employ techniques from that Part to tackle more ‘computer experiments’. In addition new techniques, involving random processes and matrix methods, are introduced and used to explore topics in physics.
Part 2 comprises the remaining 5 of the 9 ‘computer experiments’ in 3/PH/SS. Each of these must be completed within two weeks. Each ‘computer experiment’ is described in a separate chapter of this manual and contains a series of exercises for the students to complete. Each student works alone and keeps a record of work in a logbook that must be submitted for assessment every two weeks. Note that the four projects 5-8 are compulsory; the remaining one is to be chosen from chapters 9, 10 or 11.
The Salford Fortran95 compiler will be used in this course, and this may be started by double-clicking on the Plato icon under the “Programming - Salford Fortran 95” program group. A “FTN95 Help” facility is supplied with this software and can be found within the same program group. This help facility includes details of the standard FORTRAN95 commands as well as the compiler-specific graphics features. All of the programs needed during this course may be downloaded from the Part III - Computational Physics page on the department’s web-server (
Web Site Information
In addition to all the chapters and programs required for this course, there are links to other useful sites including a description of programming style; a description of computational science in general and FORTRAN programming in particular; a tutorial for FORTRAN90; and a description of object-oriented programming.
References
Fortran 90 for Scientists and Engineers
By Brian Hahn
Published by Arnold £19.99
Fortran 95 Handbook
By Jeanne Adams, Walt Brainerd, Jeanne Martin, Brian Smith, and Jerry Wagener
Published by MIT Press, 1997. 710 pages.
$55.00
Fortran 90/95 Explained
By Michael Metcalf and John Reid
Oxford University Press ISBN 0-19-851888-9
$35.00
Fortran 90/95 for Scientists and Engineers
By Stephen J. Chapman
McGraw-Hill, 1998 ISBN 0-07-011938-4
$68.00
Numerical Recipes in Fortran90
By William Press, Saul Teukolsky, William Vetterling, and Brian Flannery
Published by Cambridge University Press, 1996. 550 pages. $49.00
A more complete list of reference texts is held at
where books can be ordered directly.
Logbooks
Each students is required to keep an accurate and complete record of the work in a logbook. In these logbook each student should answer all the questions asked in the text, include copies of the programs with explanations of how they work, and record details of the program inputs and of the output created by the programs. On completion of each chapter a students should write an abstract which summarises what has been achieved in the project..
We are often asked what should go in the logbook. It is difficult to give a simple answer to this since each computer experiment is different but as a guide you should have sufficient information for someone else
(i)to understand what you were doing and why; and
(ii)to be able to repeat the experiment.
Unit Assessment
The logbooks must be submitted for assessment to the School Office (Physics 217) by 1 pm on the Wednesday of weeks 2, 4, 6, 8 and 10. Late submissions will only be assessed in extenuating circumstances. Guidelines on the assessment are given below.
Feedback
In addition to comments written in the students’ laboratory notebooks by the assessors during marking, feedback on the projects will be provided by a class discussion and, when appropriate, by individual discussion between lecturer and student.
Assessment Guidelines
This unit is assessed completely by continuous assessment. Each project (which corresponds to one chapter of the manual) is assessed as follows.
The depth of understanding and level of achievement will be assessed taking into account the following three categories:
1. / Completion of the exercises (0 – 17 marks)- Completeness of the record
- Description and justification of all actions
- Following the documented instructions, answering questions and performing derivations etc
2. / Summary (0 – 3 marks)
Review of objectives
Summary of achievements
Retrospective comments on the effectiveness of the exercises
3. / Bonus for extra work (0 - 3 marks)
Any computational work beyond the requirements stated
An exceptional depth of analysis
An outstanding physical insight
The total mark of a project cannot exceed 20. Unfinished work will be marked pro rata, unless there are extenuated circumstances. The numerical mark will be translated to a grade:
The marks awarded for the best 8 of the 9 computer experiments will collectively form the overall assessment for this unit.
If you are unable to attend the laboratory session you should inform the School Records Administrator (Mrs S A Moon) by telephone 8543 (0118 931 8543) or by email
CONTENTS (PART 2 AUTUMN TERM)
- Analysis of Waveforms
- Random Processes
- Numerical Integration: Application to some Problems in Quantum Statistics
- Equilibrium and Temperature
- Eigenvalues and Eigenvectors of Matrices: Application to the Vibrational Modes of Linear Chains (option)
- Multilayer Dielectric Stacks (option)
- The Random Walk (option)
For the final ‘computer experiment’ select one from 9, 10 and 11.
1
5. ANALYSIS OF WAVEFORMS
Objectives
In this chapter the main elements of Fourier Analysis are reviewed and the methods are applied to some basic wave forms. On completion of this chapter students will have utilised simple numerical techniques for performing Fourier Analysis; studied the convergence of Fourier series and Gibb’s overshoot phenomenon; and investigated the best choice of Fourier coefficients in finite series. More advanced students will also apply the techniques to the mechanical vibration of a taut string.
Fourier Analysis
Fourier’s theorem states that any well-behaved (or physical) periodic wave form s(x) with period L may be expressed as the series
, / (5.1)where the Fourier coefficients are given by
/ (5.2)and
. / (5.3)Here the integrals are over one complete period (e.g. x=0 to x=L).
The nth component of the sum in (5.1) corresponds to a harmonic wave with spatial frequency and hence a wavelength .
In numerical work we can only deal with series with a finite number of terms. Suppose the first m terms of the series, are used to approximate the function, .
, / (5.4)Then the mean-square-error can be expressed as
/ (5.5)(5.6)
Thus the mean square error is a monotonically decreasing function of m (ie Em+1 Em) and sm(x) converges to s(x) as more terms are added to the series. It can also be shown that the mean-square-error Em (for any fixed m) is minimised by using the Fourier coefficients as calculated through equations (5.2) and (5.3).
The Numerical Methods
We now consider how to calculate the coefficients an and bn.
In general we can approximate an integral by means of the trapezoidal rule. The essence of this is shown in the diagram:
The integrand is divided into N equal intervals of size L/N and the integrand is approximated by a sequence of straight-line segments. The function s(x) is evaluated at the positions xk = kL/N.
The result of this procedure applied to (5.2) and (5.3) is
(5.7)
(5.8)
We can simplify this by defining an array
(5.9)
Using this array gives the result for the coefficients as
(5.10)
(5.11)
This procedure predicts coefficients for n < N/2. It fails to correctly predict coefficients for n>N/2. That is it fails to predict the Fourier components with spatial frequency greater than N/2L and wavelengths less than 2L/N.
In fact if the function s(x) has no spatial frequency greater than N/2L, the Sampling Theorem tells us that the expressions (5.10) and (5.11) are exact.
Now attempt the exercises.
Exercises
1. / [10 Marks](a) / Obtain a copy of the program Fourier.f95 and run it. The program sets up the various waveforms and plots them into the graphics window. Make sure you understand how it works, and record its main features in your log-book. Remember that the graph-plotting routine also writes to the clip-board, so that copies of the graphs can be pasted into Word and printed.
(b) / Write a subroutine to calculate the coefficients an and bn using equations (5.10) and (5.11), writing the results to a data file. Compare your numerical results with the analytical solutions to equations (5.2) and (5.3) for one of the wave forms.
(c) / What symmetry do the square and triangular wave forms have, and how is this related to the values of the Fourier coefficients? What about the ramp wave?
(d) / Use your numerical values of an and bn to reconstruct the approximation sm in equation (5.4), and plot this curve alongside the original wave form. The user should be prompted to select how many terms they want to include in the summation (i.e. choose the value of m). Note that the subroutine ‘graph’ provided will plot all the curves stored in y(npoints, ngraphs).
(e) / Investigate the convergence of the series; does that of the triangular wave converge faster than those of the square or ramp waves? Also investigate the Gibbs overshoot phenomenon observed in Fourier series for discontinuous wave forms.
(f) / Calculate the mean-square-error in equation (5.5); whole-array commands such as SUM will be useful here. Check that Em is indeed minimised by the Fourier coefficients calculated above in part (b), i.e. change one coefficient from its calculated value and show that the mean-square-error increases.
(g) / Show that Em monotonically decreases with m, and compare your numerical results with the analytical formula of equation (5.6).
Remember to keep an accurate record of your work in your log-book.
2. / [7 Marks]
Use your Fourier series for the triangular wave form to investigate how a plucked taut string vibrates. The region (0, L/2) corresponds to the ‘string’ and the region (L/2,L) corresponds to a ‘reflection of the string’. Sinusoidal standing waves with frequencies are supported on the string, where c is the phase (and group) velocity of the waves. To do this you must include the time dependence of each sinusoidal wave in the Fourier series, given that the string starts from rest with the triangular displacement.
What happens if damping is present?
Remember that you must finish your work on this chapter by writing a summary in your laboratory note books. This should summarize in about 300 words what you have learnt and whether the objectives of this chapter have been met.
6. RANDOM PROCESSES
OBJECTIVES
This chapter provides an introduction to random processes in physics. On completion, you will be familiar with the random number generator in FORTRAN 95, and will have gained experience in using it in two applications. You will also be ready to tackle later chapters that develop computational studies of random systems.
INTRODUCTION
It is convenient to describe models of physical processes as either deterministic or random. An obvious example of the former is planetary motion and its description via Newton’s equations of motion: given the position and momenta of the particles in the system at time t, we can predict the values of all the positions and momenta at a later time t’. Even the solution of Schrödinger’s equation is in a sense deterministic: we can predict the time evolution of the wave function in a deterministic way even though the wave function itself carries a much more restricted amount of information than a classical picture provides.
The obvious example in physics of a theory based on randomness at the microscopic level is statistical mechanics. There may well be deterministic processes taking place, but they do not concern us because we can only observe the net effect of a vast number of such processes, and this is much more amenable to description on a statistical basis. But a statistical basis does not only concern statistical mechanical (and thermodynamic) systems. Many physical systems are inherently disordered and defy a simple deterministic analysis: the passage of a liquid through a porous membrane (oil through shale, for example), electrical breakdown in dielectrics, the intertwining of polymer chains, and galaxy formation are some examples of random processes.
Statistical mechanics uses concepts like entropy, partition functions, Boltzmann, Fermi or Bose statistics and so on to describe the net effect of random processes. In computer simulations, one actually models the microscopic random processes themselves. To model randomness, we need to have something to provide the element of chance - like a coin to toss, or a dice to throw. Of course, in computing, we do not use coins or dice but rather random number generators to inject the statistics of chance, and we start by seeing how they work.
RANDOM NUMBER GENERATORS
The basic algorithm
Random number generators are more precisely known as pseudo-random number generators. The sequence of numbers they produce can be predicted if the algorithm is known, but there should be no correlations between the numbers along the sequence. In practice, the sequence will repeat itself but the period of the cycle should be longer than the equivalent scale of the process one wants to simulate. Random number generators are based on the algorithm
where xn+1 and xn are respectively the (n+1)th and nth numbers in the sequence. All the quantities in the expression, including the numbers themselves and the constants a, c, m are integers. y = z mod m means that y is the remainder left after dividing z by m. For example 27 mod 5 equals 2.
Now go to Exercise 1(a)
You will see, from the Exercise, that with a judicious choice of parameters we can produce a pseudo-random sequence of integers i such that 0 i < m (note 0 appears in the sequence but m does not). Usually real random numbers r between 0 and 1 are generated. This can be done using the algorithm with a final step r = REAL(i)/REAL(m), such that 0 r < 1.
Intrinsic Subroutine
All computing systems have a built-in random number generator (usually based on more than one of the basic generators just studied) that has been optimised (hopefully), and one would normally use that. For FORTRAN 95, the simplest use is as follows:
CALL RANDOM_NUMBER(r)
r is the generated random number (with 0r<1), and r must be declared as REAL. You can also declare r as a one dimensional real array; in this case the subroutine returns a random number in each element of the array.
If you repeat a run of a program containing this call, the same set of random numbers is produced. It is a good principle to seed the random number generator at the start of the program using the system clock as in the following example, which produces a different set of random numbers each time it is run (starts at a different place in the sequence).
INTEGER :: count
REAL :: r
INTEGER, DIMENSION(1) :: seed
CALL SYSTEM_CLOCK(count)
seed = count
CALL RANDOM_SEED(PUT = seed)
WRITE(*,*)’Random number seed = ‘,count
CALL RANDOM_NUMBER(r)
Note, count (the current value of the system clock) is integer and seed (used to seed the random number generator) is a rank-one integer array. It is sensible to write the seed value in case you wish to rerun the program with the same set of random numbers. If you do, you can use the same piece of code with the system_clock statement replaced by READ(*,*)count.
Different number ranges
Often we want to generate random numbers over a range different from 0 to 1. This is straightforward. If we want a real random number x between –2 and +2, for example, this can be obtained from the random number generator output r using the statement: x = -2.0 + 4.0*r. Generally, use the expression x = a + (b – a)*r if the range is a to b.
We have to be a little more careful with integers. Suppose we were simulating the throw of a dice (and needed to generate a random integer from the set 1, 2, ….,6). If d is the required random integer, we can use the statement: d = INT(6*r) + 1. Calculation of INT(6*r) gives one of the integers 0, 1, ….,5.
Now go to Exercises 1(b) and (c)
MONTE CARLO INTEGRATION
There is a wide range of calculations in Computational Physics that rely on the use of a random number generator. Generally they are known as Monte Carlo techniques. We will start by looking at the Monte Carlo method applied to integration. There are two approaches to choose from: the ‘hit and miss’ method and the ‘sampling’ method.
Hit and Miss Method
Here is analogy to help explain the method. An experimental way to measure the area of the treble twenty on a dart board is to throw the darts at the board at random; if N60 is the number hitting the treble twenty and N is the total number of darts thrown, then the area A60 of the treble twenty is given by the equation A60 =A*N60/N, where A is the total area of the board.
Now go to Exercises 2(a) and (b)
Sampling Method
This method can be summarised by the equation
Choose N random numbers xi in the range a<xi<b, calculate f(xi) for each, take the average, and then multiply by the integration range (b-a). It is similar to Simpson’s rule but in that case the values of xi are evenly distributed. We have seen already how to generate random numbers in the range a to b.
The extension to a multidimensional integral is easy. For example, in two dimensions, write
In this case choose N pairs of random numbers (xi,yi) and go through a similar procedure. The extension to an arbitrary number of dimensions is straightforward.
Generally, from the point of view of accuracy, it is better to use ‘conventional’ methods like Simpson’s rule for integrals in low dimensions and Monte Carlo methods for high dimensions. If you have body with an awkward shape, however, Monte Carlo methods are useful even at low dimensionality.