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


Simple program for differentiation (in Fortran)             (version in C)

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)    (in C)

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 4th 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, 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 {xi}.

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


Comparison of different runs.


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(xi)f(xi+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))

Sample matrix for above

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


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 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/r2 force) with different angular momenta; (2) gravitational force with different step sizes in theta;  (3) different central forces (1/rn) 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

Improved 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: x5 -5  and tanh(x) . Newton-Raphson converges for the former but diverges for the latter.

Last modified Nov 15, 2017