Quantum Tunneling, Part 1

In this Python project, we will look at quantum tunneling. Even if you are completely new to quantum mechanics, you must have heard of this mysterious, quantum property of matter. Particles can, apparently, move through solid barriers! What is this quantum woo all about? This, we will find out by solving a standard textbook problem — that of a one-dimensional potential barrier. Our goal will be to derive the tunneling probability. We will then implement a simulation of a particle near such a potential barrier and look at the numerical solution of the problem.

The Schrödinger equation

What do we do when we want to study the evolution of a macroscopic physical system — a bullet, a ping-pong ball, the Solar system — in time? We reach for an equation of motion, such as the second law of Newton. To describe a system in classical mechanics we only need the equation of motion and two pieces of data: the initial position, and the initial velocity of all parts of the system. Armed with these, we can predict the state of the system at any arbitrary point in time.

Quantum mechanics differs from classical mechanics in the equation of motion and the required initial conditions. A quantum particle is described not by its position and velocity, but rather by a complex function — called the wave function — spread out over all of space. Its motion (or change), in turn, is dictated by Schrödinger’s equation of motion:

\hat{H}\Psi(\mathbf{r}, t) = i \hbar \frac{\partial}{\partial\,t} \Psi

The capital H with the fancy hat is the Hamiltonian operator. It is related to the Hamiltonian we discussed some time ago in relation to the classical Hamiltonian description of the double pendulum. The difference here is that instead of just a scalar number (which can be interpreted as the total energy of the system in certain cases), it is an operator — a mathematical entity that acts on a wave function and changes it, producing a different wave function.

Furthermore, the above equation is clearly a first-order partial differential equation in time. Hence, only the initial state of the system is required as the initial condition.

Finally, it would be nice to know the shape of the Hamiltonian operator. Generally, one can go from classical mechanics to quantum mechanics just by putting hats on the quantities (yes, it’s often that easy). For example, a classical particle in a spatially changing potential (think electric potential energy or gravitational potential) has the following Hamiltonian:

H = \frac{\mathbf{p}^2}{2m} + V(\mathbf{r})

The corresponding quantum-mechanical Hamiltonian operator is:

\hat{H} = \frac{\hat{\mathbf{p}}^2}{2m} + \hat{V}(\mathbf{r})

Now we just plug in the expressions of the position and momentum operators in position space (read regular space): \hat{\mathbf{r}} = r and \hat{\mathbf{p}} = -i \hbar \nabla. We find:

\hat{H} = - \frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r})

Our Schrödinger equation therefore looks like this:

\left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r})\right]\Psi = i \hbar \frac{\partial}{\partial\,t} \Psi

What are its solutions \Psi? Since the left hand side of the equation is independent of time, we could try to look for a separated solution. Moreover, the right hand side would then represent a simple first-order differential equation in time. It makes sense, then, to look for a solution in the form:

\Psi(\mathbf{r}, t) = \psi(\mathbf{r})e^{-i\frac{E}{\hbar}t}

Plugging this Ansatz into our equation we have that:

e^{-i\frac{E}{\hbar}t}\left[\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r})\right]\psi(\mathbf{r}) = i \hbar \psi(\mathbf{r})  \frac{\partial}{\partial\,t} e^{-i\frac{E}{\hbar}t}

Carrying out the time derivation and dividing by the complex exponential we get:

\left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r})\right]\psi(\mathbf{r}) = E \psi(\mathbf{r})

This is called the time independent Schrödinger equation. Its solutions are called eigenstates — they represent wave functions that, upon action of the Hamiltonian operator, remain identical up to a scalar multiplicative coefficient: the eigenvalue. This scalar value can in turn be identified with the total energy corresponding to the state. In simpler terms, think of the eigenstates as of the basic building blocks of the solution, while the eigenvalues are the corresponding energies. Let us now solve the Schrödinger equation for a model system illustrating the essence of quantum tunneling!

Rectangular Potential Barrier

Consider a one-dimensional potential barrier described by a potential V(x):

V(x) = \begin{cases} 0\ &\text{if}\ x < 0 \\ V_0\ &\text{if}\ 0 \leq x \leq a \\ 0\ &\text{if}\ x > a\end{cases}

This naturally splits the physical system into three parts. Let us designate the solution for all x < 0 by \psi_{\text{I}}(x), for 0 \leq x \leq a by \psi_{\text{II}}(x), and for all even greater x by \psi_{\text{III}}(x).

The time independent Schrödinger equation for the region \text{I} is:

-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\psi_{\text{I}}(x) = E \psi_{\text{I}}(x)

the solution of which is trivial. In the absence of a potential, the solution is a linear combination of sines and cosines. However, we might instead prefer the complex exponential form:

\psi_{\text{I}}(x) = A_r e^{ikx} + A_l e^{-ikx}

Plugging this into the left hand side of our equation of motion we find that:

-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\left( A_r e^{ikx} + A_l e^{-ikx} \right) = -\frac{\hbar^2}{2m} \left(-k^2A_re^{ikx} - k^2 A_l e^{-ikx}\right) = \frac{\hbar^2k^2}{2m} \psi_{\text{I}}(x)

and thus:

\frac{\hbar^2k^2}{2m} \psi_{\text{I}}(x)  = E\psi_{\text{I}}(x) \quad \Rightarrow \quad k = \sqrt{\frac{2mE}{\hbar^2}}


The solution in region \text{III} is, naturally, very similar:

\psi_{\text{III}}(x) = C_r e^{ikx} + C_l e^{-ikx}


Finally, we can look at the region of space inside the potential barrier. While the potential here is still constant, it is non-zero. This, however, only results in a change of the wavenumber; the solution itself is mathematically the same as in the other regions. We have that:

\psi_{\text{II}}(x) = B_r e^{iqx} + B_l e^{-iqx}

where:

q = \sqrt{\frac{2m(E - V_0)}{\hbar^2}}


How can we interpret the above solution? In the regions to both sides of the potential barrier, the solution represents a superposition of harmonic waves travelling in either direction. Inside the potential barrier, the solution depends on the energy of the particle compared to the barrier height. If E > V_0, the particle “flies above” the barrier. This has a direct analogy in classical physics, like a ball flying above a physical wall. However, if E < V_0, the wavenumber q is imaginary, and the solution becomes a superposition of real exponentials. While exponentials are quickly decreasing functions, they still lead to a non-zero wave function (and thus probability density) inside the potential barrier! Particles flying through walls!


The solution presented just now includes several free parameters. In total, there are seven of them: A_r,\, A_l,\, B_r,\, B_l,\, C_r,\, C_l, and the energy of the particle E. On the other hand side, the wave function has to satisfy certain boundary conditions. For one thing, it must be a smoothly varying function of the spatial coordinate. Hence, we get two continuity conditions on either side of the barrier. For another, the derivative of a wave function must also be continuous. Both the continuity of the wave function, and that of its derivative, stem from the mathematical structure of the Schrödinger equation — it is a second order differential equation in space. As long as the potential is finite, both continuity conditions must hold. Long story short, this gives us another two boundary conditions. Finally, we have the normalization condition. Unlike in the case of spatially confined systems, where the integral of the wave function is a finite number we can normalize, here we must use a different normalization criterion. However, the exact choice does not matter. The important point is that we get a fifth constraint.

This leaves us with 7 parameters and 5 constraints. Hence, the energy spectrum of the solution is continuous and doubly degenerate. In other words: for any choice of energy E, there exist two solutions, representing a particle moving in either direction.


Now, let’s see where we stand. The solution to the problem is:

\begin{aligned}\psi_{\text{I}}(x) &= A_r e^{ikx} + A_l e^{-ikx} \\  \psi_{\text{II}}(x) &= B_r e^{iqx} + B_l e^{-iqx}  \\  \psi_{\text{III}}(x) &= C_r e^{ikx} + C_l e^{-ikx} \end{aligned}

The continuity of the wave function and of its derivative yields four boundary conditions:

\begin{aligned} \psi_{\text{I}}(0) &= \psi_{\text{II}}(0)   \\ \left. \frac{\partial\psi_{\text{I}}}{\partial x}\right|_0 &= \left.  \frac{\partial \psi_{\text{II}}}{\partial x}\right|_0 \\ \psi_{\text{II}}(a) &= \psi_{\text{III}}(a) \\  \left. \frac{\partial \psi_{\text{II}} }{\partial x}\right|_a  &=  \left. \frac{\partial \psi_{\text{III}} }{\partial x}\right|_a    \end{aligned}

Plugging in the expressions for the wave functions \psi we get:

\begin{aligned}A_r  + A_l  &= B_r + B_l \\  i k \left(A_r - A_l\right) &= iq\left(B_r - B_l\right) \\ B_r e^{i q a} + B_l e^{-i q a} &= C_r e^{i k a} + C_l e^{-i k a} \\ iq\left(B_r e^{iqa} - B_l e^{-iqa}\right)&=  ik\left(C_r e^{ika} - C_l e^{-ika}\right)  \end{aligned}

This is the system of equations we need to solve. However, this represents a very general case.

Let us consider a particle moving in from the left with a unit amplitude A_r \equiv 1. We will neglect the case of a particle moving towards the barrier from the right by setting C_l \equiv 0. We now rename the parameter A_l \equiv r, as it represents the reflected amplitude, and C_r \equiv t, since this corresponds to the transmitted part. We now have:

\begin{aligned}1  + r &= B_r + B_l \\  i k \left(1 - r\right) &= iq\left(B_r - B_l\right) \\ B_r e^{i q a} + B_l e^{-i q a} &= t e^{i k a} \\ iq\left(B_r e^{iqa} - B_l e^{-iqa}\right)&=  ikt e^{ika}  \end{aligned}

Four equations, four parameters. We can get rid of B_r and B_l, and find the expressions for the reflected and transmitted amplitudes. Their squares, |r|^2 and |t|^2, give us the reflection and transmission coefficients, respectively.


As this post is getting quite long, let me stop here. In the next instalment we shall solve the above system of equations, find the two parameters r and t, and calculate the probability with which the particle tunnels through our potential barrier.

Until next time!

Advertisements

Radial Distribution Function

For my research, I study semiconductor / water interfaces, much like this one:

One of the methods I use is ab-initio molecular dynamics. The idea is rather simple: the nuclei are treated as classical particles, propagated in space and time using the Verlet algorithm. At each time step, the electronic structure of the current atomic configuration is calculated using density functional theory (DFT), which, in principle, should yield the exact electronic density of the ground state of the system. Next, the force on the nuclei is calculated taking the newly-computed distribution of electrons. And finally, the force is used to propagate the nuclei another step forward in time.

Due to the inherently high computational cost of first-principles methods such as DFT, we can afford to simulate only so many atoms. Careful tests have to be carried out to ensure that our simulation box is large enough to faithfully replicate the much larger system we are trying to model — at least from the perspective of the studied phenomenon. Among other things, we need to verify that the layer of water representing the solvent is sufficiently thick so as to allow for the water to behave as a liquid. Too few water molecules cramped between the periodically repeated interfaces, and the water remains frozen, bound at the interface.

In today’s post, we will write a script to investigate the radial distribution function of the oxygen atoms in water. This is a key indicator of the liquid-like quality of the simulated water layer. Let’s hop into it!

Radial Distribution Function

So what is this radial distribution function? As it often so happens, a picture is worth more than a thousand words:

The radial distribution function (RDF), pair correlation function, or often just g of r, describes the probability of finding a particle at a given distance from a reference particle. In the picture above, the oxygen-oxygen RDF for liquid water is shown. Imagine you are the oxygen atom of any given H2O molecule. If you observe your surroundings over some period of time, how far away, on average, are other water molecules (or their oxygen atoms, to be more specific)? At distances below 2.5 angstrom (0.25 nanometers), you are extremely unlikely to see any other oxygen atom. Water molecules simply do not want to be that close, the reason for the remarkable incompressibility of liquid water. At around 2.7 Å, a sharp peak in the RDF is observed. That’s the first solvation shell. While water molecules are strongly repulsive at short distances, a powerful network of hydrogen bonds binds them together. If you average out over each angle, and over the neighborhood of every water molecule, these hydrogen bonds manifest themselves in a sharp increase in the probability of finding water molecules at a distance of about 2.7 Å. Above this distance, the RDF drops, and then oscillates a bit, before finally converging to a value of 1. The interpretation is as follows. If \rho is the macroscopic mass density of water (1 g / cm), then \rho g(r) is the density observed at a distance r from the reference point. It makes sense that in a disordered system like liquid water, the local water density at large distances should approach the macroscopic value, hence g(r) \rightarrow 1.

The Algorithm

The goal here is simple. Given a molecular dynamics trajectory — a series of atomic configurations at each time step — compute the oxygen-oxygen RDF. Comparison of the calculated RDF with that of liquid water will tell us how close (or not) we are to simulating the proper behaviour of the water layer. How do we go about it?

First, let’s start with an oxygen atom of our choice as the reference. We will iterate over every other water-molecule-oxygen atom and calculate the distance. We will then create a histogram of the oxygen atom count as a function of distance. Repeat this for every water oxygen atom, and average over every frame. This will give us the distribution of oxygen atoms as a function of their distance. Naturally, the greater the distance, the larger the volume of the spherical shell corresponding to the corresponding bin of the histogram. We need to divide the count of the oxygen atoms by the volume of these spherical shells to achieve something akin to a numerical density. From there, it is a matter of simple normalization to arrive at the final radial distribution function. The following series of graphics illustrates the process:

There are two points to keep in mind. First of all, the simulation cell is periodically repeated. Hence, when evaluating distances between two atoms, we have to take the cell parameters into account. Two atoms at the opposing ends of a simulation box are, technically, quite close to each other (or rather, close to their respective periodic image). Second, when calculating the volume of the spherical shell corresponding to a given distance, we need to account for the finite size of the water region. Since water molecules are present in a portion of the simulation box only, we should just consider the volume of the part of the spherical shell that falls within this subspace.

The Script

Let’s start simple with a function to calculate the distance between two particles represented by arrays of three coordinates:

import matplotlib.pyplot as plt 
import scipy as sp

def distance(a, b):
    dx = abs(a[0] - b[0])
    x = min(dx, abs(A - dx))
    
    dy = abs(a[1] - b[1])
    y = min(dy, abs(B - dy))
    
    dz = abs(a[2] - b[2])
    z = min(dz, abs(C - dz))

    return sp.sqrt(x**2 + y**2 + z**2)

This way, we correctly account for the periodic images of the atoms we sample. Note that A, B, and C are the cell parameters of the computational cell, which we assume is orthogonal. I know, I know, global variables are evil. Deal with it.


Next, we implement a function to calculate the volume of a sphere of a given radius. Recall, however, that we need to take the size of the water region into account. We will assume that the water layer is sandwiched between two horizontal semiconductor interfaces. Only a certain range of vertical coordinates is thus admissible. When we calculate the volume of the sphere, we only restrict ourselves to this range of the vertical coordinate. Anything below or above will get cut off. Go look up the volume of a spherical cap to understand the calculations in this one:

def volume(z, r, z_bot, z_top):
    """ volume of a sphere of radius r located at height z """
    volume = 4.0 / 3.0 * sp.pi * r**3
    
    """ cut off a spherical cap if the sphere extends below z_bot """
    if z - r < z_bot:
        h = z_bot - (z - r)
        volume -= sp.pi * h**2 * (r - h / 3.0)

    """ cut off a spherical cap if the sphere extends above z_top """
    if z + r > z_top:
        h = z + r - z_top
        volume -= sp.pi * h**2 * (r - h / 3.0)
    
    return volume


With these two helper functions down, we go on to implement the Trajectory class:

class Trajectory:
    def __init__(self, filename, skip, z_bot_interface, z_top_interface,
                 interface_offset=4.0, resolution=200):

        self.filename = filename
        self.skip = skip
        self.z_top = z_top_interface - interface_offset
        self.z_bot = z_bot_interface + interface_offset
        self.resolution = resolution
        
        self.parse_input()
        self.compute_volume_per_h2o()

    ...

The input file, passed to our class through the filename argument, is expected to be of the .xyz format. It’s rather simple. The .xyz coordinate file is a human-readable file format. The first line holds the number of atoms. The second line is reserved for comments. Then, for each atom, there is one line that consists of the chemical symbol of the atom type, and three floating point numbers representing the x, y, and z coordinate of the atom in angstrom. Knowing this, we will implement the parse_input() method as follows:

class Trajectory:
    ...
    def parse_input(self):
        with open(self.filename, 'r') as f:
            data = f.readlines()

        self.n_atoms = int(data[0].split()[0])
        self.n_steps_total = int(len(data) / (self.n_atoms + 2))
        self.n_steps = self.n_steps_total // self.skip
        
        self.atom_list = [line.split()[0] for line in 
                          data[2 : self.n_atoms + 2]]

        self.coordinates = sp.zeros((self.n_steps, self.n_atoms, 3))
        for step in range(self.n_steps):
            coords = sp.zeros((self.n_atoms, 3))
            
            i = step * self.skip * (self.n_atoms + 2)
            for j, line in enumerate(data[i + 2 : i + self.n_atoms + 2]):
                coords[j, :] = [float(value) for value in line.split()[1:]]
            
            self.coordinates[step] = coords

    ...

Let’s go step by step here. First, we read the whole file into data. While this is approach is rather convenient, the input trajectories can be very large, on the order of hundreds of megabytes — even several gigabytes. In such a case, it would be much more efficient to parse the file line by line on the go. We are not going to worry about such things now.

Next, we read in the number of atoms. We do that by taking the first line of the trajectory, splitting it into words (to get rid of the end-of-line character), taking the first word, and converting it to an integer. We then calculate the number of steps by dividing the length of the file by the number of atoms, not forgetting the 2 header lines per each configuration. Finally, we define the n_steps member variable as a subset of the total number of steps. Since the the distribution of water molecules does not change from one step to another (remember, we need a short time step to capture the oscillations of the light hydrogen nuclei), we will get a perfectly neat RDF by considering only every 10th, or even every 100th snapshot. If we need full accuracy and / or statistics, we can decrease skip to 1.

Finally, we parse the actual data itself. We need to read the chemical symbols of the atoms just once. As for the coordinates, we create a member variable in the form of a three-dimensional array holding the x, y, and z coordinates of each atom at each step of our subset of the trajectory. Some list comprehension magicks are happening here, but by this point you should be able to comprehend what is going on.


With the coordinates parsed, we can go on to calculate an important quantity: the volume (in angstrom) per water molecule. We will use this number to normalize our RDF such that it approaches a value of one for large distances. Moreover, we can compare the computed value with the experimental one to determine how close our simulated water is to the real deal:

class Trajectory:
    ...
    
    def compute_volume_per_h2o(self):
        n_h2o = 0
        for step in range(self.n_steps):
            for i, atom in enumerate(self.coordinates[step]):
                z = atom[2]
                if self.atom_list[i] == "O" and self.z_bot < z < self.z_top:
                    n_h2o += 1
        
        volume = A * B * (self.z_top - self.z_bot)
        average_n_h2o = n_h2o / self.n_steps
        self.volume_per_h2o = volume / average_n_h2o

    ...

We simply loop over every snapshot and count oxygen atoms that lie within our water region. By the way — did you know you can do comparisons of the type a < x < b in Python? Now you know.


Now that the preliminaries are sorted out, we can implement the algorithm to compute the oxygen-oxygen RDF described above:

class Trajectory:
    ...
    
    def compute_radial_distribution(self):
        """ no reason to go above half of the smallest lattice parameter
            as mirror images start to be double-counted """
        r_cutoff = min(A, B) / 2.0

        dr = r_cutoff / self.resolution
        self.radii = sp.linspace(0.0, self.resolution * dr, self.resolution)

        volumes = sp.zeros(self.resolution)
        self.g_of_r = sp.zeros(self.resolution)
        
        for step in range(self.n_steps):
            print('{:4d} : {:4d}'.format(step, self.n_steps))
            
            """ isolate all liquid water molecules based on their oxygens """
            data_oxygen = []
            for i, atom in enumerate(self.coordinates[step]):
                if self.atom_list[i] == "O":
                    if self.z_bot < atom[2] < self.z_top:
                        data_oxygen.append(atom)
            data_oxygen = sp.array(data_oxygen)
            
            for i, oxygen1 in enumerate(data_oxygen):
                """ calculate the volume of the cut spherical shells """
                for j in range(self.resolution):
                    r1 = j * dr
                    r2 = r1 + dr
                    v1 = volume(oxygen1[2], r1, self.z_bot, self.z_top)
                    v2 = volume(oxygen1[2], r2, self.z_bot, self.z_top)
                    volumes[j] += v2 - v1

                """ loop over pairs of O atoms, each pair counts as 2 """
                for oxygen2 in data_oxygen[i:]:
                    dist = distance(oxygen1, oxygen2)
                    index = int(dist / dr)
                    if 0 < index < self.resolution:
                        self.g_of_r[index] += 2.0
        
        for i, value in enumerate(self.g_of_r):
            self.g_of_r[i] = value * self.volume_per_h2o / volumes[i]

    ...

First, we set to zero arrays representing the volumes of the spherical shells, the horizontal axis of the RDF, and the RDF itself. We only calculate the RDF up to a cutoff distance. We could consider distances as large as the lattice parameters. Above this value, pairs of periodic images would pollute the RDF and lead to nonsensical artifacts. However, already a distance larger than one half of the smallest lattice parameter leads to the unphysical double-counting of periodic images of some oxygen atoms. We will thus restrict ourselves to distances that correspond to relevant data, and choose the conservative cutoff of one half the lattice parameter.

Next, we loop over all oxygen atoms and isolate those that lie within the sampled bulk water layer. For each of these oxygen atoms we calculate the volumes of the (cut) spherical shells corresponding to all sampled distances.

We then loop over every oxygen atom pair, calculate their distance, and build up a histogram of the count vs. distance. We can speed up the algorithm by only looping over unique pairs, and counting each pair as two.

Finally, normalizing by the volumes of the spherical shells, and the volume per water molecule computed earlier, we arrive at the radial distribution function!


Now we just write a quick method to plot the result:

class Trajectory:
    ...

    def plot(self, filename=""):
        """ plots the radial distribution function
        if filename is specified, prints it to file as a png """
        
        if not self.g_of_r.any():
            print('compute the radial distribution function first\n')
            return
        
        plt.xlabel('r (Å)')
        plt.ylabel('g$_{OO}$(r)')
        plt.plot(self.radii, self.g_of_r)
        
        if filename:
            plt.savefig('radial_distribution_function.png', dpi=300,
                        bbox='tight', format='png')
        
        plt.show()

… and run the script like this:

A = 12.991
B = 11.832
C = 30.577
bottom_interface = 23.0
top_interface = 40.0

TiO2_H2O = Trajectory('TiO2_H2O-pos-1.xyz', 10, 
                      bottom_interface, top_interface)
TiO2_H2O.compute_radial_distribution()
TiO2_H2O.plot()

Et voilà:

Sure, the curve is a bit noisy. I used only around a thousand snapshots of a rather short molecular dynamics trajectory. More data means a smoother curve — hopefully approaching the oxygen pair distribution of liquid water even closer. But it’s apparent that the main features of the simulated water layer and those of bulk liquid water are very similar: the principal peak is at the correct distance and of the same height, and even the long range behavior is near-identical. Success!!

Until next time!

Double Pendulum, Part 3

In the previous instalment, we slogged through the derivation of Hamilton’s equations of motion for the double pendulum. We did this because these equations are a set of first order differential equations in a form that is very conducive to a numerical solution. In this part, we will adopt a few simplifications, write a Python script to simulate the system, and animate the motion of the double pendulum. Let’s go!


As we have found last time, the system is described by a set of four first order differential equations, two per each degree of freedom:

\dot{\theta}_1 = \dfrac{l_2p_1 - l_1p_2\cos(\theta_1-\theta_2)}{l_1^2l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

\dot{\theta}_2 = \dfrac{l_1\left(m_1 + m_2\right)p_2 - l_2m_2p_1\cos(\theta_1-\theta_2)}{l_1l_2^2m_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

\dot{p}_1 = -\left(m_1 + m_2\right)gl_1\sin(\theta_1) - A + B

\dot{p}_2 = -m_2gl_2\sin(\theta_2) + A - B

where

A \equiv \dfrac{p_1p_2\sin(\theta_1 - \theta_2)}{l_1l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

B \equiv \dfrac{l_2^2m_2p_1^2 + l_1^2\left(m_1 + m_2\right)p_2^2 - l_1l_2m_2p_1p_2\cos(\theta_1 - \theta_2)}{2l_1^2l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]^2}\sin\left[2(\theta_1 - \theta_2)\right]

While we could directly simulate the above equations, we will not. Let us do some assumptions about the lengths and masses of the pendulums in order to simplify the system. The implementation of the full, general case will be left as an exercise for the motivated reader (I used to hate this sentence with a passion, but I am slowly starting to see the point :).

We will assume that the lengths of the two pendulums are identical and equal to one, l_1 = l_2 = 1. We will also set the masses to m_1 = m_2 = 1. This greatly simplifies our equations:

\dot{\theta}_1 = \dfrac{p_1 - p_2\cos(\theta_1-\theta_2)}{1 + \sin^2(\theta_1 - \theta_2)}

\dot{\theta}_2 = \dfrac{2p_2 - p_1\cos(\theta_1-\theta_2)}{1 + \sin^2(\theta_1 - \theta_2)}

\dot{p}_1 = -2g\sin(\theta_1) -  \dfrac{p_1p_2\sin(\theta_1 - \theta_2)}{1 + \sin^2(\theta_1 - \theta_2)}  +   \dfrac{p_1^2 + 2p_2^2 - p_1p_2\cos(\theta_1 - \theta_2)}{2\left[1 + \sin^2(\theta_1 - \theta_2)\right]^2}\sin\left[2(\theta_1 - \theta_2)\right]

\dot{p}_2 = -g\sin(\theta_2) +   \dfrac{p_1p_2\sin(\theta_1 - \theta_2)}{1 + \sin^2(\theta_1 - \theta_2)}  - \dfrac{p_1^2 + 2p_2^2 - p_1p_2\cos(\theta_1 - \theta_2)}{2\left[1 + \sin^2(\theta_1 - \theta_2)\right]^2}\sin\left[2(\theta_1 - \theta_2)\right]

There, much better.

Python Script

We will now write a Python program that will solve the above equations numerically and plot the motion of the double pendulum. I will be using Python 3, but the code can be adapted for Python 2 with only minor changes.

Since our pendulum, as a physical system, is associated with both some data (parameters such as the lengths and masses) and some behavior, it makes sense to define a Pendulum class:

import matplotlib
matplotlib.use('Qt4Agg') # 'tkAgg' if Qt not present 
import matplotlib.pyplot as plt 
import scipy as sp
import matplotlib.animation as animation

class Pendulum:
    def __init__(self):
        ...

What parameters do we need to supply in order to initialize the pendulum? From a physical point of view, there are four first-order differential equations to be solved, which means we require four initial conditions. However, we may choose to always start from a state in which the initial velocities are zero. This would correspond to picking up the pendulum, and releasing it without giving it any nudge. Hence, we will restrict ourselves to solely the initial values of the two angles \theta_1 and \theta_2 — or theta1 and theta2 — as we will refer to them in our code.

Another parameter that we need to supply is the length of the time step $latex\text{d}t$. The smaller the time step, the more accurate the numerical integration of the equations of motion will be. But we need a value large enough to simulate a meaningful portion of the trajectory in a given amount of real time.

We will also set the values of some other parameters, such as the gravitational constant g and the pendulum lengths. However, the former is not expected to change, and the latter we already decided to put equal to 1.

Finally, we will append our initial position of the pendulum masses to our trajectory. Note, however, that for the purpose of plotting, we should keep track of the trajectory in Cartesian coordinates. To that end, we will use the method polar_to_cartesian():

import matplotlib
matplotlib.use('Qt4Agg') # 'tkAgg' if Qt not present 
import matplotlib.pyplot as plt 
import scipy as sp
import matplotlib.animation as animation

class Pendulum:
    def __init__(self):
        self.theta1 = theta1
        self.theta2 = theta2
        
        self.p1 = 0.0
        self.p2 = 0.0
        
        self.dt = dt
        
        self.g = 9.81
        self.length = 1.0
        
        self.trajectory = [self.polar_to_cartesian()]

    def polar_to_cartesian(self):
        x1 =  self.length * sp.sin(self.theta1)        
        y1 = -self.length * sp.cos(self.theta1)
        
        x2 = x1 + self.length * sp.sin(self.theta2)
        y2 = y1 - self.length * sp.cos(self.theta2)
        
        return sp.array([[0.0, 0.0], [x1, y1], [x2, y2]])
    ...

The polar_to_cartesian() method takes the current position of the two masses in terms of angles and lengths, and transforms them into Cartesian coordinates (horizontal and vertical positions). The method returns an array of the coordinates corresponding to the pivot, the first, and the second mass of the pendulum, respectively.

Finally, we implement a method that solves the equations of motion and advances the state of the pendulum by the time step $latex\text{d}t$. While we could use more advanced integration techniques, the simplest implementation will work fine for our purposes, as neither performance nor accuracy are of importance here. Given a differential equation of the form:

\dfrac{\text{d}x}{\text{d}t} = f(x)

we can estimate the value of the variable x at a time $t + \text{d} t$ from the value at time t by simply approximating the value of the derivative by a fraction:

\dfrac{\text{d}x}{\text{d}t} \approx \dfrac{x(t +  \text{d} t) - x(t)}{ \text{d}t} = f(x)

from where:

x(t + \text{d}t) = x(t) +  \text{d} t\, \cdot \, f(x)

The generalization to the case of multiple variables is trivial, and a possible implementation of this procedure for our equations of motion may look like the following:

import matplotlib
matplotlib.use('Qt4Agg') # 'tkAgg' if Qt not present 
import matplotlib.pyplot as plt 
import scipy as sp
import matplotlib.animation as animation

class Pendulum:
    def __init__(self):
        self.theta1 = theta1
        self.theta2 = theta2
        
        self.p1 = 0.0
        self.p2 = 0.0
        
        self.dt = dt
        
        self.g = 9.81
        self.length = 1.0
        
        self.trajectory = [self.polar_to_cartesian()]

    def polar_to_cartesian(self):
        x1 =  self.length * sp.sin(self.theta1)        
        y1 = -self.length * sp.cos(self.theta1)
        
        x2 = x1 + self.length * sp.sin(self.theta2)
        y2 = y1 - self.length * sp.cos(self.theta2)
        
        return sp.array([[0.0, 0.0], [x1, y1], [x2, y2]])
    
    def evolve(self):
        theta1 = self.theta1
        theta2 = self.theta2
        p1 = self.p1
        p2 = self.p2
        g = self.g
        l = self.length

        expr1 = sp.cos(theta1 - theta2)
        expr2 = sp.sin(theta1 - theta2)
        expr3 = (1 + expr2**2)
        expr4 = p1 * p2 * expr2 / expr3
        expr5 = (p1**2 + 2 * p2**2 - p1 * p2 * expr1) \
        * sp.sin(2 * (theta1 - theta2)) / 2 / expr3**2
        expr6 = expr4 - expr5

        theta1 += self.dt * (p1 - p2 * expr1) / expr3
        theta2 += self.dt * (2 * p2 - p1 * expr1) / expr3
        p1 += self.dt * (-2 * g * l * sp.sin(theta1) - expr6)
        p2 += self.dt * (    -g * l * sp.sin(theta2) + expr6)

        new_position = self.polar_to_cartesian()
        self.trajectory.append(new_position)
        return new_position

Note that here we define the variables theta1, theta1, p1, … , as well as the expressions labeled expr1 through expr6 , for the sake of being concise only (and also to avoid writing self a bazillion times). Once we update the angles and the angular momenta, we append the new position — in Cartesian coordinates — to the trajectory. We also make the method return the new positions for reasons that will be apparent later.


Now that we can simulate the motion of our pendulum, it is time we create an animation thereof. As the code snippets above already gave away, we will use the FuncAnimation() function from the animation module of the matplotlib package.

We introduce a new class, Animator, the instance of which will be responsible for plotting the movement of a pendulum. We will implement the following methods:

  • __init__(), the constructor, which takes an instance of Pendulum as an argument, and sets up the parameters of the figure,
  • advance_time_step(), which does exactly what it says — advances the time step and updates the state of the pendulum,
  • update(), which is responsible for updating the data displayed in the next frame of the animation,
  • and finally, animate(), which does the actual animation by calling the FuncAnimation() function.

Let’s start with the constructor:

class Animator:
    def __init__(self, pendulum, draw_trace=False):
        self.pendulum = pendulum
        self.draw_trace = draw_trace
        self.time = 0.0

        # set up the figure
        self.fig, self.ax = plt.subplots()
        self.ax.set_ylim(-2.5, 2.5)
        self.ax.set_xlim(-2.5, 2.5)

        # prepare a text window for the timer
        self.time_text = self.ax.text(0.05, 0.95, '', 
            horizontalalignment='left', 
            verticalalignment='top', 
            transform=self.ax.transAxes)

        # initialize by plotting the last position of the trajectory
        self.line, = self.ax.plot(
            self.pendulum.trajectory[-1][:, 0], 
            self.pendulum.trajectory[-1][:, 1], 
            marker='o')
        
        # trace the whole trajectory of the second pendulum mass
        if self.draw_trace:
            self.trace, = self.ax.plot(
                [a[2, 0] for a in self.pendulum.trajectory],
                [a[2, 1] for a in self.pendulum.trajectory])
        
    ...

That’s a lot of code, so let me explain what we are doing line by line. First, we initialize the member variables, including the pendulum we want to simulate the motion of, and a flag representing the optional drawing of a trace. We also define a time variable. One could argue that the pendulum could track its own time instead of delegating it to the Animator class. However, I wanted to leave the pendulum class as simple as possible. We only keep track of time to display it in the animation, hence the Animator seemed like the most rational place for storing this value.

Next, we set up a figure and an axes object, which will hold our actual plot(s). Keeping in mind that the lengths of both pendulums are equal to one, we choose appropriate dimensions for the window.

We then create a text window in the upper left corner to display the time. Coincidentally, we also align the text to the top left corner of the text window.

The initial configuration is plotted as a simple line plot with round markers. Recall that we save the trajectory as a list of scipy arrays. There are three arrays, each holding the horizontal and vertical coordinate for the pivot, the first mass, and the second mass of the pendulum, respectively. The first argument of the plt.plot() function is a list of horizontal coordinates, hence we take the last entry in the trajectory (in case we start plotting a pendulum at some later time), and extract the horizontal coordinate of the pivot and the two masses. We do the same with the vertical coordinates for the second argument of the plotting function.

Finally, in case we want to highlight the trajectory of the pendulum by drawing a trace, we set up a line plot (no markers here) of the whole trajectory of the second pendulum mass. Here, those not familiar with Python’s list comprehensions might be baffled at first. But the possibly cryptic one-liner:

x_coords = [a[2, 0] for a in self.pendulum.trajectory]

can be rewritten as:

x_coords = list()
for step in self.pendulum.trajectory:
    second_pendulum = step[2]
    x = second_pendulum[0]
    x_coords.append(x)

List comprehension is a powerful and versatile tool that often allows us to condense code while retaining — or even improving — clarity. Many times it’s easier to comprehend a single line than it is to take in a whole for-loop with several temporary variables thrown into the mix.


The second method to implement is advance_time_step():

class Animator:
    ...
    def advance_time_step(self):
        while True:
            self.time += self.pendulum.dt
            yield self.pendulum.evolve()
    ...

I will come back to the yield statement later. We then implement a method to update the data in the plots:

class Animator:
    ...
    def update(self, data):
        self.time_text.set_text('Elapsed time: {:6.2f} s'.format(self.time))
	
        self.line.set_ydata(data[:, 1])
        self.line.set_xdata(data[:, 0])
        
        if self.draw_trace:
            self.trace.set_xdata([a[2, 0] for a in self.pendulum.trajectory])
            self.trace.set_ydata([a[2, 1] for a in self.pendulum.trajectory])
        return self.line,
    ...

which should be pretty self-explanatory.


Finally we tie everything together in the animate() method:

class Animator:
    ...
    def animate():
        self.animation = animation.FuncAnimation(self.fig, self.update,
                         self.advance_time_step, interval=25, blit=False)

Let’s discuss this FuncAnimation() function. The first argument is the figure on which the animation will draw each frame. Our figure comes with one, possibly two plots (depending on the presence or absence of the trace), and a text window.

The second argument is the function which updates the data in those plots. If we look at the update() method above, we see that it updates the data in the plots through data which is passed to this method as an argument.

The third argument is supposed to hold, or calculate and return, the actual values to be passed to the update function of the second argument. We, of course, use advance_time_step(). Why do we use the yield statement instead of the more recognizable return statement?

The reason is that the FuncAnimation() function expects an iterator. An iterator is an object which can be iterated upon. More specifically, it implements a next() method that returns consecutive elements of the object. In practice, this means that the object can be used in e.g. a for-loop. Now, the yield statement is used to turn a regular function (one that uses the return statement) into a generator. A generator function will calculate and a single return value when needed, and then wait until the next value is required. In effect, the function itself acts as an iterator!

Iterators and generators might seem quite abstract — and they are supposed to. In effect, iterators allow one to iterate over any custom data structure. Hence, those that develop functions do not have to worry about the data structures the user will try to use them with — whether they be lists, tuples, arrays, dictionaries, custom data structures, or other classes. As long as they act as iterators (have a correctly implemented next()), they will work fine.

And generators, in turn, make it possible to create iterators out of functions. Not only does this provide even more freedom when writing code, it has some practical benefits as well. For example, if you need to generate a very long list of values, but only ever need one at a time, creating a list of all values and storing it is needlessly slow, not to mention it could potentially require a lot of memory for storage. In contrast, a generator only yields a single value at a time, avoiding costly storage and calculation of extensive amounts of data.

Neat! The fourth argument is the time, in milliseconds, to display each frame; and the fifth argument refers to a technique, called blitting, which enables one to only redraw the parts of the figure that are changing between each frame, thus saving resources. If, for any reason, your system does not support blitting, pass the value blit=False.


Now that we are done with our two class definitions, let us write a quick script that instantiates a pendulum and an animator, and animates the movement of the pendulum:


import matplotlib
matplotlib.use('Qt4Agg') # 'tkAgg' if Qt not present 
import matplotlib.pyplot as plt 
import scipy as sp
import matplotlib.animation as animation

class Pendulum:
    ...

class Animator:
    ...

pendulum = Pendulum(theta1=sp.pi, theta2=sp.pi-0.01, dt=0.01)
animator = Animator(pendulum=pendulum, draw_trace=True)
animator.animate()
plt.show()

And the output:


Splendid! We started with a seemingly trivial problem of a double pendulum. We managed to derive the equations of motion for the two pendulum masses, both in the Lagrange and in the Hamiltonian formalism. We then wrote a Python program to integrate Hamilton’s equations of motion and simulated the movement of the pendulum. Mission accomplished!

Note the chaotic behavior of the system. The motion starts with the pendulum being oriented vertically. As it drops and starts flailing around, the trajectory looks kind of symmetric, but not completely. Not only did we start from an initial configuration that was slightly leaning to one side, we also introduce tiny errors at every time step. We could reduce \text{d}t arbitrarily, but due to the very nature of the system, all these tiny errors would still lead to an unpredictable, vastly different trajectory. That’s the hallmark of chaos. Small changes in initial conditions result in completely different solutions. I plan to talk about chaos in some more detail in a future series — so if you are interested in this topic, stay tuned!

Until next time!

Abelian Sandpile: a C++ Project, Part 3

In the previous part, we managed to implement the toppling of the tower of sand grains on any given tile. In today’s instalment, we will finish the algorithm that will repeatedly iterate over the whole grid until every single pile is toppled. A chapter will be dedicated to some modern C++ topics, such as r-values and move semantics. Finally, we will run our code and print out the resulting configuration — in plain ASCII only, for now.


Implementing the topple_all() Method

File “Sandpile.h”

#pragma once

class Sandpile
{
private:
    Grid <int> m_grid;
    int m_height;

public:
    Sandpile(int height);
    ~Sandpile();

    topple();
    topple_all();
    draw();
};

With the topple() method finished, we can now implement the topple_all() method responsible for repeatedly toppling each pile until no grid tile contains more than three grains of sand.

We introduce the boolean variable is_larger_than_four which keeps track of exactly what it says: initialized to false, it will become true if any pile in the grid is more than three grains high. We will iterate over all grid tiles, counting the number of sand grains. If we ever find more than four, we switch our boolean to true and topple the pile by calling the topple() method. Note that we are not done even after the whole grid is traversed, as there may be piles that reached a height of more than three grains due to their neighbours toppling. Hence, as long as the boolean is true, we must keep iterating over the grid. We will switch
is_larger_than_four to false prior to each loop. This way, the algorithm stops once there are no more piles to topple. Let’s code it up!

File “Sandpile.cpp”

void Sandpile::topple_all()
{
    bool is_larger_than_four = true;
    while (is_larger_than_four)
    {
        is_larger_than_four = false;
        
        for (int i = 1; i < m_size - 1; ++i)
        {
            for (int j = 1; j < m_size - 1; ++j)
            {
                int height = m_grid.at(i, j);
                if (height > 3)
                {
                    is_larger_than_four = true;
                    topple(i, j, height);
                }
            }
        }
    }
    return;
}

Let us remark that we iterate starting at the second tile in each row (and column), and finish at the penultimate grid element. The boundary of our simulation grid acts as a sort of buffer enabling us to ignore boundary conditions for the most part. More on this has been written in the previous part of the series.

Constructor (and Destructor)

We will implement a constructor that takes one argument of type integer, representing the height of the initial pile placed at the center of our grid. We also need to choose the size of the simulated grid. One could go for a fixed number, like 1000 times 1000, but what if we want to just produce a small picture? Or what if we use too many grains of sand and the toppling piles “overflow” the boundary of the grid?

It therefore seems wiser to estimate the required size of the grid using the initial height the sandpile. We know that there are only four possible states for each tile in the final configuration: it’s either empty, or holds 1, 2, or 3 grains. Assuming a uniform distribution of the number of grains in occupied tiles, we find a value of 2.0 (one could argue that the average should be 1.5, but strangely enough, there are not that many empty tiles “inside” the occupied region, skewing the average number of sand grains towards a higher value). The images I included in the first part show a rather circular pattern. So, if we put N grains on a grid and let it topple, and the final configuration is close to a circle with 1.5 sand grains on average, we have:

N = 2\pi a^2 \quad \Rightarrow \quad a = \sqrt{\dfrac{N}{2\pi}}

Let us not forget that we only simulate one quadrant, which means that a square grid of edge length equal to the radius of the circular image is enough. However, we also need to have two extra rows and columns as boundaries, and we may include some extra padding just to be safe. Hence, an estimate for the grid size a as a function of the number of grains N is:

a =   \sqrt{\dfrac{N}{2\pi}} + 4

Using this estimate, our constructor will look something like this:

File “Sandpile.cpp”

Sandpile::Sandpile(int height) : m_height(height)
{
    m_size = (int)std::sqrt(height / 2.0 / 3.14159265359) + 4;
    m_grid = Grid<int>(m_size, m_size);
    m_grid.at(1, 1) = height;
}

This won’t compile as is — the assignment operator is not implemented for the Grid class. We could easily write one, but we would encounter a problem. We can’t just take one instance of the Grid class and copy its members to another. The actual grid data is allocated on the heap and a std::unique_ptr points to it. And a std::unique_ptr can not be copied! That’s the whole point of having unique pointers: they may not be randomly copied and passed around.

Hence, instead of copy assignment we must implement move assignment. Move assignment may be necessary in this case, but it is often desirable as well. It avoids allocating new memory, copying content, and deallocating old memory by simply passing the whole resource from one instance to another. So, how does move assignment look like in practice?

File “Grid.h”

#pragma once
#include 

template 
class Grid
{
private:
    std::unique_ptr<T[]> m_grid{ nullptr };
    int m_size_x{ 0 };
    int m_size_y{ 0 };

public:
    Grid<T> & operator=(Grid<T> && other)
    {
        if (this != &other)
        {
            m_size_x = other.m_size_x;
            m_size_y = other.m_size_y;
            m_grid = std::move(other.m_grid);
        }
    return *this;
    }

    ...
};

First, the return type is a reference to a Grid. This makes it possible to chain the assignment operator multiple times in a single statement. Honestly, we will probably never do such a thing, but that’s how the assignment operator is commonly overloaded. Second, notice the use of two consecutive ampersands. This is another modern C++ concept: an r-value reference. Without going into too much detail, one can split variables into l-values and r-values. Roughly speaking, while l-values refer to a specific place somewhere in memory, r-values don’t. Literals are examples of r-values — if you write int x = 10, a place in memory is assigned the handle a, but the literal 10 doesn’t have a location in memory, it only exists as some transient value in some registers on the CPU.

When we tell the compiler to execute the command:

m_grid = Grid<int>(m_size, m_size);

we don’t really need a Grid to be created somewhere in memory just so that its contents are a moment later moved to m_grid. So we use two ampersands to pass it as an r-value reference to the assignment operator. Next, we avoid self-assignment by verifying that both operands are different instances. We then set the grid dimensions by simply copying them. We use std::move to move the std::unique_ptr without copying the grid data between the two instances. The left hand side instance simply transfers its pointer to the right hand side one. Finally, we return the now updated instance by returning the dereferenced this pointer.


Whenever we implement move assignment, it is good practice to implement the move constructor as well, and vice versa. I will leave this as an exercise for the motivated reader, as the implementation of the move constructor looks very similar to the above assignment operator.

Printing the Results

Before we get to write the draw() method, let us write a quick print() function that iterates over the toppled sandpile and prints out the amount of sand grains remaining on each grid element. That way we can finally test whether our algorithm works!

File “Sandpile.h”

#pragma once

class Sandpile
{
private:
    Grid <int> m_grid;
    int m_height;

public:
    Sandpile(int height);
    ~Sandpile();

    topple();
    topple_all();
    print();
    draw();
};

File “Sandpile.cpp”

void Sandpile::print()
{
    for (int i = 1; i < m_size - 1; ++i)
    {
        for (int j = 1; j < m_size - 1; ++j)
        {
            int height = m_grid.at(i, j);
            std::cout << " " << height;
        }
        std::cout << std::endl;
    }
    return;
}


We have everything set up to run a simulation of the Abelian sandpile. Let’s modify our main() function to create a pile containing a given number of sand grains at the origin, topple it, and print the resulting configuration!

File “Abelian_Sandpile.cpp”

#include "Sandpile.h"

int main()
{
    Sandpile pile(40);
    pile.topple_all();
    pile.print();

    return 0;
}

Let’s go!

That’s it. Well, I don’t know what you expected, but it works. We don’t mirror the simulated part of the grid to regain the full picture, which we definitely should do once we generate images of the resulting configuration. But for now, it is reassuring to see that our algorithm does exactly what we want it to do!


Now that our simulation works, we can shift our focus to creating engaging visuals. While the ASCII representation is functional, it is nothing to write home about. In the last part we will write our own bitmap images, and run a simulation with tens of millions of grains to produce a high-resolution fractal image.

Until next time!

Double Pendulum, Part 2

Recently, we talked about different ways how to formulate a classic problem — the double pendulum. We finally arrived at the Lagrangian method. Today, we will write down the Lagrangian of the system and derive the Euler-Lagrange equations of motion. We will also take a look at the Hamiltonian method, yet another approach to solve physics problems. The resulting equations of motion can not be solved analytically, they are too complicated for that. They can, however, be solved numerically, which we will do in the last part on this topic.

Lagrangian Formalism

To recap, we defined a Lagrangian as the difference between the kinetic energy and the potential energy of the system:

L = T - V

The trajectory of the system must then satisfy the Euler-Lagrange equation of motion (EoM):

\dfrac{\textrm{d}}{\textrm{d}t}\dfrac{\partial L}{\partial \dot{x}} = \dfrac{\partial L}{\partial x}

First, let us write down the horizontal x and vertical y cartesian coordinates in terms of the relevant degrees of freedom. For the first pendulum we have:

x_1 = l_1 \sin(\theta_1)\qquad y_1 = -l_1 \cos(\theta_1)

This should be self-evident. The position of the second pendulum will look similar, it will just be offset by the position of the first pendulum:

x_2 = l_1\sin(\theta_1) +  l_2 \sin(\theta_2) \qquad y_2 = -l_1\cos(\theta_1) - l_2 \cos(\theta_2)

To express the Lagrangian, we need the potential energy, which depends on the vertical coordinates of the two pendulums, and the kinetic energy that is proportional to the square of the velocity. So, in order to get the velocity, we derive the coordinates above in time:

v_{x\,1} =  l_1 \cos(\theta_1)\dot{\theta}_1 \qquad v_{y\,1} = l_1 \sin(\theta_1)\dot{\theta}_1

v_{x\,2} =  l_1 \cos(\theta_1)\dot{\theta}_1 + l_2 \cos(\theta_2)\dot{\theta}_2\qquad v_{y\,2} = l_1 \sin(\theta_1)\dot{\theta}_1 + l_2 \sin(\theta_2)\dot{\theta}_2

The square of the velocity of the first pendulum is then:

v^2_1 = l_1^2 \sin^2(\theta_1)\dot{\theta}_1^2 + l_1^2 \cos^2(\theta_1)\dot{\theta}_1^2 = l_1^2 \dot{\theta}_1^2 \left[\sin^2(\theta_1) + \cos^2(\theta_1)\right] = \left(l_1 \dot{\theta}_1\right)^2

The velocity of the second one is a bit trickier:

v^2_2 = l_1^2 \cos^2(\theta_1)\dot{\theta}_1^2 +2l_1l_2\cos(\theta_1)\cos(\theta_2)\dot{\theta}_1\dot{\theta}_2+l_2^2\cos^2(\theta_2)\dot{\theta}_2^2 +

\, \qquad +\  l_1^2 \sin^2(\theta_1)\dot{\theta}_1^2 +2l_1l_2\sin(\theta_1)\sin(\theta_2)\dot{\theta}_1\dot{\theta}_2+l_2^2\sin^2(\theta_2)\dot{\theta}_2^2

Here we can again use the sum rule for the squares of trigonometric functions. This leaves us with:

v^2_2 = \left(l_1\dot{\theta}_1\right)^2 + \left(l_2\dot{\theta}_2\right)^2 +2l_1l_2\dot{\theta}_1\dot{\theta}_2\Big[\sin(\theta_1)\sin(\theta_2) + \cos(\theta_1)\cos(\theta_2)\Big]

The expression inside the square bracket is awfully reminiscent of a trigonometric identity — a quick Wikipedia visit later we can write:

v^2_2 = \left(l_1\dot{\theta}_1\right)^2 + \left(l_2\dot{\theta}_2\right)^2 +2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2)


Now that we have rewritten our coordinates and velocities in terms of the two angular degrees of freedom of the system, we should be able to express the kinetic and potential energy of the system easily:

T = \dfrac{1}{2}m_1\left(l_1 \dot{\theta}_1\right)^2 + \dfrac{1}{2}m_2\left[\left(l_1\dot{\theta}_1\right)^2  + \left(l_2 \dot{\theta}_2\right)^2  + 2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2)\right]

V = -gl_1\left(m_1 + m_2 \right)\cos(\theta_1)  - m_2gl_2\cos(\theta_2)

The Lagrangian is then just the difference of these two. We can now take this Lagrangian and plug it into the Euler-Lagrange equation of motion(s). We expect two equations, one for each angular coordinate. Notice that while the kinetic energy only depends on both the velocities and position, the potential energy is solely a function of the coordinates themselves. This simplifies the next steps a tiny bit.

For \theta_1 we can write:

\dfrac{\textrm{d}}{\textrm{d}t}\dfrac{\partial L}{\partial \dot{\theta}_1} =  \dfrac{\textrm{d}}{\textrm{d}t} \left[m_1l_1^2\dot{\theta}_1 +  m_2l_1^2\dot{\theta}_1 + m_2l_1l_2\dot{\theta}_2\cos(\theta_1 - \theta_2)\right]

\qquad = \left(m_1 + m_2\right)l_1^2\ddot{\theta}_1 + m_2l_1l_2\ddot{\theta}_2\cos(\theta_1 - \theta_2) - m_2l_1l_2\dot{\theta}_2\sin(\theta_1 - \theta_2)\,(\dot{\theta}_1 - \dot{\theta}_2 )

The right hand side of the Euler-Lagrange EoM looks more tame in comparison:

\dfrac{\partial L}{\partial \theta_1} = -gl_1\left(m_1 + m_2\right)\sin(\theta_1) - m_2l_1l_2 \dot{\theta}_1\dot{\theta}_2\sin(\theta_1 - \theta_2)

We can do the same for the other coordinate, \theta_2:

\dfrac{\textrm{d}}{\textrm{d}t}\dfrac{\partial L}{\partial \dot{\theta}_2} = \dfrac{\textrm{d}}{\textrm{d}t} \left[m_2l_2^2\dot{\theta}_2  + m_2l_1l_2\dot{\theta}_1\cos(\theta_1 - \theta_2)\right]

\qquad =  m_2l_2^2\ddot{\theta}_2 + m_2l_1l_2\ddot{\theta}_1\cos(\theta_1 - \theta_2) - m_2l_1l_2\dot{\theta}_1\sin(\theta_1-\theta_2)\,(\dot{\theta}_1 - \dot{\theta}_2)

and

\dfrac{\partial L}{\partial \theta_2} = -gm_2l_2\sin(\theta_2) + m_2l_1l_2 \dot{\theta}_1\dot{\theta}_2\sin(\theta_1 - \theta_2)


Finally, we can combine all the above and write down the actual Euler-Lagrange EoM that determine the behavior of the system [drumroll intensifies]:

(1)\qquad \left(m_1 + m_2\right)l_1\ddot{\theta}_1 + m_2l_2\ddot{\theta}_2\cos(\theta_1 - \theta_2) - m_2l_2\dot{\theta}_2\sin(\theta_1 - \theta_2)\,(\dot{\theta}_1 - \dot{\theta}_2 ) =

-g\left(m_1 + m_2\right)\sin(\theta_1) - m_2l_2 \dot{\theta}_1\dot{\theta}_2\sin(\theta_1 - \theta_2)

(2) \qquad m_2l_2\ddot{\theta}_2 + m_2l_1\ddot{\theta}_1\cos(\theta_1 - \theta_2) - m_2l_1\dot{\theta}^2_1\sin(\theta_1-\theta_2) =-gm_2\sin(\theta_2)

These equations can now be solved to determine the trajectory of the point masses. However, the mathematical nature of these coupled, second order differential equations makes any attempt at a solution very painful. In general, first order differential equations are much nicer to work with.

Hamiltonian Formalism

There is another way to describe physical problems. One that often makes use of the Lagrangian developed above, but which results in a set of two first order differential equations for each degree of freedom, instead of a single second order equation. The Lagrangian is defined as the difference between the kinetic and potential energy of a system. But the sum of these might also be of interest — after all, it is the total energy. And the total energy of an isolated system should be conserved, something that might prove helpful in our analysis.

Let us perform a Legendre transformation on the Lagrangian, moving from coordinates q_i and the corresponding velocities \dot{q}_i, to coordinates and the corresponding generalized momenta p_i. The Legendre transformation of the Lagrangian is the Hamiltonian:

H\left(p_i, q_i\right) \equiv \sum_i p_i \dot{q}_i - L

and in the absence of an explicit time dependence of the Lagrangian, it represents the total energy of the system. The generalized momenta are nothing other than:

p_i \equiv \dfrac{\partial L}{\partial \dot{q}_i}

This is something we already calculated for our coordinates \theta_1 and \theta_2. Once we write down the Hamiltonian (as a function of coordinates and their momenta only), we can compute the trajectory of the system using Hamilton’s equations of motion:

\dot{q}_i = \dfrac{\partial H}{\partial p_i} \qquad \textrm{and} \qquad \dot{p}_i = \dfrac{\partial H}{\partial q_i}


We first express the generalized momenta:

(3) \qquad \dfrac{\partial L}{\partial \dot{\theta}_1}   = p_1 = \left(m_1 + m_2\right)l_1^2\dot{\theta}_1 + m_2l_1l_2\dot{\theta}_2\cos(\theta_1 - \theta_2)

(4) \qquad \dfrac{\partial L}{\partial \dot{\theta}_2} = p_2 = m_2l_2^2\dot{\theta}_2 + m_2l_1l_2\dot{\theta}_1\cos(\theta_1 - \theta_2)

which means that:

\sum_i p_i \dot{q}_i = \left(m_1 + m_2\right)l_1^2\dot{\theta}_1^2 +  m_2l_2^2\dot{\theta}_2^2 + 2 m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2)

With the Lagrangian being:

L = \dfrac{1}{2}\left(m_1 + m_2\right) l_1^2 \dot{\theta}_1^2 + \dfrac{1}{2}m_2l_2^2\dot{\theta}_2^2    + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2) +

+ gl_1\left(m_1 + m_2 \right)\cos(\theta_1) + m_2gl_2\cos(\theta_2)

the Hamiltonian must be:

H \equiv \sum_i p_i \dot{q}_i - L = \dfrac{1}{2}\left(m_1 + m_2\right) l_1^2 \dot{\theta}_1^2 + \dfrac{1}{2}m_2l_2^2\dot{\theta}_2^2 + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2) +

- gl_1\left(m_1 + m_2 \right)\cos(\theta_1) - m_2gl_2\cos(\theta_2)

Unfortunately, we are not done yet. The Hamiltonian is a function of the coordinates \theta_1, \theta_2, p_1, and p_2. We must express \dot{\theta}_1 and \dot{\theta}_2 as functions of generalized positions and momenta. To that end, we isolate \dot{\theta}_2 in eq. (4):

p_2 = m_2l_2^2\dot{\theta}_2 + m_2l_1l_2\dot{\theta}_1\cos(\theta_1 - \theta_2)

m_2l_2^2\dot{\theta}_2 = p_2 - m_2l_1l_2\dot{\theta}_1\cos(\theta_1 - \theta_2)

\dot{\theta}_2 = \dfrac{p_2}{m_2l_2^2} - \dfrac{l_1}{l_2}\dot{\theta}_1\cos(\theta_1 - \theta_2)

and plug into eq. (3):

p_1 = \left(m_1 + m_2\right)l_1^2\dot{\theta}_1 + m_2l_1l_2\cos(\theta_1 - \theta_2)\left[\dfrac{p_2}{m_2l_2^2} - \dfrac{l_1}{l_2}\dot{\theta}_1\cos(\theta_1 - \theta_2)\right]

Let us gather all terms that contain \dot{\theta}_1:

\dot{\theta}_1\left[\left(m_1 + m_2\right)l_1^2 -m_2l_1^2\cos^2(\theta_1 - \theta_2)\right] = p_1 - p_2\dfrac{l_1}{l_2}\cos(\theta_1 - \theta_2)

whence:

\dot{\theta}_1 = \dfrac{p_1 - p_2\dfrac{l_1}{l_2}\cos(\theta_1 - \theta_2)}{\left(m_1 + m_2\right)l_1^2 -m_2l_1^2\cos^2(\theta_1 - \theta_2)}

Multiplying both the numerator and the denominator by l_2, and rearranging some terms, we get:

\dot{\theta}_1 = \dfrac{l_2p_1 - l_1p_2\cos(\theta_1 - \theta_2)}{l_1^2l_2\left[m_1 + m_2 -m_2\cos^2(\theta_1 - \theta_2)\right]}

where we use the good ol’ trigonometric identity \sin^2(x) + \cos^2(x) = 1 to find the final expression:

\dot{\theta}_1 = \dfrac{l_2p_1 - l_1p_2\cos(\theta_1 - \theta_2)}{l_1^2l_2\left[m_1 +m_2\sin^2(\theta_1 - \theta_2)\right]}

We could perform a similar procedure, expressing \dot{\theta}_1 from eq. (3), plugging into eq. (4) and solving for \dot{\theta}_2 to find:

\dot{\theta}_2 = \dfrac{l_1\left(m_1 + m_2\right)p_2 - l_2m_2p_1\cos(\theta_1-\theta_2)}{l_1l_2^2m_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

These formulas represent two of the four Hamilton’s equations we want to derive. What we need to do now is plug the expressions for \dot{\theta}_1 and \dot{\theta}_2 into the Hamiltonian, and then derive with respect to the generalized positions \theta_1 and \theta_2 in order to find the two other Hamilton’s equations. This is a rather lengthy calculation, one I would like to avoid if at all possible. If time permits, I might later update this post with the full derivation. I hope you will forgive me my laziness.

Acknowledging this caveat, we finally arrive at Hamilton’s equations of motion:

\dot{\theta}_1 = \dfrac{l_2p_1 - l_1p_2\cos(\theta_1-\theta_2)}{l_1^2l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

\dot{\theta}_2 = \dfrac{l_1\left(m_1 + m_2\right)p_2 - l_2m_2p_1\cos(\theta_1-\theta_2)}{l_1l_2^2m_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

\dot{p}_1 = -\left(m_1 + m_2\right)gl_1\sin(\theta_1) - A + B

\dot{p}_2 = -m_2gl_2\sin(\theta_2) + A - B

where

A \equiv \dfrac{p_1p_2\sin(\theta_1 - \theta_2)}{l_1l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}

B \equiv \dfrac{l_2^2m_2p_1^2 + l_1^2\left(m_1 + m_2\right)p_2^2 - l_1l_2m_2p_1p_2\cos(\theta_1 - \theta_2)}{2l_1^2l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]^2}\sin\left[2(\theta_1 - \theta_2)\right]


There, done. We eventually succeeded in writing down the equations of motion of the two masses in a double pendulum. And we did this in a form that is very conducive to numerical simulations. In the next part, we will finally get to write a Python script to simulate the time evolution of the system.

Until next time!

Abelian Sandpile: a C++ Project, Part 2

In the previous instalment, we have introduced the Abelian sandpile model, and discussed our approach to simulate it. Today, we are going to implement the container class that will represent our square grid. We will also write the topple() method, which is going to be in charge of toppling the towers of sand in each square tile. A clever trick is going to save us millions of operations. Let’s get started!


To recap, this is the skeleton of the Sandpile class that we’ve drawn up last time:

File “Sandpile.h”

#pragma once

class Sandpile
{
private:
    Grid <int> m_grid;
    int m_height;

public:
    Sandpile(int height);
    ~Sandpile();
    topple();
    topple_all();
    draw();
};

First things first — what does #pragma once do? Albeit not standard, it is a widely supported C++11 replacement of the well-known construct

#ifndef CLASS_H
#define CLASS_H
...
#endif

It makes sure that a given piece of code is only included once in the final program. But if you absolutely prefer the standard way, feel free to use it instead!

The Grid Class

With this being said, let us focus our attention on the Grid<int> m_grid line. We’ve established that m_grid is going to represent our 2D square grid. Whatever the internal implementation of
Grid is going to be, it will need to dynamically allocate an array to store the amounts of sand grains on each tile. There are several ways how to dynamically allocate a 2D array.

We could instantiate a pointer pointing to an array of pointers pointing to an array of integers. This would certainly work, but is a bit awkward. For one thing, we will need to traverse the array over and over. But this method doesn’t guarantee us that the individual integer arrays (representing, for example, the rows of our grid) will be allocated in a contiguous block of memory. Moreover, and more importantly, using raw pointers like this forces us to pay very close attention to all the ways a memory leak could happen.

Alternatively, we could use some container in the standard library.
std::array needs to be initialized with a constant value. This would make it necessary to hard-code the size of the square grid, which might not be what we want to do. What if we want to choose the size of the image at runtime? Or let the size of the simulated square grid adapt to the number of grains? std::vector comes to mind, but we don’t really need to resize the grid once declared.

C++ allows us to use smart pointers to facilitate smart memory management. Can we make a 2D array work using smart pointers, while avoiding the creation of arrays containing more pointers? Yes, we can! The idea is to really declare a simple 1D array of integers that will store all our data in a single, contiguous block of memory. We then write a quick access function which makes it behave like a 2D array.

Since we may want to reuse this container in the future, we should write a templated version of our class. This just means that instead of specifying the type of elements stored in our 2D array, we use an anonymous, placeholder type (usually dubbed T). Once the compiler sees we are trying to instantiate a Grid of a certain type, it copies the templated implementation and automatically creates a class of the required type. Due to how the compiling of template functions works, it is generally the easiest to implement everything in the header file:

File “Grid.h”

#pragma once
#include <memory>

template <class T>
class Grid
{
private:
    std::unique_ptr<T[]> m_grid{ nullptr };
    int m_size_x{};
    int m_size_y{};

public:
    ...
};

The Grid class shall have three member variables. Two integers, m_size_x and m_size_y, which specify the size of the grid. We are going to use a square grid, but the option of a rectangular grid might make this class more reusable in the future. The final member is a pointer to an array of type T (note the use of square brackets to signify an array of Ts instead of just a single T). We are using a smart pointer living in the <memory> header. A unique pointer ensures that the allocated memory will duly be deallocated — hence no need for an explicit destructor. We initialize it to nulltpr.

Next, we need to write the constructor. It should take two integer parameters, initialize the vertical and horizontal size of the grid, and allocate memory for the grid itself. We initialize all values to 0. Finally, we will need an accessor to the grid elements. The trick here is to fold the single allocated 1D array in such a way that it behaves like a 2D grid. Since the at(int x, int y) is the only public member (the only interface with the user), we achieve our goal by using two integers to specify “rows” and “columns”, then mapping them onto a single array index. We want this function to return a reference to the array element in order to be able to change the latter. The resulting Grid class should look like this:

File “Grid.h”

#pragma once
#include <memory>

template <class T>
class Grid
{
private:
    std::unique_ptr<T[]> m_grid{ nullptr };
    int m_size_x{};
    int m_size_y{};

public:
    Grid(int size_x, int size_y)
        : m_size_x(size_x), 
        m_size_y(size_y), 
        m_grid(new T[size_x*size_y]{ 0 })
    {
    }
    
    T& at(int x, int y)
    {
        return m_grid[x * m_size_y + y];
    }
};

Implementing the Topple Method

We will now turn our attention to the Sandpile class. As we discussed in Part 1, the toppling will be implemented via two methods. One function will iterate over all grid elements, evaluate the number of sand grains, topple, and repeat this until all squares contain less than four grains. The other one will actually implement the toppling.

The Abelian Sandpile is defined by the following toppling rule: if any tile holds 4 or more sand grains, this pile will topple, reducing the number of grains on this tile by 4 and sending 1 grain to each of the four adjacent tiles.

We can improve this by performing all the toppling at once. Assuming a number N of grains of sand before the toppling, we will keep N\!\mod 4 after, and each immediate neighbour will receive N/4.

We also said we would exploit the square symmetry of the whole setup. How can we do it in practice? Consider numbering the tiles that represent one quadrant of the “infinite” grid starting from 0, placing the center of the pile in the tile with indices (0, 0). The boundaries with the mirror images should thus extend along lines defined by either of the two indices equal to zero. These boundary tiles should not topple in the direction of the mirror image. Likewise, whenever the row or column next to the boundary topples, due to symmetry, the mirrored row or column on the other side should also topple. The boundary tiles should thus receive twice the amound of sand grains whenever a tile in the neighbouring row or column topples:

The naive implementation of the topple() method would look something like this:

File “Sandpile.cpp”

void Sandpile::topple(int x, int y)
{
    // mirror image boundary for x = 0
    if ((y > 1) && (y < m_size) && (x == 0))
    {
    	m_grid.at(x, y) -= 4;
    	m_grid.at(x, y + 1) += 1;
    	m_grid.at(x, y - 1) += 1;
    	m_grid.at(x + 1, y) += 1;
        return;
    }
    
    // similarly for all y = 0
    ...
    
    // the normal case
    m_grid.at(x, y) -= 4;
    m_grid.at(x + 1, y) += 1;
    m_grid.at(x - 1, y) += 1;
    m_grid.at(x, y + 1) += 1;
    m_grid.at(x, y - 1) += 1;
    
    // boundary condition for x = 1
    if (x == 1)
        m_grid.at(x - 1, y) += 1;
    
    // boundary condition for y = 1
    if (y == 1)
        m_grid.at(x, y - 1) += 1;

    return;
}

This should do what we want it to do, but not very efficiently. Notice the code reuse all over the place. Not to mention, for virtually most cases, two if-conditions will have to be evaluated before the toppling begins — even though the boundary is insignificant relative to the size of the whole grid, it slows down the simulation for every single tile. We can do better.

First of all, we only toppled four grains. Let us take as an argument the height h of the sandpile in the given tile. We will reduce the height of the pile to h \mod 4. Each neighbouring tile receives h/4 sand grains (note we require integer division).

Secondly, and more importantly, we can use a neat trick to prevent multiple boundary checks at each iteration. We simply increase the size of our grid by one in both dimensions. This creates a buffer which we can indiscriminately fill without worrying about accessing non-existent grid elements. This buffer will not be toppled, and hence we don’t incur any performance penalty. Our numbering scheme is shifted by one, the boundary with the mirror images now extends along the lines defined by either index being equal to one. In effect, we trade a tiny amount of additional memory for a noticeably faster algorithm!

With these changes, the topple() method should be:

File “Sandpile.cpp”

void Sandpile::topple(int x, int y, int height)
{
    m_grid.at(x, y) = height % 4;
    m_grid.at(x, y + 1) += height / 4;
    m_grid.at(x, y - 1) += height / 4;
    m_grid.at(x + 1, y) += height / 4;
    m_grid.at(x - 1, y) += height / 4;
    
    // boundary elements receive twice the grains
    // numbering shifted by one due to buffer
    if (x == 2)
        m_grid.at(x - 1, y) += height / 4;
 
    if (y == 2)
        m_grid.at(x, y - 1) += height / 4;
    
    return;
}

There, that is much more concise — and faster as well.


We have implemented a class to represent our square grid, and written an algorithm for toppling which makes use of the square symmetry of the system. In the next instalment, we will implement the whole toppling process responsible for iterating over all grid tiles repeatedly until all piles are toppled. We will also print out the resulting configuration for testing purposes.

Until next time!

Abelian Sandpile: a C++ Project, Part 1

Let’s imagine an infinite square grid. In one of the squares we place a certain amount of sand grains. We simulate this “sandpile” as follows: if there are more than four grains on any given tile, the tower of sand topples. Its height is reduced by four, while the tiles to the North, West, South, and East of the toppled square receive each one grain. Repeat this toppling until no square field contains more than three grains. Let’s color the final tiles depending on how many grains of sand they hold (0, 1, 2, or 3). What do we get?

Pretty random. Let’s zoom out!
There seems to be some pattern emerging. Let’s keep going!
Whoa…
You gotta be kiddin’ me!

As with many other fractal structures, this beauty is also the result of a deceivingly simple set of rules. However, the nature of the problem makes the simulation of the sandpile extremely computationally demanding. If we want to produce HD images like those above, we will need all the speed a language like C++ can muster.


In this series, we will write our own C++ program which will eventually draw beautiful images of sandpiles containing millions of grains. Before we get there, though, we have a lot of ground to cover. Note that this is my first real C++ project. I will try to explain as much as I can along the way. This doesn’t mean, however, that I will always present the optimal solution. If at any point you feel like you know better, don’t hesitate to do it your own way (and if possible, let me know in the comments what I could improve)!

I recommend those who are new to C++ (or new to the functionalities added in C++11 and higher) to have their favourite C++ textbook at hand in case some point or other is not entirely clear. Personally, I found the free book at https://www.learncpp.com/ a very good modern C++ primer. The Udemy course https://www.udemy.com/learn-advanced-c-programming/ was exceedingly helpful as well. If you are interested in learning more about the language, I warmly recommend the above resources.

First Step

Before we delve into coding, let us first analyse the problem at hand, identify our requirements, choose an appropriate algorithm, and think of ways how to reduce the computational effort needed.

We will need an array of integer values representing the square grid. While in the end all values will range from 0 to 3, during the toppling process some squares may contain millions of grains — necessitating an int. Since we want to have a nice, high-definition image, we need the size of the array to be on the order of 1000 by 1000 elements. With millions of integer values, this array may potentially be too big for the stack. Hence, we will allocate it dynamically on the heap.

We instantiate our array with some initial amount of sand grains in one of the squares. What do we do next? Instead of toppling the uppermost four grains again and gain, we can topple the whole pile at once, sending one fourth of the toppled amount to each neighbour, and leaving us with the original height modulo four on the central tile. Now we need to iterate over all tiles repeatedly, toppling the ever smaller towers of sand. Note that even though we topple all the grains in a given square at once, the same tile might need toppling again if enough sand grains fall on it from its four neighbours.


Can we say something about the complexity of the task? Let us designate by N the number of initial sand grains. In the final picture, each square contains zero, one, two, or three grains. With a constant average of 1.5, this means that the number of tiles that may topple is linearly proportional to the initial N. If we double N, we double the number of tiles that need to be toppled. However, each toppling may cause a cascade, an avalanche of topples that could spread repeatedly through the whole grid. At best, the number of topples is proportional to the number of tiles, and with a direct simulation approach, we get a complexity of, at best, \mathcal{O}\left(N^2\right). There are probably better algorithms out there which make use of repeating sequences, clever tricks, or some deep mathematical observations — but I am trying to keep this project as simple as I can, so we will ignore “the smarter approach” and use our simple, brute force algorithm.

Nevertheless, we will make use of symmetry. The square grid and the symmetric toppling rule ensure that the final pattern will have mirror symmetry in both the horizontal and vertical direction. Hence, we only really need to simulate one quadrant. Once we have the final configuration, we can simply mirror it to display the final image. A trivial trick that will potentially save us billions of operations.

The Sandpile Class

Before we conclude this first part, let us at least set up the scaffolding of our program. Having a clear idea of what we want/need to do upfront makes it so much easier to implement it!

The simulation will consist of a dynamically allocated array on which we will operate with our toppling algorithm. Since this seems pretty self-contained, an object-oriented approach feels like a good idea. Let us create a class named Sandpile. As a member, it will hold an array m_grid representing the square grid. It should also have the initial height m_heightas a member variable. We will implement the toppling process in a topple_all() method. This method will iterate repeatedly over all elements. If they contain more than four sand grains, it will call the topple() method. This way we make the algorithm a bit easier to read and understand. Finally, we will draw the resulting distribution of sand grains via the draw() method:

File “Sandpile.h”

#pragma once

class Sandpile
{
private:
    Grid <int> m_grid;
    int m_height;

public:
    Sandpile(int height);
    ~Sandpile();

    topple();
    topple_all();
    draw();
};

In the next part, we will implement the templated container class for our square grid, which I named Grid. We will also implement the topple() method, keeping in mind our exploit of the square symmetry.

Till next time!