Physics 329: Introduction to Computational Physics: Course Notes
This document will be updated throughout the course
Index
- Unix (Aug 28, Aug 30, Sep 4)
- Maple
- Lecture 1 (Sep 6)
- Worksheet (PS)
showing calculations in Chapter 2 of Maple V Learning Guide
by Heal et al
- Lecture 2 (Sep 9)
- (1) Useful Maple commands, (2) Maple V.4 features discussed in
Maple V Programming Guide which are not available in our
version, (3) Common mistakes, things to watch for.
(PS)
- Simple Maple source file (text)
- Sample Maple procedure: ladd which adds elements of numeric
list (text)
- Sample Maple procedure: series_op which computes expansion
coefficients for various power-series operations
(text)
- Building series_op interactively.
(PS)
- Lecture 3 (Sep 11)
- More on the series_op example
(PS)
- Maple procedure polyinterp for computing Lagrange
interpolating polynomial
(PS)
- polyinterp examples
(PS)
- Lecture 4 (Sep 13)
- polyinterp ``sin-fit'' example: Maple session
(PS)
and plot showing sin(x) and polynomial fit
(PS). Select 'Landscape'
from 'Orientation' menu if viewing plot with ghostview.
- Additional Material
- Annotated maple
source file to accompany Chapter 1 of Maple V Programming Guids
- Fortran Programming
- Most of the source code covered in class can be found on
the phy329 account on einstein in
~phy329/f77/ex1, ~phy329/f77/ex2, ... (organized by lecture).
- IMPORTANT: Refer to on-line (postscript) version of SGI
Fortran 77 Language Reference Manual
for full details of Fortran syntax.
- Lecture 1 (Sep 16)
- Lecture 2 (Sep 18)
-
fdemo1.f: Program demonstrating some essential
Fortran 77 language elements.
-
fdemo1_output: Output from fdemo1
- Lecture 3 (Sep 20)
- Lecture 4 (Sep 23)
- Lecture 5 (Sep 25)
- Sample output
from nth, a filter (shell-script) which
selects specified columns (separated by white space)
from standard output and writes them to standard
output.
-
arraydemo.f Demonstrates a general technique
for writing Fortran programs which define and manipulate
multi-dimensional arrays whose bounds are determined at
run time.
Program output
- Lecture 6 (Sep 27)
- Lecture 7 (Sep 30)
- dmroutines.f:
Two dimensional array (matrix) routines:
- dmmmult: Matrix-matrix multiply (square matrices only)
- dmfrom: Read matrix from file or stdin
- dmto: Write matrix to file or stdout
- tdm.f: Test
driver for dmroutines.f.
Program output on SGIs.
Program output on Cray J90.
Note
that this last document also contains general instructions re
how to 'port' a Fortran code to the Cray.
- Linear Systems
- Most of the source code covered in class can be found on
the phy329 account on einstein in
~phy329/linsys/ex1, ... (organized by lecture).
- General References
- Lecture 1 (Oct 2)
- dgesv.f:
LAPACK driver routine for solving general linear system via
LU decomposition with partial pivoting.
- tdgesv1.f: Test
program illustrating use of dgesv.
Makefile
and output on SGIs.
Output on Cray J90.
- Maple worksheet (PS)
which verifies above test.
- Lecture 2 and 3 (Oct 7 and 9)
- dgtsv.f:
LAPACK driver routine for solving general tridiagonal system with
no pivoting.
- bvp1d.f: Test
program illustrating use of dgtsv to solve a 1d
boundary-value problem using second-order finite-difference techniques.
Makefile
and sample output on SGIs.
Plots showing level-8
solution and
error for levels 5, 6 and
7.
- Lecture 4 (Oct 11)
- Lecture Notes:
Instructor version (PS) and
student version (PS).
- dgbsv.f:
LAPACK driver routine for solving general banded system via
LU decomposition with partial pivoting.
- bvp1d4.f: Test
program illustrating use of dgbsv to solve a 1d
boundary-value problem using mixed second and fourth order
finite-difference techniques.
Makefile
and sample output on SGIs.
Plots showing level-4
solution and
error for levels 4, 5 and
6.
- Finite Difference Methods
- Most of the source code covered in class can be found on
the phy329 account on einstein in
~phy329/fd/ex2 and
~phy329/fd/ex5.
- General References
- Lecture 1 (Oct 14)
- Lectures 2 and 3 (Oct 16)
- rexbvp1d.f:
Program illustrating use of Richardson extrapolation to
produce O(h^4) results from base O(h^2) approximation to
simple boundary-value problem.
Makefile
and sample output on SGIs.
Plot showing level-4
solution.
- Lecture 5 (Oct 23)
- Correction to lecture notes concerning determining
initial data for finite-difference version of wave
equation (PS).
- gpwave.f:
Program illustrating output of time series of one-dimensional
profiles (i.e. f(x,t) for subsequent plotting (surface
plot) using gnuplot.
Sample output, gnuplot
command file, and
plot (PS). Use help
splot in gnuplot for more information on making
surface plots.
- Root-Finding and Non Linear Equations
- Source code covered in class can be found on
the phy329 account on einstein in
~phy329/nonlin/ex1, ... (organized by lecture).
- General References
- Lecture 1 (Oct 25)
- Lecture 2 (Oct 28)
- newtsqrt.f:
Fortran program for computing square roots (solving x^2 - a = 0)
using Newton's method.
Sample output on SGIs.
- Lecture 3 (Oct 30)
- newt2.f:
Fortran program for computing root of simple set of non-linear
equations (d=2) discussed in class.
Sample output on SGIs
and a check on the
computation using Maple's numerical root finding capabilities
(fsolve).
- Cellular Automata
- General References
- Statistical Mechanics of Cellular Automata,
Stephen Wolfram, Rev. Mod. Phys, 55, 601-644 (1983)
- Cellular Automata and Complexity, Stephen Wolfram,
(collected papers), Addison-Wesley (1994)
- Source code covered in class can be found on
the phy329 account on einstein in
~phy329/ca/ex1
- Lectures 1 and 2 (Nov 1)
- ca1dr1.f:
Fortran program for simulation of two-state, 1d cellular automata
(CA) with range-1 update rule. Sample
output on SGIs.
- ca1dr2t.f:
Fortran program for simulation of two-state, 1d cellular automata
(CA) with range-2 totalistic update rule. Sample
output on SGIs.
Higher resolution
image
of simulation using Rule 12: graphics via SGI gl library.
- Supplementary Material
- Particle Simulations
- General Reference
- Computer Simulation Using Particles,
R. W. Hockney and J. W. Eastwood, IOP Publishing, (1988)
- Source code covered in class can be found on
the phy329 account on einstein in
~phy329/part/ex1
- Lecture 1 (Nov 4)
- Sample output illustrating
use of some utility commands: dvmesh, nf and
paste, which be useful for doing certain homework problems.
- Lecture 3 (Nov 13)
- nbody.f: Solves
gravitational n-body problem using second-order finite difference
techniques, and constant global time step.
Sample output on SGIs.
PS Plot of particle traces
for evolution of equal-mass binary in circular orbit
(about 1/3 of an orbit).
Initial data for binary
orbit.
- Supplementary Material (Nov 20)
- Solution of ODEs
- General Reference
- Source code covered in class can be found on
the phy329 account on einstein in
~phy329/ode/ex1 and ~phy329/ode/ex2
- Lectures 1 and 2 (Nov 18 and 20)
- lsoda.f ODEPACK
routine for solving general systems of ODEs.
- tlsoda.f Sample
driver routine which uses lsoda to integrate
u''(x) = -u(x). Sample
usage on SGIs.
Plot of computed
and exact solution on x = [0 .. 15] with all error tolerances
set to 1.0d-6. Plot of
error in computed solution.
- Lecture 3 (Nov 22)
- integral.f:
Driver routine illustrating use of LSODA (ODE integrator)
to compute definite integral.
fcn.f:
User function which defines integrand.
Makefile
- 2body.f:
Driver routine illustrating use of LSODA to solve restricted
2-body gravitational problem (one body non-gravitating).
fcn.f:
User function which implements equations of motion.
Makefile
- Lecture 4 (Nov 25)
- Equations of motion
for orbiting dumbbell problem.
- dumb.f:
Driver routine for solution of orbiting-dumbbell problem.
fcn.f:
User function which implements equations of motion.
Makefile
- Plots of omega (angular frequency of dumbbell about its
center of mass) versus time for a
circular orbit
and for an
elliptical orbit.
Note the "chaotic" nature of omega in the latter case.
Detail
of omega for the elliptical orbit.
- Plots illustrating energy conversation for elliptical
orbit (chaotic case).
Energy quantities
computed with LSODA tolerance 1.0e-6.
Top to bottom: Translational
kinetic energy, rotational kinetic energy, total energy,
gravitational potential energy. In this case, "non-conversation"
of total mechanical energy is a good sign that we need to
make the error tolerance more stringent.
Energy quantities
computed with LSODA tolerance = 1.0e-12. Note that
energy conservation is improved in this case.
Plot of
rotational kinetic energy.
- Monte Carlo Techniques: Random Numbers and Stochastic Simulations
- General Reference
- Source code covered in class can be found on
the phy329 account on einstein in
~phy329/rand/ex1 and ~phy329/rand/ex2
- Lecture 1 (Dec 2)
- tsrand.f:
Demonstration program showing use of srand() to set
seed for uniform random number generator rand().
Plot
of first 50 random numbers generated using tsrand
with seeds 0, 1, 2 and 3.
Note that there may be excess correlation between sequences generated
with seeds separated by some constant amount, so seeds themselves
should generally be chosen randomly. As pointed out in Numerical
Recipes, it is healthy to always view system-supplied random number
generators with a certain amount of suspicion, particularly
if you are going to generate a lot of random numbers.
However, I feel that the SGI supplied routine will suffice for our
purposes.
- nurand.f: Fortran
routine for generating non-uniform deviates given user-supplied
probablility-distribution function (PDF). Sample PDFs,
pdfs.f
for generating uniform and gaussian-distributed deviates.
tnurand.f:
Driver program for nurand.f.
Plot
showing normalized counts (1000 bins) for gaussian-distributed
deviates---blue: n = 100 000 random numbers, red: n = 10 000 000
random numbers, dashed black line: peak of analytic PDF.
Doit:
C-shell script which performs runs then
invokes sm with command file
sm_nurand
to create plot.
- Lecture 2 (Dec 4)
- dla.f: Program
which performs Diffusion Limited Aggregation (DLA) simulation
in 2D with variable non-stochastic central force (as discussed
in class). This central "bias" speeds cluster formation but
alters the grown cluster's "fractal properties" (for example,
the cluster's radius grows more slowly with particle number
for any non-zero bias than for zero bias). Sample
usage on SGIs.
- Plot
of final state in pure DLA simulation with 10 000 particles.
Detail of above simulation.
Same
simulation
after an additional 50 000 particles have
been added. Here, the effect of launching all particles
from a fixed, finite radius is apparent.
Bias = 0.0625
simulation after 10 000 steps. Note how structure is
already slightly less "distributed" than the pure DLA
case.
Bias = 0.1
and
Bias = 0.9
simulations, also after 10 000 steps. It is left as an
exercise to the reader to investigate and/or predict run
time as a function of bias.
- A nice DLA
specimen
from the Laboratory of Computer Science at MIT.
- An "embedded"
specimen
from
Eric Weeks'
(UT grad student) collection of
algorithmically-generated pictures.