15.9.
Radioactive carbon 11 has a decay rate k of 0.0338 per minute—that is, a particular C$^{11}$ atom has a 3.38% chance of decaying in any one minute. Suppose we start with 100 such atoms. We want to simulate their fate over a period of, say, 100 minutes, ending up with a bar graph showing how many atoms remain undecayed after 1, 2, ..., 100 minutes.
We need to simulate when each of the 100 atoms decays. This can be done, for each atom, by generating a random number $r$ for each of the 100 minutes until either $r > k$
100
80
Undecayed atoms
60
40
20
0
20
40
60
80
100
120
Time (minutes)
FIGURE 15.1
Radioactive decay of carbon 11: simulated and theoretical.
(that atom decays) or the 100 minutes are up. If the atom decays at time $t < 100$, increment the frequency distribution $f(t)$ by 1. $f(t)$ will be the number of atoms decaying at time $t$ minutes.
Now convert the number $f(t)$ decaying each minute to the number $R(t)$ remaining each minute. If there are $n$ atoms to start with, after one minute the number $R(1)$ remaining will be $n - f(1)$ since $f(1)$ is the number decaying during the first minute. The number $R(2)$ remaining after two minutes will be $n - f(1) - f(2)$. In general, the number remaining after $t$ minutes will be (in MATLAB notation)
$R(t) = n - \sum(f(1:t))$
Write a script to compute $R(t)$ and plot its bar graph. Superimpose on the graph the theoretical result, which is
$R(t) = 100 \exp^{-kt}$
Typical results are shown in Figure 15.1.
