Monte Carlo Simulations
Many processes in physics are stochastic meaning they are determined randomly. Computational physics methods that involve stochastic processes are called Monte Carlo methods. All Monte Carlo methods rely on being able to generate random numbers. The numpy.random module has methods that generate random numbers from various probability distributions. One of the simplest to understand is rand(N). It creates an array of length N and populates it with a random sample from a uniform distribution over the interval [0, 1). Leaving out the argument N generates a single random number in the interval. The function randn(N) generates a random set of N values drawn from a univariate "normal" (Gaussian) distribution of mean 0 and variance 1. We'll see how we can use both of these functions in the examples below.
Uniform Distribution of Random Numbers
Uniformly distributed random numbers are random numbers in a specified range, for which any one number in the range is just as likely as any other. The numpy.random module has a function called rand() that generates random numbers drawn from a uniform distribution in the interval [0,1). You can get uniformly distributed random numbers in any other range by multiplication or division. For example 2*rand() would generate a number in the interval [0,2) and rand() - 0.5 would generate numbers in the interval [-0.5,0.5).
Simulating radioactive decay
Suppose we have a collection of \(N\) radioactive particles with a half life \(t_{1/2}\). We know that we will have about \(\frac{1}{2}N\) particles after a time equal to \(t_{1/2}\), about \(\frac{1}{4}N\), after \(2t_{1/2}\), and so on. However, we can't predict exactly how many particles we will have at any given time because process is stochastic. We can, however, simulate the decay1, but before we explore the simulation lets review the radioactive decay process.
For radioactive decay, the probability of a single particle decaying per unit time is constant and is called the decay constant. Suppose we have \(N\) radioactive particles at time \(t\) and that the decay constant is \(\lambda\). The probability of any one particle decaying in a time interval \(\Delta t\) is \(P = \lambda \Delta t\). The number of particles that would decay in the time interval is \(P N = \lambda \Delta t N\) and hence change in the total number of radioactive particles is $$ \Delta N = - \lambda \Delta t N. $$ The minus sign indicates that the number of radioactive particles is decreasing as they decay. Rearranging this equation gives the decay rate $$ \frac{\Delta N}{\Delta t} = - \lambda N $$ and in the limit as \(\Delta t \rightarrow 0\) $$ \frac{dN}{dt} = - \lambda N. $$ Solving this differential equation gives $$ N(t) = N_0 e^{-\lambda t}, $$ where \(N_0\) is the number of particles at time \(t = 0\). The lifetime of the particle is defined as \(\tau = 1/ \lambda\), and the half-life \(t_{1/2} = \ln(2)/\lambda\).
The program radioactive.py simulates radioactive decay, plots the result, and compares it to the analytic calculation above. It's easiest to understand how the simulation works by looking at the program. The code from radioactive.py that is used to do the simulation is shown below. P_decay is the probability that a single particle will decay in a given time interval, P_decay equals \(P = \lambda \Delta t\). This code snippet starts by setting the number of particles left equal to the initial number of particles and creating an array to hold the results. It then proceeds to decide if each of the particles has decayed in the first time interval. The code determines if a given particle has decayed by getting a random number between 0 and 1 by using rand(). It then compares this random number to P_decay. If the random number is less than P_decay then the particle has decayed and the program reduces the number of particles left by one. It then repeats this process for the next time interval until all the particles are gone.
# Set the initial value of the number of particles left
Nleft = N0
# Creat an array to hold the number of particles left at the end of each time interval
N_t=np.array([Nleft])
# Repeat until all the particles have decayed
while Nleft != 0:
# Set the number of particles to be the number at the end of the last time
# interval.
Npart = Nleft
# Go through all the particles to see if they decay. If a particle decays,
# reduce the number left by one.
for N in range(Npart): # Go through all the particles
r = rand()
if (r < P_decay):
Nleft = Nleft-1
# Append the number of particles left after this time interval to the array
# containing the number of partles left from the previous time intervals.
N_t = np.append(N_t,Nleft)
Exercise 1
Download radioactive.py from the downloads page. This program does a simulation of radioactive decay with \(\lambda = 0.01\) s\(^{-1}\) and \(\Delta t = 0.01\) s. Run it for several values of \(N_0\) from 10 to 1000. For which values of \(N_0\) does the exponential law give results most consistent with the simulation.
Exercise 2
Write your own program to do ten different decay simulations with the same \(\lambda\) and \(\Delta t\) values used in radioactive.py. Plot all of these simulations along with the exponential decay law, and the average of the ten simulations. Use N0 = 15 for each simulation. How do the individual simulations compare with the exponential decay law? How does the average of the ten simulations compare with the exponential decay law? Is this what you expected? Explain. (Feel free to steal liberally from radioactive.py when writing your program.)