__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 (in Fortran)

Same
program
but somewhat better written (in Fortran)
(in C)

Compare 2-point and
3-point differentiation as a function of x (in Fortran)
(in
C) Graph
for
dx = 0.3 for different x

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

Another code comparing 2-point
and 3-point derivatives (in Fortran)*
*(in
C)

*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 (in fortran)

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

(to compile, do
gfortran
–o qtest.x quadmaster.f90 quadraturesubs.f90 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)

Alternate
version using standard fortran function *random_number*

For C/C++ use rand() in <stdlib.h>

*Exercise*: modify the
example calling routine above but instead of computing the average and
std deviation of a uniform distribution (ran3 or random_number) 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.

#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, also these notes
on the Metropolis algorithm).

Both codes use
ranlib.f and metrolib.f90.

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}}.

The subroutine metropolis
(in metrolib.f90) 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>^{2}

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

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

**Introduction to
parallel computing.
**

**Linear algebra**

Matrix multiplication (for benchmarking speed and optimization) uses Fortran inbuilt function random_number()

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

LU decomposition in C (will also need nrutil.h )
: ludcmp.c and lubksb.c

Notes on "nearly painless conjugate gradient methods"

Notes on Broyden's method by Dan Mickelsen;
for an example application see Phys
Rev. C **78**, 014318 (2008).

For a derivative-free minimizer POUNDerS see section
III in Phys. Rev C **82**, 024313 (2010).

**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 )

Here are the same codes in C (actually C++ despite the suffix): speedeigen.c tred2.c tqli.c ranlib.c nrutil.h

**Ordinary differential equations**

Solving a simple, first-order DE:

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

Outbound solutions for different dr ; difference between outbound soln and exact (homogeneous soln)

2nd order Runge-Kutta and comparison with 1st order Euler Graph of results

4th order Runge-Kutta for 2 dependent variables: rungekuttalib.f ; application to simple differential equation and resulting graph

More 4th-order Runge-Kutta: a central
force
problem (must link to rungekuttalib.f)
which calculates orbit (r vs. theta, or x vs y). Some results: (1)
"gravitational" force (1/r^{2}
force) with different angular momenta; (2) gravitational
force with different step sizes in theta; (3) different
central forces (1/r^{n}) notice the perihelion shift.

**Partial differential equations**

Sample code for Neumann boundary condition du/dx = 0 at x = 0, L

Example of naive discretization: Forward time, centered space (FTCS) and sample output showing instabilities: example 1 and example 2

Lax solution via numerical dissapation/viscosity, and example output

Implicit inversion and example output

Crank-Nicholson "unitary" time evolution and example output

Relaxation for simple inhomogenous differential equations:

Dirichlet boundary conditions with simple relaxation sweep

Graph of output and with "overrelaxation"

Combined Dirichlet and Neumann boundary conditions (with improved relaxation sweep) and output graph

**Root finding**

Newton-Raphson root-finding for two functions:
x^{5} -5 and tanh(x)
. Newton-Raphson converges for the former but diverges for the latter.

*Last modified Nov 15, 2017
*