session1

Session 1: Radioactive Decay, Part 1

Lecture slides for this session can are here.

Uranium 235 ( ²³⁵U ) decays to Lead 207 ( ²⁰⁷Pb ), with a half-life τ=704τ=704 Myrs. In other words, in 704 Myrs, half of the available ²³⁵U decays to ²⁰⁷Pb . More generally, each unit of time (say each million years) a constant fraction of ²³⁵U decays to ²⁰⁷Pb . The removal of ²³⁵U (to produce ²⁰⁷Pb ) is therefore directly proportional to the available ²³⁵U itself (check this for yourself!). Let’s simplify notation and call ²³⁵U simply NN . Mathematically, this radioactive decay is then written as:

dNdt=−λNdNdt=−λN, with λλ being the ‘decay constant’.

NOTE: An equation like this one, containing a variable and its derivative, is called a differential equation (henceforth referred to as DE) It nicely describes the physics of (in this case) radioactive decay, but it doesn’t provide an immediate solution for the unknown variable NN; we have to solve it first. In some cases, including this one, it is possible to solve the DE analytically. In other cases, a numerical approximation must be used.

Quantities (or rather concentrations) of ²³⁵U and ²⁰⁷Pb are usually measured with respect to the concentration of a stable isotope of the daughter element, in this case ²⁰⁴Pb. At the start of Earth’s life at 4.567 Ga, the concentrations in the mantle of ²³⁵U relative to ²⁰⁴Pb was Nini=²³⁵U ²⁰⁴Pb=5.5.

1.1. Run the following script rad235U.py in your Python editor. Study it briefly, and run the code to see what it does. Make sure you understand what the code is doing.

# Subitop modelling course, Edinburgh 2017
# Jeroen van Hunen, March 2017
# rad235U.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     = np.zeros(nt);            # Array to store 235U/204Pb values in
N_ini = 5.5                      # Initial 235U/204Pb
N[0]  = N_ini                    # Fill 1st position in array with N_ini
tau   = 0.704                    # halflife in Gyrs
lam   = np.log(2)/tau            # decay constant = ln(2)/tau 

for it in range(1,nt):   
    N[it] = N[it-1] * (1-dt*lam) # Numerical solution using timestepping

plt.figure(1)                    # Plot the solution
plt.clf()
plt.plot(time, N, 'o-', label='numerical')
plt.legend()
plt.xlabel('t(Gyrs)')
plt.ylabel('235U / 204Pb')
plt.show()

1.2) The analytical solution to the DE for radioactive decay is given by N=Niniexp(−λt)N=Niniexp⁡(−λt). Add the analytical solution from the lecture notes to the code and plot.

1.3. How does the numerical model perform?

1.4. How does the size of the timestep affect your results? Try different timesteps to check that the numerical solution approaches the analytical one for decreasing timesteps. This is called a ‘resolution test’, which is an important test for all codes you will be building in the future.

Session 1: Radioactive Decay, Part 2

1.5) Add the Backward Euler method to your existing code, and test it again.

1.6) Add the Crank-Nicholson method to your existing code, and test it again.

1.7) (Extra) If time permits, calculate the error between the analytical and numerical solutions for ²³⁵U for each of the three methods. Plot the error through time (if necessary in a separate subplot). Which of the three methods has the smallest error for a given size timestep and why?

A model solution for this session is available here.