Decay chain differential equations: Solutions through matrix analysis

doi:10.1016/j.cpc.2018.07.011

Computer Physics Communications

Volume 234, January 2019, Pages 188-194

https://doi.org/10.1016/j.cpc.2018.07.011 Get rights and content

Abstract

This paper demonstrates how the solutions to conventional radioactive decay equations can be derived using results from matrix analysis. In particular, we draw on results for the matrix exponential function when the matrix is triangular. By applying key theorems, the paper explains how the solutions to these equations can be presented in algorithmic form and in terms of divided differences of the exponential function. Furthermore, with little additional effort, the approach yields solutions to more general variations of these equations.

Introduction

The production rates of a network of nuclides undergoing spontaneous decay is described by the simultaneous linear differential equation system introduced by Bateman [1]. When nuclear reactions are induced in a neutron flux, Rubinson [2] provided a modification to include transformation by neutron absorption in addition to spontaneous decay. Although this model and its solution are well known through applying the methods of Laplace transforms, integrating factors or matrix algebra, this paper shows that there is much to be gained by taking an alternative, and perhaps more natural, approach using the matrix exponential function.

Following a description of the decay chain differential equations we introduce the matrix exponential function. The subsequent sections demonstrate the easy discovery of the Bateman solution and how important extensions to the basic model may be evaluated using this approach.

Section snippets

The decay chain equations

We begin with the description of the basic model as presented by Amaku et al. [3]. A serial decay chain of $n$ nuclides is one where the $i$ th nuclide of the chain decays to the $(i + 1)$ th nuclide of the chain. We denote the quantity of the $i$ th nuclide at time $t$ by $N_{i} (t)$ and its decay constant by $λ_{i}$ ( $\sec^{- 1}$ ). If we allow for nuclei immersed in a constant neutron flux $ϕ$ (neutrons ${cm}^{- 2} \sec^{- 1}$ ) and total neutron reaction cross section $σ_{i}$ ( ${cm}^{2}$ ) for the $i$ th nuclide, then the differential equation system may

The matrix exponential function

Apostol [6] provides a helpful introduction to the development of the matrix exponential function and its use in solving systems of differential equations with constant coefficients.

For a real number $t$ and fixed $n \times n$ matrix $A$ , with typical element $a_{i j}$ , let $S_{m} (t)$ denote the partial sum $S_{m} (t) = I + A t + \frac{1}{2!} A^{2} t^{2} + \frac{1}{3!} A^{3} t^{3} + \dots + \frac{1}{m!} A^{m} t^{m} .$ Then, subject to a suitable definition for the $norm$ of a matrix, e.g. $‖ A ‖ = \sum_{i, j} | a_{i j} |$ , it can be shown (see [6], Sec. 7.5) that the sequence of matrices { $S_{1} (t), S_{2} (t), \dots$ } has a

Solving the decay equations

In this section, we use the matrix exponential function to solve the system (3) and exploit its triangular form. We show how an explicit expression for the general solution can be derived in an economic manner. An earlier paper that applied this approach is Attaya [7].

As (3) is a homogeneous system with constant coefficients, we may write the unique solution to (3) as $N (t) = e^{A t} N_{0} .$ (See [6], Theorem 7.7.) Differentiating (5) with respect to $t$ gives $N^{'} (t) = A N (t)$ and noting (5) yields $N (0) = N_{0}$ at $t = 0$

Divided differences

We begin by introducing Newton’s divided differences. These are typically encountered in classical approximation theory and provide the coefficients for a polynomial interpolating a given function $f (.)$ whose value is known at a finite set of data points. Most textbooks in numerical analysis will provide a thorough coverage of this concept and one good choice is Atkinson [12].

For a function $f (x)$ defined at points $x_{0}, \dots, x_{k}$ , the $k$ th-order divided difference can be defined by the recurrence relation

Explicit solutions and the Bateman solution

Using the divided differencedefinition in (7), with $f (x) = e^{x t}$ , the solutions for $f_{i j}$ in (6) for our 3-membered branching chain example can now be expressed as follows: $\begin{matrix} f_{21} & = b_{21} f [- κ_{1}, - κ_{2}] \\ f_{32} & = b_{32} f [- κ_{2}, - κ_{3}] \\ f_{31} & = b_{31} f [- κ_{1}, - κ_{3}] + b_{32} b_{21} f [- κ_{1}, - κ_{2}, - κ_{3}] . \end{matrix}$ Extending the solution for a 4-membered chain is straightforward. The additional needed terms may be deduced as follows: $\begin{matrix} f_{43} & = b_{43} f [- κ_{3}, - κ_{4}] \\ f_{42} & = b_{42} f [- κ_{2}, - κ_{4}] + b_{43} b_{32} f [- κ_{2}, - κ_{3}, - κ_{4}] \\ f_{41} & = b_{41} f [- κ_{1}, - κ_{4}] + b_{42} b_{21} f [- κ_{1}, - κ_{2}, - κ_{4}] \\ + b_{43} b_{31} f [- κ_{1}, - κ_{3}, - κ_{4}] \\ + b_{43} b_{32} b_{21} f [- κ_{1}, - κ_{2}, - κ_{3}, - \end{matrix}$

Including an independent production source

Consider the problem of solving the non-homogeneous system $N^{'} (t) = A N (t) + q (t), N (0) = N_{0}$ where $q (t)$ is a $n \times 1$ vector function with $i$ th element representing an independent rate of production for the $i$ th nuclide. Using the formula for the derivative of $e^{A t}$ and the chain rule for differentiation, observe that $\begin{matrix} \frac{d}{d t} (e^{- A t} N (t)) & = e^{- A t} N^{'} (t) - A e^{- A t} N (t) \\ = e^{- A t} [N^{'} (t) - A N (t)], \end{matrix}$ making use of the fact that $e^{- A t}$ and $A$ commute. Hence if we premultiply (12) by $e^{- A t}$ (and rearrange), we can write $\frac{d}{d t} (e^{- A t} N (t)) = e^{- A t} q (t)$ and so $e^{- A}$

Repeats in the decay constants

In the instance when one or more of the decay constants are repeated, then the expression for $f_{i j}$ in Theorem 1 when $i > j$ will involve division by zero for one or more terms. Theorem 1 can be modified to cope with repeated eigenvalues (i.e. repeats in the diagonal terms as these are the eigenvalues for a triangular matrix) by solving the system considered as being in block triangular form. However, some care is needed if the blocks themselves have eigenvalues in common (see Parlett [16] and Davis

Conclusion

We have shown how the matrix exponential function can be used to develop an economic and more natural approach to solving decay chain differential equation systems. The approach not only delivers a concise expression for the well known Bateman solution but also informs how these solutions should be modified when one or more decay constants are repeated. It is also shown how richer forms of decay activities, such as branching and inclusion of independent production sources, can be easily handled

Acknowledgment

The author thanks Martin Owen for his comments on an earlier draft of this paper.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflict of interest

The author declares that there is no conflict of interest in this paper.

References (19)

AmakuM. et al.
Comput. Phys. Comm.
(2010)
AttayaH.
Fusion Eng. Des.
(1995)
DavisC.
Linear Algebra Appl.
(1973)
ParlettB.N.
Linear Algebra Appl.
(1976)
DreherR.
Ann. Nucl. Energy
(2013)
BatemanH.
Proc. Camb. Phil. Soc. Math. Phys. Sci.
(1910)
RubinsonW.
J. Chem. Phys.
(1949)
PressyanovD.S.
Amer. J. Phys.
(2002)
MoralL. et al.
Amer. J. Phys.
(2003)

There are more references available in the full text version of this article.

Cited by (0)

View full text