Brandon Ling

Computational Physics


I once took a class on computational physics in grad school, but the class was taught using Fortran 90, which is usually used for intense numerical applications in academia and maybe industry, but not much else. I wanted to learn some Python and also review and learn more about computational techniques, so I decided to work my way through Mark Newman’s Computational Physics. My solutions to the exercises can be found here.


The most interesting thing I learned was the idea of Richardson extrapolation. The book introduces this through Romberg integration and the Bulirsch-Stoer method of solving ordinary differential equations.

Here’s the basic idea. Suppose you have a method to approximate some quantity AA, and that this method depends on a small parameter hh in such a way that the approximations get better and better the smaller hh is. As a concrete example, think of computing a definite integral and hh as the step size in a numerical integration method such as the trapezoidal rule.

The idea of Richardson extrapolation is that by approximating AA with two different (appropriately related) values of hh, you can combine them to get a new estimate with a much smaller error estimate. How those values of hh should be related and how you combine them depend on the specifics of what you’re computing.

For example, in Romberg integration you compute an estimate I1I_1 with step size hh, then compute a second estimate I2I_2 with a step size h/2h/2. You get a new, better estimate by adding to I2I_2 a term proportional to I2I1I_2 - I_1.


Adaptive Bulirsch-Stoer

Oscillating chemical reaction
An oscillating chemical reaction modeled by a pair of coupled ODEs.

In the example above, we solve the coupled, nonlinear ordinary differential equations

dxdt=1(b+1)x+ax2y,dydt=bxax2y\frac{dx}{dt} = 1 - (b + 1) x + a x^2 y, \quad \frac{dy}{dt} = bx - a x^2 y

with a=1,b=3a=1, b=3 with initial conditions x(0)=y(0)=0x(0) = y(0) = 0 and accuracy of 101010^{-10} per unit time. We solved this system using an adaptive Bulirsch-Stoer method. See my solution here.

Adaptive trapezoidal rule

Adaptive trapezoidal rule example
The dots represent the endpoints of the integration slices used in the trapezoidal rule. Notice that there are fewer slices in regions where f(x) changes less rapidly.

Another example is using an adaptive trapezoidal rule to estimate the value of

010(sinxx)2 dx\int_0^{10} \left(\frac{\sin x}{x}\right)^2 \ dx

We use the trapezoidal rule to compute the integral, but the step size is automatically adjusted so that the error estimate for each integration slice is no greater than some target value (so that the overall error is bounded). See my solution here.

2D Ising model

Of course we also simulate the well-known 2D ferromagnetic Ising model. We perform a Markov chain Monte Carlo simulation of the Ising model on a 100 × 100 square lattice with zero magnetic field and nearest-neighbor interactions.

system of spins on 100 by 100 lattice
The system after ten million simulation steps. Note we are near the critical temperature. The blue stands for a down spin, and orange stands for an up spin.
plot of magnetization
The average magnetization per spin of the system during each step of the simulation. Note we are near the critical temperature.