Radioactivity (Work in progress)

Imagine you have a bunch (an N₀-numbered bunch) of radioactive nuclei of an element X. These nuclei, sadly, will start to decay and other, more stable nuclei will be formed, along with some ionising radiation. Our aim of the game is to work out how many X nuclei remain after a certain time period. Radioactive decay is a random process, but the probability of decay (P) per unit time of a specific nuclei decaying is fixed (as all nuclei are of the same element X) - this is λ, the decay constant.

We can use this to derive our model for decay of this bunch over time. Let's start with a differential equation. The rate of change of undecayed nuclei (N) is proportional to N (as N increases, the number of decays per period will increase, leading to a greater negative change in N). We also require a constant of proportionality, which we'll call k. This gives equation B. Manipulating this to get the fractional change of N ($\frac{\mathrm{dN}}{N}$) in an instant (a period of duration dt) yields equation C. Probability of decay is, on average, equal to the (number of decays)/(total number of nuclei) over each dt period, which is -$\frac{\mathrm{dN}}{N}$. As N approaches infinity, -$\frac{\mathrm{dN}}{N}$ will approach P. So, we can approximate, and say that -$\frac{\mathrm{dN}}{N}$ = P, and so we can equate A and C to yield D. D shows that our coefficient k is equal to our decay constant λ.

Now let's solve D to get a function for N in terms of t. From this, we will be able to model how many nuclei still remain after a certain time period. There's a delta term on both sides, and so we can integrate on both sides to remove them. Using E, and plugging in N = N₀ at t = 0, yields that N₀ = e^c. We can now reach our final form for our model in F. F is the green plot on the graph, and you can play with the sliders to see how the theoretical plot changes when N₀ and λ change.

Let's move onto the main event: the radioactivity simulation. I've simulated the random nature of the decay of radioactive nuclei by this algorithm:

• Start with N nuclei
• For each timestep (dt increment), run a random number generator (between 0 and 1) N times (simulating each nucleus' decay probability)
• For each RNG run, if the value returned is greater than P (=λ*dt), add 1 to the new nuclei count for the next timestep
• Repeat for as many timesteps as you want (I've set dt to 0.1s so the graph simulation appears as smooth as possible without being laggy)

This simulation is the shown by the red plot on the graph. It can be seen that the red simulation plot is not always particularly close to the green theory plot. This is due to the random nature of the sim. In order to prove that radioactive decay can be modelled by an exponential, I have implemented another plot which takes an average of multiple runs - this is the blue plot. To compare between the random plot and the averaged plot, a chi-squared (Χ²) test has been deployed, displaying the "goodness of fits" of both simulations with the theory. We can indeed now prove that the exponential model is very accurate, as increasing the number of samples for the average plot decreases the Χ² value significantly, indicating a better fit between theory and average.

$\begin{array}{cccc}P& =& \lambda \mathrm{dt}& \left(A\right)\\ \frac{\mathrm{dN}}{\mathrm{dt}}& =& -kN& \left(B\right)\\ \frac{\mathrm{dN}}{N}& =& -k\mathrm{dt}& \left(C\right)\\ \frac{\mathrm{dN}}{N}& =& -\lambda \mathrm{dt}& \left(D\right)\\ \int \frac{1}{N}\mathrm{dN}& =& \int -\lambda \mathrm{dt}\\ lnN& =& -\lambda t+c\\ N& =& e^{-\lambda t+c}\\ N& =& e^c·e^{-\lambda t}& \left(E\right)\\ N& =& N_0{e}^{-\lambda t}& \left(F\right)\end{array}$