PHYS 580 COMPUTATIONAL PHYSICS SPRING 2005

PHYS 580 COMPUTATIONAL PHYSICS

Sample programs for numerical algorithms

Clink on the links below; you can save these files to disk and compile them.

Note: very important! To view properly, you must view as page source (open the file, and then ctrl+u) You can also right-click on the link and choose "Save Target As..."

Differentiation

Simple program for differentiation

Same program but somewhat better written

Compare 2-point and 3-point differentiation as a function of x Graph for dx = 0.3 for different x

Compare 2-point and 3-point differentiation as a function of dx Graph (of “error”) for x = 1.0 for different dx

Another code comparing 2-point and 3-point derivatives

Exercise: use the above programs (or write your own) and using xmgrace or some other graphing program and try to make graphs similar to those above. Try to write as either a postscript file or as a jpeg file (hint: you have to use “print” you cannot simply “save”).

Extrapolation on dx assuming a power law (dx)^k

Quadrature (integration of functions)

Trapezoidal rule

Comparison of trapezoidal and Simpson’s rules: quadmaster.f and quadraturesubs.f

(to compile, do f77 –o qtest.x quadmaster.f quadraturesubs.f or use your favorite compiler)

Comparison of trapezoidal and Simpson’s for integrating x*exp(-x) on a linear scale on a log scale

Exercise: Modify the above program but compare for a different function (you choose). Make graphs.

Exercise: Write a subroutine for Bode’s rule and add to above program. Make graphs.

Gauss-Legendre quadrature: calling program and Gauss-Legendre subroutine (compile together)

Random number library (ranlib.f) and example calling routine

Plot of randomly distributed numbers: uniform (ran3) and Gaussian (gaussvar)

Exercise: modify the example calling routine above but instead of computing the average and std deviation of a uniform distribution (ran3) use the Gaussian distribution gaussvar . Do you get the average and (importantly) the std dev you expect? For experts: Compute the 4^th moment, both analytically and numerically.

Monte Carlo quadrature: (all these must be linked to ranlib.f)

#1 Simple averaging and calculation of error

#2 Weight function through change of variable

#3 Weight function through von Neumann rejection

Metropolis algorithm: here are 2 codes that compute via the Metropolis algorithm (see the recommended texts for details). Both codes use ranlib.f and metrolib.f.

Version 0: uses exp(-x**2) as a weight function. Here is the distribution of {x_i}.

Version 1: uses exp(-x**6) as a weight function. Here is the distribution of {x_i}.

Comparison of different runs.

The subroutine metropolis (in metrolib.f) writes the functions FUNC0 and FUNC1 to files, so that one can look at the correlation function C(k) = < f(x_i)f(x_i+k) > - <f>²

< x**2 exp(-x**2) > = < FUNC1 >

< exp(-x**2) > = < FUNC0>

For version 0, FUNC0 =1 and correlation function is trivial. Correlation function for FUNC1.

For version 1, correlation functions for both FUNC1 and FUNC0

(The program I used to compute the correlation function is here, but it is not well documented)

More notes on the Metropolis algorithm

Another simple Monte Carlo example -- computing volume of n-dimensional sphere

Linear algebra

Matrix multiplication (for benchmarking speed and optimization) requires library ranlib.f

Inversion of a matrix via LU decomposition (requires the library matinv.f (which includes needed LU decomposition and substitution routines))

Eigenvalues

Finding eigenvalues (requires the library eiglib.f )

Benchmark: Comparison of solving the same tridiagonal matrix through (1) full Householder and (2) through just calling TQLI (both require eiglib.f and ranlib.f )

Ordinary differential equations

Solving a simple, first-order DE:

Explicit Euler

Explicit Euler + Explicit Taylor series comparison

Explicit Euler + Explicit Taylor series + implicit comparison

Example: baseball trajectory + air resistance (Euler vs Taylor) comparison detailed comparison

Example: stellar profile for polytrope eqn of state (Euler only)

Figure: density vs R another figure

Figure: M vs R for different polytrope indices

Numerov solution of poisson eqn-- integrating inward

output for different box size L for L = 20, different dr

Numerov solution of Poisson eqn-- integrating outward