TCD School of Mathematics Module MA3469

Practical Numerical Simulations

Click here for a sample exam paper

Part I. Ordinary Differential Equations.

1. 1D motion: dimensionless radial equation for Kepler problem. Types of motion and trajectories in the phase space. Numerical analysis of the problem:

1a. Euler method and its accuracy.
1b. 2nd-order Runge-Kutta method and its accuracy.
1c. Euler-Cromer/leapfrog method and its accuracy.

Program kepler1d.cpp.

Demonstration how to compile it, run and make plots.

Assignment 1, due 12th October 2010.

2. 2D motion: full Kepler problem. Analysis of the scaling behaviour of the problem and introduction of appropriate scales for time and distance. Transition to dimensionless variables. Numerical analysis of the problem:
2a. Euler-Cromer method and fundamental benefits of symmetric methods.
2b. Reduction of the problem to a vectorised system of 1st-order ODEs. Importance of the translation table between the physical variables and the vector of unknowns.
2c. Vectorised Euler, 2nd and 4th-order Runge-Kutta codes for the problem.

Programs kepler2d.cpp, kepler2fun.cpp, kepler_RK2.cpp.
A couple of hints on compilation and making plots.

2d. Importance of error control in numerical ODE solvers. Two approaches to estimating the errors: embedded Runge-Kutta methods and Richardson extrapolation.
2e. Error Lemma. Obtaining the step error and improved solution using a simple 2-point Richardson extrapolation.
2f. Practical strategies of adaptive step control.

Program kepvar_RK2.cpp and the relevant plots.

Assignment 2, due 26th October 2010.

2g. Richardson extrapolation over N+1 points.
2h. Error expansion for symmetric methods. Modified midpoint method and polynomial Bulirsch-Stoer approach.
2i. Implementation of the polynomial Bulirsch-Stoer method. Neville algorithm. Practical strategies for adaptive step and order control.

Fixed-order Bulirsch-Stoer program kepler_BS_fixedN.cpp.

Adaptive-step fixed-order Bulirsch-Stoer program kepler_BS_fixN_varH.cpp.

Assignment 3, due 16th November 2010.

Part II. Partial Differential Equations.

3. One-dimensional Schroedinger equation as an example of a parabolic PDE. Analysis of the scaling behaviour of the problem and introduction of appropriate scales for time and distance. Schroedinger equation in dimensionless variables.

3a. Spatial discretisation of the problem. Centred expression for the 2nd-order space derivate and its error estimate. Discretised Hamiltonian in matrix notation and the role of boundary conditions.
3b. Discretisation of the time evolution of the problem. Right (Forward), Centred and Left (Backward) expressions for the 1st-order time derivative and their error estimates. Forward-Time Centred-Space (FTCS) differential scheme for Schroedinger equation, its space-time diagram and matrix notation.
3c. Numerical study of a one-dimensional quantum harmonic oscillator. Relevant physical scales and the appropriate choice of the coupling constant, lattice spacing and time step.

A little horror movie "Explosion In A Potential Well" has been created using the program ftcs.cpp and two gnuplot scripts, movie.gpl and loop.gpl:
           part I (ht/h2=1, 240KB),
           part II (ht/h2=0.1, 760KB),
           part III (ht/h2=0.01, 1.6MB),
           part IV (ht/h2=0.001, 1.7MB).

3d. Von Neumann stability analysis for FTCS in the imaginary (Schroedinger equation) and real (diffusion equation) time.
3e. Backward-Time Centred-Space (BTCS) differential scheme as a simplest example of an implicit scheme for Schroedinger equation. Space-time diagram and matrix notation of BTCS. Proof of its unconditional stability both in real and imaginary time.
3f. Crank-Nicolson differential scheme and its accuracy in comparison to FTCS and BTCS. Space-time diagram and matrix notation for Crank-Nicolson scheme. Stability analysis in imaginary and real time.
3g. Practical implementation of implicit differential schemes using tridiagonal Gaussian elimination and backsubstitution.

Example: implicit schemes such as BTCS (btcs.cpp) and Crank-Nicolson (cn_map.cpp) are, unlike FTCS, unconditionally stable. The latter program generates datafiles suitable for making 2D probability maps with gnuplot, see the script map.gpl.

Assignment 4, due 30th November 2010.

3h. Unitarity of the continuous and discretised Schroedinger equations. Violation of unitarity by FTCS and BTCS schemes. Proof of unitarity of Crank-Nicolson scheme.

Part III. Monte-Carlo methods.

4. Langevin and reduced Langevin equations for Brownian motion. Time-discretised reduced Langevin equation and its interpretation in terms of random walks.

4a. Statistical properties of the random force.
4b. Random number generators and statistical properties of uniform random numbers. Obtaining random numbers of any given distribution. Box-Muller algorithm for generation of normally distributed random numbers.

Demonstration how to use uniformly distributed random numbers. The programs fixed_state.cpp, random_state.cpp, uni_hist.cpp and the relevant plots.

4c. Random walks in the double-well potential. Analysis of the physical scales of the problem and introduction of appropriate dimensionless variables. Choice of the time step and effective temperature.
4d. Physical observables in double-well random walks. Transition time over the potential barrier and its effect on particle distribution between the wells. Rapid thermalisation of particle distribution within a single well.

Two movies about how the particles distribution in a two-well potential thermalises starting from a symmetric and asymmetric initial state. The movie has been generated with walks.cpp program and the two usual (but slightly modified) gnuplot scripts, movie.gpl and loop.gpl.

4e. Emergence of Boltzmann distribution in many-particle random walks. Equivalence of single-particle and many-particle distributions. Role of normalisation factor for Boltzmann distribution. Finite-time-step artifacts in numerically observed distributions.

Similar movies with two distributions at different temperatures, in the linear and logarithmic scales.

4f. Generation of a thermalised particle ensemble using the Metropolis algorithm. Its applications to numerical studies of statistical systems. Simulated annealing method of solving minimisation problems.