Half-Life Calculation of Radioactive Nuclei

Topic: Radioactive Decay

Subject: Physics

Tool: Scilab

By: Gani Comia, Jan. 2025

Half-life of the radioactive substance is defined as the time it takes for half of the radioactive nuclei in a sample to decay. This can be calculated from the initial value condition (IVC) of the ordinary differential equation (ODE) that describes radioactive decay as shown in Equation (1).

$$\frac{dN}{dt}\;=\;-\lambda\;N \tag{1}$$

Where:

\(N\) - amount of radioactive nuclei at time \(t\),

\(\lambda\) - decay constant of a nuclei in per unit time,

And the initial conditions are,

$$\left(\;t_0 = 0 \;,\; N_0 = 1 \right) \tag{2}$$

This article will present the calculation of the half-life and plot the decay profile of the sample nuclei \(^{226}Ra_{88}\) for illustration of the solution to ODE as shown in Figure (1) using Scilab scripting.

Figure 1. Radioactive Decay and Half-Life of \(^{226}Ra_{88}\).

Scilab Script:

The solution will be based on the IVC of the ODE representing radioactive decay. Its ODE from Equation (1) can be written in Scilab script.

clear;clc;
// ordinary differential equation with initial value condition
// dN/dt = f(t,N) with IVC N(0)=1.

// (1) differential equation for radioactive decay
function dNdt=f(t, N)
    // lambda - per year, decay constant of radioactive substance
    lambda = (1.4e-11)*(60*60*24*365); // for 226 Ra_88 
    dNdt = -lambda.*N;
endfunction

Assume an initial amount of nuclei as \(N_0 = 1\) at the initial time, \(t_0 = 0\). If the amount of radioactive nuclei reached \(N_h = 0.5\), which is \(\frac{1}{2}\) of the initial condition, then that is the time the given amount reached its half-life. Our simulation will be run from 0 to 5000 years, as we still do not know the exact time of its half-life. Here is the code snippet for the initial condition and time domain.

// (2) initial condition and time domain
t0 = 0;
N0 = 1;     // unit, arbitrary amount, 
            // if N = 0.5, it reaches the substance's half-life
t = 0:5000; // years, time domain

Scilab has a built-in function, \(ode()\), to provide a solution to first-order IVC for ODEs.

// (3) ODE solution using the Scilab ode() function
N = ode(N0,t0,t,f);

Plotting of the solution for the amount of nuclei, \(N\), for the given time domain can be done using the below script.

// (4) plotting of solution
clf;
plot(t, N,"b-","linewidth",3)
title("$\large\text{Radioactive Decay of}\;{Ra_{88}^{226}}$")
xlabel("Time [ years ]","fontsize",3)
ylabel("Amount of Nuclei [ unit ]","fontsize",3)
legend("$\LARGE\frac{dN}{dt}\;=\;-\,\lambda\,N$",with_box = %F)
xstring(3000,0.8,"https://gani-mech-toolbox.blogspot.com")
xgrid(3,1)
Solving for the half-life, Th, for the given amount, Nh = 0.5, can be done with the use of the Scilab interpolation function, interp1().
// (5) interpolation to find the half-life, Th, for Nh = 0.5 unit
Nh = 0.5;
Th = interp1(N,t,Nh);
disp(Th)
The information about the half-life will be plotted using the following script. 
// (6) plotting of half-life lines
verX = [Th Th];
verY = [0 Nh];
plot(verX, verY, "r--","linewidth",1.5)
horX = [0 Th];
horY = [Nh Nh];
plot(horX, horY, "r--","linewidth",1.5)
plot(Th,Nh,"ro");

format(7);
note1 = ["T (yrs) = "+string(Th)+"; N (unit) = "+string(Nh)];
note2 = ["$\Large \textbf{Half - Life}$"];
xstring(Th,Nh,note1)
xstring(Th,Nh+0.05,note2)
a = gca();
a.children.children.mark_background = 5;
The Scilab script provided here can be used for other substances with different decay constants to calculate the half-life and decay profile.
Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the script file.