Radioactivity

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.integrate import odeint
from scipy.integrate import solve_ivp

Radioactive Decay by solving ODE

Rate of decay of a radioactive nuclei is $$\frac{dN}{dt} = -\lambda N $$

This is a first order differential equation which can be solved by any standard method. But, we must mention initial condition of the decay, i.e. at time $t=0$, what is number of radioactive nuclei $N(0)$ present in the substance? Consider, $$\text{at}~~~ t=0, \hspace{1cm}N(0)=N_0$$

N0 = 1000
T_half = 4.5*60
lamda = 1/ T_half

Creating time array

t_init, t_final, step_size = 0, 1000, 1
t = np.arange(t_init, t_final, step_size)

Model of radioactivity

def model(N, t):
  return -lamda*N

By odeint method from scipy.integrate

sol = odeint(model, N0, t)
N = sol[:, 0]

plt.plot(t, N, ls='--', color='red', label='Parent nucleus')
plt.plot(t, N0 - N, ls='--', color='blue', label='Daughter nucleus')
plt.legend()
plt.show()

By solve.ivp from scipy.integrate

def model(t, N):   # order of the argument is very important; first independent and then dependent variable
  return -lamda*N

sol = solve_ivp(model, t_span=(t[0], t[-1]), y0=([N0]), t_eval= t)

t = sol.t
N = sol.y.T

plt.plot(t, N, ls='--', color='green', label='Parent nucleus')
plt.plot(t, N0 - N, ls='--', color='orange', label='Daughter nucleus')
plt.legend()
plt.show()

Simulation of Radioactive Decay by Monte Carlo method

P = []                # Parent Nucleus
D = []                # Daughter Nucleus
T_half = 4.5*60       # Half life 
lamda = 1/ T_half     # Decay constant
N = 1000              # Initial number of Radioactive nuclei
Np, Nd = N, 0
time = np.arange(0, 1000, 1)

for t in time:
  P.append(Np)
  D.append(Nd)
  decay = 0
  for n in range(Np):
    p = random.random()
    if p < lamda:
      decay += 1
  Np -= decay
  Nd += decay
  
plt.plot(P, linestyle='-', color='red', label='Parent Nucleus')
plt.plot(D, linestyle='-', color='blue', label='Daughter Nucleus')
plt.legend()
plt.xlabel('t', fontsize=14)
plt.ylabel('N(t)', fontsize=14)
plt.show()