1. Introduction
The system of linear, constant real coefficient ordinary differential equations, referred to as the linear dynamic system in this article, appears in various disciplines. The linear dynamic systems can describe reactor depletion and radioactive decay dynamics and many related literature results concerning the radioactive decay linear dynamic system are available spanning over a century [18].
Many nuclear engineering and science applications require numerous iterations of radioactive decay calculations. An example requiring numerous iterative radioactive decay calculations is molten salt reactors (MSR) [9,10] where the reactor discharges a part of the molten salt inventory to the used fuel processor (UFP) for chemical treatment. The UFP treats the discharged used fuel to separate the fission products and feeds the processed fuel back to the reactor along with the exogenous fuel feed to make up the discharged salt from the reactor. The performance of the MSR system requiring online used fuel processing is tightly coupled with the performance of the UFP as the daily makeup fuel feed comes from the used fuel process operation. The evaluation of the MSR system performance requires numerous iterations of reactor depletion and used fuel process simulations. Another example is a pyroprocessing facility where radioactive materials are accumulated, divided, moved, and decayed daily [1114]. The evaluation of the pyroprocessing facility via modeling and simulation also involves tens of thousands of radioactive material decay calculations. In these applications, computationally efficient radioactive decay calculations are crucial to avoid a prolonged overall calculation time. The radioactive decay calculation involves matrix exponential that can be addressed with various numerical and mathematical techniques. A review on this topic with extensive references is available in [15]. The foci have been given to find accurate, robust, and efficient algorithms for matrix exponential.
The intent of this article is not in contributing to the prescribed aspects. Rather, specifically having the radioactive decay application in mind, the focus is on identifying the structural information of the graph induced from the radioactive decay dynamics in order to fragment the problem into the set of smaller size subproblems. The computational gain depends on the degree of the problem fragmentation and the radioactive decay application exhibits a high degree of fragmentation allowing a multiple order of magnitude computational improvement. To do so, section 2 gives relevant preliminaries regarding the linear dynamic systems and graph theories. Section 3 presents two graphbased algorithms for solving the linear dynamic system. An application to the radioactive decay calculation is detailed in Section 4. Section 5 gives concluding remarks.
2. Preliminaries
Consider the system of linear, constant real coefficient ordinary differential equations given below.
Denote the set of real numbers as R. With matrix and vector notations, a compact representation of the above equation set can be written as below.
where
The constraint is that the solution X(t) satisfies the initial state X_{0}, i.e., X(0) = X_{0}. Then, the solution is expressed with matrix multiplication of the initial state X_{0} and the matrix exponential e^{At} as below.
In this article, we refer matrix A as the rate matrix and e^{At} as the transition matrix.
A directed graph G is an ordered pair (V, E) where V is the set of vertices and E⊆V×V is the set of directed edges. A directed graph G = (V, E) is induced from rate matrix A where V is the set of vertices representing the index set of equation (1) and E is the set of edges reflecting the influence structure of rate matrix A. That is,
The set {V_{i}} is called the partition of a set V if
A path of length k in a directed graph G = (V, E) is the sequence of vertex {v_{0}, v_{1}, …, v_{k}} such that (v_{i}_{–1}, v_{i})∈E, for all i = 1, …, k. With this, we call that u∈V is reachable from v∈V if there exists a positive integer k such that {v_{0}, v_{1}, …, v_{k}} forms a path of length k in G = (V, E) with v = v_{0}, and u = v_{k}. Let’s define the reverse edge set of E as below.
We call that a directed graph G = (V, E) is weakly connected, for every pair of vertices u, v∈V, if u is reachable from v in G' = (V, E∪E^{–1}). A graph G_{s} = (V_{s}, E_{s}) is called a subgraph of G = (V, E) if V_{s}⊆V and E_{s}⊆E. A graph G_{s} = (V_{s}, E_{s}) is called a strict subgraph of G = (V, E) if V_{s}⊂V or E_{s}⊂E. A subgraph G_{a} = (V_{a}, E_{a}) of a directed graph G = (V, E) is called a weakly connected component of a directed graph G if G_{a} is weakly connected and there is no subgraph G_{b} of G such that G_{b} is weakly connected and G_{a} is a strict subgraph of G_{b}. With this, we define the weakly connected components of a directed graph G = (V, E) denoted by {G_{c,i} := (V_{c,i}, E_{c,i})} such that for all i, G_{c,i} is a weakly connected component and {V_{c,i}} and {E_{c,i}} are the partition of V and E, respectively.
Below gives an example illustrating the definitions presented above.
Example 1. Consider the set of linear differential equations below.
Rate matrix A in (3) can be constructed from the linear dynamic system in (5) as below.
Digraph G = (V, E) induced from rate matrix A is given in (6); Fig. 1 gives its graphical representation.
There are two weakly connected components: G_{c}_{,1} = ({1,2,3}, {(2,1),(3,2),(3,3)}) and G_{c}_{,2} =({4}, Ø).
The induced digraph facilitates an intuitive examination of the influence structure of the involved variables. Typically, the Depth First Search (DFS) over the undirected graph is used for finding the weakly connected components with O(V+E) complexity. A textbook treatise of an algorithm for finding the weakly connected components is available in [16].
3. Subdynamics Fragmentation
Denote the initial state imposed on the system given in Equation (1) by
Consider the following index set collecting the index of the nonzero elements in Equation (7). We call this the initially activated vertices.
The reach of i∈Î in the induced graph G = (V, E) can be defined as below.
The union of R_{i} gives the union of the reachable vertices from the initially activated vertices, denoted by R_{Î}. That is,
Then, we can construct the subgraph induced with R_{Î}, ${G}_{{R}_{\widehat{I}}}=\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$, where
The subgraph ${G}_{{R}_{\widehat{I}}}=\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$ of the induced graph G = (V, E) can be constructed with an application of DFS with O(V+E) complexity.
The set of unreachable vertices from Î, denoted by ${R}_{\widehat{I}}^{c}$, is the complement of R_{Î}. That is, ${R}_{\widehat{I}}^{c}\left\{1,\hspace{0.17em}\cdots ,N\right\}\backslash {R}_{\widehat{I}}$. By grouping the reachable and unreachable vertices, Equation (2) can be rewritten as below.
where ${X}_{{R}_{\widehat{I}}}\left(t\right)$ is the collection of the states corresponding to the reach of the initially activated vertices. On the other hand, ${X}_{{R}_{\widehat{I}}^{c}}\left(t\right)$ represents the remaining vertices that are not in the reach of the initially activated vertices. Note that
and
As ${R}_{\widehat{I}}^{c}$ is not reachable from R_{Î}, we have ${A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}^{c}}=0$. Then, Equation (9) becomes
Thus, the unreachable part of the linear dynamic system becomes
With ${X}_{{R}_{\widehat{I}}^{c}}\left(t\right)=0$, Equation (8) is further simplified to
The solution of Equation (10) is
Given the subgraph ${G}_{{R}_{\widehat{I}}}=\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$, building the rate matrix ${A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}}$ takes O(R_{Î}^{2}). The matrix exponential ${e}^{{A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}}t}$ can be calculated with O(R_{Î}^{3}) [15]. The multiplication, ${e}^{{A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}}t}{X}_{{R}_{\widehat{I}}}\left(0\right)$, can be done in O(R_{Î}^{2.3737}) with the algorithm given in [17]. Thus, overall complexity is governed by the matrix exponential computation step and O(R_{Î}^{3}).
3.1 Algorithm with Weakly Connected Component (WCC)
Evaluating the activated part of the system with graph reachability analysis before deriving the solution provides an opportunity to reduce the size of the problem as in Equation (10). In this subsection, for the problem formulated in Equation (10), we attempt to divide the problem further into smaller subproblems with the application of the weakly connected components graph analysis described in Section 2, which can be performed with $O\left(\left{R}_{\widehat{I}}\right+\left{E}_{{R}_{\widehat{I}}}\right\right)$ complexity. In a rough sense, two vertices of a directed graph, G, are weakly connected if a path can be found from one to the other in the undirected graph created by ignoring the direction of the edges in the directed graph, G.
First, denote the i^{th} weakly connected component of ${G}_{{R}_{\widehat{I}}}\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$ by G_{c,i} := (V_{c,i}, E_{i}) where V_{c,i}⊆R_{Î} and E_{c,i} is the edges defined over V_{c,i}. That is,
Suppose that there are M weakly connected components of G_{Î}. By rearranging the index, one can get the following block diagonal rate matrix where rate matrix A_{c,i} is constructed from the connected component G_{c,i} := (V_{c,i}, E_{i}). We note that the value of offdiagonals elements representing the connection between the weakly connected components is zero. That is,
With the above rate matrix, the corresponding linear dynamic system can be described as below.
Then, for i = 1, …, M,
Rather than solving the linear system with Equations (4) or (11) in a monolithic manner, the solution can be derived for each diagonal block induced from the weakly connected components.
Note that synthesizing rate matrix A_{c,i} takes O(G_{i}^{2}) and building all diagonal block matrices takes ${\sum}_{i}^{M}O\left({\left{G}_{i}\right}^{2}\right)$. As $\left{R}_{\widehat{I}}\right={\displaystyle {\sum}_{i}^{M}\left{G}_{i}\right}$, one can see that
Equation (13) implies that the overall complexity of fragmenting the whole dynamics into the weakly connected component based subdynamics is not worse than building the rate matrix ${A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}}$ from ${G}_{{R}_{\widehat{I}}}=\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$ Also, the overall complexity of building matrix exponentials of all subdynamics is ${\sum}_{i}^{M}O\left({\left{G}_{i}\right}^{3}\right)$, which is not worse than O(R_{Î}^{3}), the complexity of computing matrix exponential ${e}^{{A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}}t}$. That is,
Other miscellaneous computation steps such as matrix multiplications and combining subdynamics solution can be performed with complexities less than O(R_{Î}^{3}). Thus, the proposed solution strategy with the weakly connected component algorithm takes O(R_{Î}^{3}). For brevity, we shall refer to the algorithm described in this section as WCC.
3.2 Algorithm with Minimum Starting Vertices (MSV)
The section assumes that ${G}_{{R}_{\widehat{I}}}=\left({R}_{\widehat{I}},{E}_{{R}_{\widehat{I}}}\right)$ is acyclic. The minimum starting vertex set of ${G}_{{R}_{\widehat{I}}}$ is the collection of vertices that do not have incoming edges. That is,
As ${G}_{{R}_{\widehat{I}}}$ is acyclic, all vertices in R_{Î} should be reachable from the minimum starting vertex set. The reachable vertices from the initially activated vertices Î can be expressed with $S\left({G}_{{R}_{\widehat{I}}}\right)$ as below.
For $i\in S\left({G}_{{R}_{\widehat{I}}}\right)$, one can construct the disjoint set of the reachable vertices as below.
Then, $\left\{{R}_{i}^{}\hspace{0.17em}:\hspace{0.17em}i\in S\left({G}_{{R}_{\widehat{I}}}\right)\right\}$ becomes a partition of R_{Î}. That is,
For A⊆{1, …, N}, the projection map P_{A} : R^{N} → R^{N}, is defined as:
where
Denote the complement of R_{Î} by ${R}_{\widehat{I}}^{c}\hspace{0.17em}:=\hspace{0.17em}\left\{1,\hspace{0.17em}\cdots ,\hspace{0.17em}N\right\}\backslash {R}_{\widehat{I}}$. As the nonzero initial condition index Î and ${R}_{\widehat{I}}^{c}$ are disjoint, we have
Note that, with any partition {V_{1}, …, V_{M}} of R_{Î}, $\left\{{V}_{1},\hspace{0.17em}\cdots ,{V}_{M},{R}_{\widehat{I}}^{c}\right\}$ becomes the partition of the index set {1, …, N}. The initial condition can then be expressed with the summation of the projected initial values with the partition $\left\{{V}_{1},\hspace{0.17em}\cdots ,{V}_{M},{R}_{\widehat{I}}^{c}\right\}$ of R_{Î}. That is,
Note that X(t) can be expressed as the sum of respective state evolution with the projected initial state X_{0} induced from any partition {V_{i}} of the index set V. That is,
As $\left\{{R}_{i}^{}\hspace{0.17em}:\hspace{0.17em}i\in S\left({G}_{{R}_{\widehat{I}}}\right)\right\}$ is a partition of R_{Î}, $\left\{{R}_{i}^{}\hspace{0.17em}:\hspace{0.17em}i\in S\left({G}_{{R}_{\widehat{I}}}\right)\right\}$ and ${R}_{\widehat{I}}^{c}$ become the partition of the index set V. Then, we have
As ${P}_{{R}_{\widehat{I}}^{c}}\left(X\left(0\right)\right)=0$, we have
Let us focus on the following dynamics below.
The solution of the dynamics given in Equation (15) is
The set of unreachable vertices of G from $i\in S\left({G}_{{R}_{\widehat{I}}}\right)$, denoted by ${R}_{\widehat{I}}^{c}$, is the complement of R_{Î}, that is, ${R}_{\widehat{I}}^{c}\left\{1,\hspace{0.17em}\cdots ,N\right\}\backslash {R}_{\widehat{I}}$. Note that ${R}_{i}^{}\subseteq {R}_{\widehat{I}}$. Thus, we have ${R}_{i}^{}\cap {R}_{\widehat{I}}^{c}=\varnothing $ and all components of X(0) belonging to ${R}_{\widehat{I}}^{c}$ become zero with the projection map ${P}_{{R}_{\widehat{I}}^{}}$. Grouping the reachable and unreachable vertices and rewriting Equation (15) result in
With ${A}_{{R}_{\widehat{I}}\to {R}_{\widehat{I}}^{c}}=0$, one can see that
Then, the solution of the unreachable part of the dynamics becomes for all t ≥ 0,
The dynamics of the reachable part is
The solution of the reachable part dynamics becomes
The above derivation implies that the state evolution activated with the initial state X_{0} can be expressed by combining subdynamics solutions described in Equation (16) that are governed by the reachable vertices from a given starting vertex and an appropriately partitioned initial value. For brevity, we shall refer to the algorithm described in this section as MSV. Below gives an example that illustrates the partition concepts described in this section.
Example 2. Consider the digraph G = (V, E) in Fig. 2. The shaded vertices in the figure are the starting vertices of the graph G. The reachable vertices from each starting vertex are R_{1} and R_{5}. The partition of the vertex set V is {{1,2,3,4,6},{5}}. The corresponding initial value partition are ${P}_{{\overline{R}}_{1}}\left({X}_{0}\right)$ and ${P}_{{\overline{R}}_{5}}\left({X}_{0}\right)$.
4. Application to Decay Chain Graph
The set of differential equations representing radioactive decay with no external neutron flux is shown in Equation (17). For k = 1, …, N
In Equation (17), x_{k}(t) represents the amount of the k^{th} nuclide at decayed to time t and ẋ_{k} is its time derivative. λ_{k} is the decay constants for the k^{th} nuclide and f_{k,m} is the moles of the k^{th} nuclide created from the decay of one mole of the m^{th} nuclide, i.e., branching fraction. Then, the rate matrix can be constructed as below:
Rate matrix A is constructed with Idaho National Laboratory Mass Tracking System Database built from ENDF/ BVII database [18]. Decay digraph G = (V, E) that is induced from rate matrix A has the following properties:

V = N = 3631 vertices

E = 2367 transitions

2063 weakly connected components

Largest weakly connected component: 1299 vertices (mostly actinides)

Smallest component size: 1

Average component size: 1.76

Acyclic
Fig. 3 below shows decay digraph G = (V, E) clustering weakly connected components. Fig. 4 gives the number of reachable vertices from each vertex. Contrary to the size of the largest weakly connected component (1299 vertices), the size of the largest reachable vertex set is substantially small (36 vertices). The transition structure of the graph allures that an algorithm building the subdynamics with the minimum starting vertices is likely to be efficient. Fig. 5 shows the largest decay chain directed graph with 36 vertices (isotopes). The diagram includes neutron as a vertex as well.
MATLAB calculations were conducted to compare the performance of the algorithms described in this article.
The monolithic calculation is performed with expm MATLAB function where a scaling and squaring algorithm with a Pade approximation [15,19,20] is used. WCC and MSV algorithms also use the same expm function for their subsystem transition matrix computation.
For the above numerical calculations, the following resources are used:
Table 1 and Fig. 6 give the average computational time of the calculation methods described in this article for 100 trials. Given the number of nonzero initial isotopes to be simulated, selection of isotopes was randomized and given the value of 1. The values in the table are relative to the monolithic calculation time with 10 random nonzero initial isotopes. One can observe that the monolithic calculation ignoring the reachable vertices shows a consistent computation time. Once the reachable part of graph is considered, a significant reduction of computational time is observed for all algorithms. All three methods accounting the reachable part of the graph show a gradual increase in computational time as the number of reachable vertices increases. MSV algorithm shows the least steep increase of computational time followed by WCC algorithm and expm with reachability consideration. The superiority of expm function with reachability consideration over WCC and MSV algorithms for the smaller size initially activated vertexes may come from the implementation inefficiency of WCC and MSV algorithms where multiple graph searches were done due to the adaptation of MATLAB generic functions, which can be consolidated with a customized graph search algorithm.
5. Concluding Remarks
An attempt to view the system of linear, constant real coefficient ordinary differential equations via a graph structure allows various graph algorithm treatments to enable efficient computation of the solution. Connected component analysis, reachability analysis, and the synthesis of the minimal starting vertex set were used to break the problem into the smaller size subproblems. From a theoretical perspective, the worstcase computational complexity does not exhibit an improvement. However, empirical results demonstrated with radioactive decay application show a significant saving of computational efforts. Furthermore, we note that the subproblems are independent and allow parallel computations.