University of Cambridge Department of Physics

Computational Physics Question Sheet

Michaelmas 2007

Computational Physics Question Sheet You should choose one of the following five questions to answer. The aim is to develop and test a program which can be used in the solution of the problem. You should use library subroutines to perform the mathematical components wherever possible. Keep a record of the analysis you do before programming, details of the library routines you use, and the tests that you make, as well as the final commented program and results. The written up form of these, (possibly together with an oral session with a head of class), will form the basis of the assessment.

The Write-up The following is intended as a guide to the structure of your report, and what should be included in it. You should follow the basic structure suggested in the booklet “Keeping Laboratory Notes and Writing Formal Reports” and the general advice given there on style, except as noted below. 1. Content The aim is to produce a report on your solution to one of the set problems. The problems are intended to be exercises on the material taught in the course and are not “projects”. What you should concentrate on in your write up, therefore, are those aspects of computational physics relevant to the problem you have chosen. The physics of some of the problems may go beyond your current knowledge, and understanding the advanced ideas behind some of the problems is not part of this exercise. However, basic physical reasoning and the physics that you have already met in Part IB and this term’s Part II courses will be assumed. You should include the following: 1. A brief introduction to the problem (taken from the question sheet mainly). 2. The analysis of the computational aspects of the problem that were necessary for you to undertake its solution. This might include a brief description of relevant computational physics, choice of approach, choice of suitable routines, scaling of the problem, etc. 3. The implementation: that is, your program. The full program listing should be included in an appendix. In the main text your description of the program should be relatively short, concentrating on your overall approach and the way you implemented the necessary algorithms to solve the problem. 4. Results and analysis. This should include discussion of the tests you performed to demonstrate that your program is working correctly. Your final results should also be presented together with a discussion of errors and any other computational problems. 2. Structure A possible structure might be the following (although details are problem specific): Abstract; Introduction; Analysis; Implementation; Results and Discussion; Conclusions. -1-

3. Marking This is, of course, an assessed piece of work. Approximately equal credit will be given to each of the following areas: 1. Analysis of the computational physics aspects of the problem (see 1.2 above). 2. Implementation of the algorithm (see 1.3 above). 3. Results, analysis of errors (if applicable), tests and discussion of the relevant computational physics (see 1.4 above). 4. Overall presentation of the report, structure, etc. A number of candidates may be selected for Viva voce (i.e. oral) examination early in the Lent Term after submission, as a matter of routine, and therefore a summons to a Viva should not be taken to indicate that there is anything amiss. You will be asked some straightforward questions on your project work, and may be asked to elaborate on the extent of discussions you have had with other students. You may also be asked to give a “live demonstration” of your program. So long as you can demonstrate that your write-up is indeed your own, your answers will not alter your projects marks.

4. Collaboration The remarks about collaboration and cheating which are contained in the course handbook, with reference to experimental work, also apply here. Discussion of the problems among students and with demonstrators is encouraged, as this will often help your understanding. However, your program and your report should be written by you, and only results produced by your program should be presented in your report. It would be unacceptable for two students to submit near identical programs. Any attempt to pass off the work of others as being produced by yourself will be regarded as cheating and may be referred to the examiners. When the report is handed in, you will be asked to sign a declaration about the degree of collaboration that went in to the work that is being submitted.

-2-

Question 1. Fourier transforms and signal processing Big Ben, the bell in the Parliament Tower, has a large number of resonant modes which are excited when it is rung. The response can be modelled very simply by treating it as the linear superposition of the responses of a number of damped harmonic oscillators. The data in the file: http://www.tcm.phy.cam.ac.uk/teaching/compphys/bigben.dat also available on the PWF Linux systems as: $PHYTEACH/part_2/bigben.dat contains a noisy 16-bit digital sample (at 8 kHz) of one ring. Write a program which reads in this data and then Fourier transforms it to obtain the power spectrum of the sample and the auto-correlation function, both of which you should plot. Investigate the effects of different filters applied to the data. Extend your program so that you can apply a variety of filters to the data in the Fourier domain and recalculate the data in the time domain. Use your program to isolate the response which corresponds to a peak in the power spectrum occurring close to 200 Hz and investigate the effect of different filters on the filtered time-domain data; you should consider appropriate high- and low-pass filters, top-hot, triangular and gaussian filters. In your report you should consider how the form of the filter effects the timedomain data. You should then be able to compare estimates of the Q-factor of Big Ben obtained in two different ways: (i) from the power spectrum in the frequency domain, (ii) by fitting an exponential function to the amplitude decay of the filtered ring in the time domain.

-3-

Question2. The Trojan asteroids Jupiter

Greeks

Trojans

π/3 Centre of mass

SUN

Accompanying the planet Jupiter in its orbit around the Sun are two groups of asteroids, named individually after the main protagonists in the Trojan wars. The two groups are constrained at points of stable equilibrium (the Lagrange points) preceding and trailing the planet at an angular distance of π/3 radians. The combination of the gravitational attraction of the Sun and Jupiter gives a resultant force on a Trojan asteroid towards the centre of mass of the two bodies such that the centripetal acceleration causes the asteroid to orbit with the same period as the planet. Write a program to solve the equations of motion for this system numerically, using an adaptive Runge-Kutta technique, to demonstrate that the asteroids will stay fixed at the Lagrange points over many orbits of Jupiter around the Sun. This should confirm the analytical solution of the equations, and show that your program is operating correctly with no numerical instabilities being propagated during the computation. The next stage is to vary the initial positions of the asteroids slightly to investigate the stability of the equilibrium positions. Plot the positions of the asteroids through a few hundred orbits, and show that they oscillate about the Lagrange points, but do not escape from the stable position. You can also derive a value for the total range of wander of the asteroid, taking the average of the leading and trailing groups. Trojan asteroids are also present in the orbit of Mars. Run your program again for a range of other planetary masses, within one or two orders of magnitude of that of Jupiter, and derive a relationship linking the range of wander and the planetary mass. Suitably scaled units for this problem are solar system units, where the mass of the Sun (M~) is taken as unit mass, unit distance is the astronomical unit (AU) and unit time is one year. For such a system the gravitational constant G is 4π2, the mass of Jupiter is 0.001 M and the average distance of Jupiter from the Sun is 5.2AU. The ~ masses of the asteroids can be regarded as negligible in this system.

-4-

Question 3. Localisation in a linear chain A set of masses mj, (j=1,N), is coupled by a collinear set of springs with spring constants kj, (j=0,N), where the extremal ends of the springs 0 and N are fixed, and the whole chain is in equilibrium as in the diagram below.

k0

m1

k1

mN -1

kN -1

mN

kN

Write a program to calculate the eigenfrequencies (wr) and normal modes (normalised eigenvectors) (q(wr)) of this system. A possible coordinate system might consist of the set of longitudinal displacements of each mass from its equilibrium position. An appropriate normalisation for the eigenvectors is: N

∑ m q (w ) q (w ) = 1 j

j =1

j

r

j

r

Plot the spectrum of eigenfrequencies and selected normal modes for some simple configurations of masses and spring constants. One possible way to display a normal mode is to plot the square of each component against its index. For the normalised eigenvectors ‘localised’ modes may be thought of as those for which only a contiguous subset of the components of q(wr) are significantly different from zero. It is expected that for such states: q j (w ) ~ exp(− 2γ l j − jd 2

)

for some jd and γl . γl can be used as a measure of the ‘localisation length’ of the mode. By taking a chain of constant masses and fixed spring constants and varying the magnitude of a central mass, extend your program to investigate the frequencies and localisation lengths of the subsequent normal modes. Perform a similar study for a ‘diatomic’ chain of alternating masses, with fixed spring constants and an isolated variable mass. The program can be developed with quite small values of N (~10), and localisation effects should be detectable for N less than 100 for suitable values of the masses and spring constants.

-5-

Question 4. Percolation on a 2D lattice The aim of the problem is to investigate percolation on a 2D cartesian lattice. The first stage of the problem is to develop a program which implements a simple percolation simulation. In practice you will need to run the simulation many times and therefore a reasonably simple and efficient algorithm needs to be developed. The second stage of the problem is to determine the percolation threshold to at least two decimal places; you need to consider carefully the statistical nature of the simulation and the effect of lattice size. Stage 1: percolation algorithm and the percolation threshold We will investigate the percolation threshold on a 2D cartesian lattice. To do this we will use a Monte Carlo approach and model a “forest-fire” type of problem as mentioned in the lectures. The specific algorithm is as follows: 1. Create a lattice of sites containing empty sites (0) and a proportion, p, of filled sites (1) which represent trees; the distribution of sites must be random. 2. Two other site types are defined: 2 = burning tree; 3 = burnt tree 3. Initialise the simulation by setting alight all trees in the first row of the simulation (site of type 1 become type 2). 4. Proceed iteratively passing through the lattice row-by-row and apply the following rule: If a lattice site is a tree (type 1) then if it has any first-nearest neighbours which were burning (type 2) during the last pass, ignite the tree; if a lattice site contains a burning tree (type 2) set it to be burnt (type 3). 5. Repeat this procedure until either the fire reaches the last row or the fire is extinguished; record the lifetime of the fire, τ, defined as the number of complete passes through the lattice for this to happen. Having developed the algorithm and tested it on a reasonably small grid (say 50 by 50 cells) you should perform some preliminary investigations of how the lifetime of the fire depends on p, and get some idea of the statistical uncertainty in this number.

Stage 2: determining the percolation threshold pc Using your preliminary investigation you should design a numerical experiment to determine pc, which in this case occurs when the lifetime of the fire, τ, becomes a maximum. You should consider carefully the errors and also the effects of the size of the lattice on your estimate of pc.

-6-

Question 5. Tidal Tails of Interacting Galaxies The aim of the problem is to develop a simple N-body simulation of interacting galaxies. The problem will be simplified by considering two masses and a number of test masses. You should setup the N-body problem to use the Verlet integration scheme discussed in lectures. You should approach the problem in two stages. Stage 1: creation of a single model galaxy The model galaxy consists of a central mass and five rings of test particles. Consider carefully the scaling of the problem. On possibility is to use a set of units such that the central mass has a mass of 1 unit and G=1. Arrange the test particles in rings about the central mass at radii 2, 3, 4, 5 and 6 “units”. An appropriate density of particles is 12, 18, 24, 30 and 36 on the five rings respectively. Since they are to be treated as test particles the only force felt by each particle is from the central mass alone. Determine appropriate velocities for the particles and setup an N-body simulation so that the particles correctly remain on circular orbits. Stage 2: interaction Stage 2 involves the introduction of a second non-test mass also of mass unity. This mass represents a perturbing galaxy orbiting around the first. Determine the initial conditions of this second galaxy so that it will pass, approximately, a distance of 9 to 12 units from the first central mass at closest approach; the orbit should be chosen to be parabolic in the first instance (i.e. Etotal = 0). Test the orbit of the two masses without the test masses. Re-introduce the test masses which now see a force due to the two non-test masses. Follow the evolution of the system by outputting the positions of all of the masses as a function of simulation time (you will not need to do this for each iteration). Plot the results.

-7-