PHYS 580 COMPUTATIONAL PHYSICS SPRING 2005

PHYS 580 COMPUTATIONAL PHYSICS FALL 2009

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”).

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

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