half-life diff eqn euler

Relation between half-life and time step

Some examples of first order differential equations are mechanical systems, electrical circuits, population models, Newton’s law of cooling, and compartemental analysis ¹. Solution of a differential equation can be approximated using many methods, where one of the oldest and easiest to use is Euler’s method ². As analytical solution, the exponential function is used for investigating temperature of heated objects, population growth and radioactive decay ³. In radioactive decay a useful measurement point is the time take for half of the initial amount of material do decay, which is known as half-life ⁴. Certain model of population growth and radioactive decay share similar differential equations, which also have similar solution.

When solving differential equation of radioactive decay with numerical approach, which is Euler method ini this case, the half-life will depend on the chosen time step. In order to obtain the right time step, a coefficient is required, where its value relates to time step value.

Suppose that a half-life $T_\frac12$ is defined, than the decay constant would be

$$\tag{1} \lambda = \frac{\ln 2}{T_\frac12} $$

from half-life and decay constant relation.

Since Euler method is derived from the simple forward difference expression ⁵, following relation between half-time and time step

$$\tag{2} T_\frac12 = \left( \frac{\sqrt[n]{2} \ln 2}{\sqrt[n]{2} - 1} \right) \Delta t. $$

can be used. Substitute Eqn (2) to Eqn (1) will give

$$\tag{3} \lambda = \frac{\sqrt[n]{2} - 1}{\sqrt[n]{2} \Delta t}, $$

which holds for $T_\frac12 = n \Delta t$ with $n = 1, 2, 3, ..$.

Using $N_0 = 1600$, $\Delta t = 0.1$, $n = 10$, total time $t_{\rm tot} = M \Delta t$ with $M = 20$ following results are obtained.

Table 1. Decay of entity $N$ from $t = 0$ to $t = 2$ with $T_\frac12 = 1$.

Results in Table 1 can be obtained using following code

N = 1600
dt = 0.1

n = 10
Thalf = n * dt

rn2 = 2 ** (1/n)
lam = (rn2 - 1) / (rn2 * dt)

M = 21
for i in range(M):
  Ns = f'{N:.2f}'
  ts = f'{i*dt:.2f}'
  print(ts, Ns)
  
  N = (1 - lam * dt) * N

which is available on OneCompiler 42uaqagv7 .

Eqn (3) can be further written in the form of

$$\tag{4} \lambda = \frac{2^{\Delta t / T_\frac12} - 1}{2^{\Delta t / T_\frac12} \Delta t}, $$

which is more applicable, especially when half-time $T_\frac12$ is a non-integer multiple of time step $\Delta t$.