rad235_solution

#!/usr/bin/env python2
#
# Subitop modelling course, Edinburgh 2017
# Jeroen van Hunen, March 2017
# rad235U_solution.py
# purpose: calculates radioactive decay of 235U to 207Pb:

import numpy as np
import pylab as plt

nt    = 10                       # nr of timesteps
time  = np.linspace(0,4.567,nt)  # Define time array over 4.567 Gyrs
dt    = time[1]-time[0]          # Calculate size of each timestep
    
N_FE   = np.zeros(nt);  # 235U / 204Pb 
N_BE   = np.zeros(nt);  # 235U / 204Pb 
N_CN   = np.zeros(nt);  # 235U / 204Pb 

N_ini  = 5.5            
N_FE[0]= N_ini
N_BE[0]= N_ini
N_CN[0]= N_ini         
lam = 0.9848         # halflife = 704 Myrs

for it in range(1,nt):   
    N_FE[it] = N_FE[it-1] * (1-dt*lam)            # Forward Euler
    N_BE[it] = N_BE[it-1] / (1+dt*lam)            # Backward Euler
    N_CN[it] = N_CN[it-1] * (2-dt*lam)/(2+dt*lam) # Crank-Nicholson

N_ana  = N_ini * np.exp(-lam*time)    # analytical solution

plt.figure(1)
plt.clf()
plt.plot(time, N_FE, 'o-', label='FE')
plt.plot(time, N_BE, 'o-', label='BE')
plt.plot(time, N_CN, 'o-',label='CN')
plt.plot(time, N_ana, label='analytical')
plt.legend()
plt.xlabel('t(Gyrs)')
plt.ylabel('235U / 204Pb')
plt.show()