Radioactive Decay¶
Faisal Qureshi
faisal.qureshi@ontariotechu.ca
http://www.vclab.ca
Copyright information¶
© Faisal Qureshi
License¶
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import random
In [19]:
l = 0.05 # Lambda, represents the average probability
# per nucleus of decay occurring per unit time
n0 = 100000 # Total number of nuclei
tmax = 60.
steps = 10000
nta = np.empty([steps,1])
times = np.linspace(0, tmax, steps)
# Analytical solution
nta = n0 * np.exp(-l * times)
# Monte Carlo solution
def nuclear_decay_monte_carlo(n0, p, steps):
print ('n0', n0)
print ('prob of decay', p)
print ('steps', steps)
nt = np.empty(steps)
nt[0] = n0
for i in range(1,steps):
decay_or_not = np.random.random(int(nt[i-1]))
num_not_decayed = np.sum(decay_or_not > p)
nt[i] = num_not_decayed
return nt
prob_of_decay_per_time_step = l * tmax / steps
nt = nuclear_decay_monte_carlo(n0, p=prob_of_decay_per_time_step, steps=steps)
n0 100000 prob of decay 0.0003 steps 10000
In [20]:
# Plotting
plt.figure(figsize=(7,7))
plt.plot(times, nta, label='Analytical Solution')
plt.plot(times, nt, 'r--', label='Monte Carlo Simulation')
plt.xlabel('Time')
plt.ylabel('Number of nuclei')
plt.legend(['Analytical Solution', 'Monte Carlo Simulation'])
plt.show()
In [ ]: