Simulating radioactive decay

The decay of an ensemble of radioactive nuclei over a period of time can be simulated as follows. Consider the time period to be divided into short, discrete intervals of duration $\Delta t \ll \tau$, where $\tau$ is the lifetime for the decay (which is related to the half-life, $t_{1/2}$, through $\tau = t_{1/2}/\ln 2$). The probability that a given nucleus will decay in time $\Delta t$ is $p = \Delta t / \tau$.

At each time-step, the simulation loops over the undecayed nuclei from the previous time-step and draws a random number from the uniform distribution on $[0,1)$: if this random number is less than $p$, the nucleus is considered to have decayed.

The code below defines a function to carry out this simulation for a set of $N_0 = 500$ $\mathrm{^{14}C}$ nuclei with half-life $t_{1/2} = 5730\;\mathrm{years}$. nsims = 10 such simulations are carried out and saved to a comma-separated file, 14C-sim.csv, with a brief, explanatory header.

import random
import numpy as np

def decay_sim(thalf, N0=500, tgrid=None, nhalflives=4):
    """Simulate the radioactive decay of N0 nuclei.

    thalf is the half-life in some units of time.
    If tgrid is provided, it should be a sequence of evenly-spaced time points
    to run the simulation on.
    If tgrid is None, it is calculated from nhalflives, the number of
    half-lives to run the simulation for.


    # Calculate the lifetime from the half-life.
    tau = thalf / np.log(2)

    if tgrid is None:
        # Create a grid of Nt time points up to tmax.
        Nt, tmax = 100, thalf * nhalflives
        tgrid, dt = np.linspace(0, tmax, Nt, retstep=True)
        # tgrid was provided: deduce Nt and the time step, dt.
        Nt = len(tgrid)
        dt = tgrid[1] - tgrid[0]

    N = np.empty(Nt, dtype=int)
    N[0] = N0
    # The probability that a given nucleus will decay in time dt.
    p = dt / tau
    for i in range(1, Nt):
        # At each time step, start with the undecayed nuclei from the previous.
        N[i] = N[i-1]
        # Consider each nucleus in turn and decide whether it decays or not.
        for j in range(N[i-1]):
            r = random.random()
            if r < p:
                # This nucleus decays.
                N[i] -= 1 
    return tgrid, N

N0 = 500
# Half life of 14C in years.
thalf = 5730

# Use Nt time steps up to tmax years.
Nt, tmax = 100, 20000
tgrid = np.linspace(0, tmax, Nt)

# Repeat the simulation "experiment" nsims times.
nsims = 10
Nsim = np.empty((Nt, nsims))
for i in range(nsims):
    _, Nsim[:, i] = decay_sim(thalf, N0, tgrid)

# Save the time grid, followed by the simulations in columns. We save integer
# values for the data and create a comma-delimited file with a two-line header.
np.savetxt('14C-sim.csv', np.hstack((tgrid[:, None], Nsim)),
    fmt = '%d', delimiter=',',
    header=f'Simulations of the radioactive decay of {N0} 14C nuclei.\n'
           f'Columns are time in years followed by {nsims} decay simulations.'

The contents of the output file, 14C-sim.csv, will resemble:

# Simulations of the radioactive decay of 500 14C nuclei.
# Columns are time in years followed by 10 decays.

This file can be read in to a NumPy array with:

arr = np.loadtxt('14C-sim.csv', delimiter=',')

See also Exercise P6.5.7.