Journal Search Engine
ISSN : 1738-1894(Print)
ISSN : 2288-5471(Online)
Journal of Nuclear Fuel Cycle and Waste Technology Vol.18 No.1 pp.19-29
DOI : https://doi.org/10.7733/jnfcwt.2020.18.1.19

# Efficient Computation of Radioactive Decay with Graph Algorithms

Tae-Sic Yoo*
Idaho National Laboratory, 775 MK Simpson Blvd, Idaho Falls, Idaho 83402, USA
Corresponding Author. Tae-Sic Yoo, Idaho National Laboratory, E-mail: tae-sic.yoo@inl.gov, Tel: +1-208-526-9394
December 4, 2019 February 27, 2020 March 25, 2020

## Abstract

This paper gives two graph-based algorithms for radioactive decay computation. The first algorithm identifies the connected components of the graph induced from the given radioactive decay dynamics to reduce the size of the problem. The solutions are derived over the precalculated connected components, respectively and independently. The second algorithm utilizes acyclic structure of radioactive decay dynamics. The algorithm evaluates the reachable vertices of the induced system graph from the initially activated vertices and finds the minimal set of starting vertices populating the entire reachable vertices. Then, the decay calculations are performed over the reachable vertices from the identified minimal starting vertices, respectively, with the partitioned initial value over the reachable vertices. Formal arguments are given to show that the proposed graph inspired divide and conquer calculation methods perform the intended radioactive decay calculation. Empirical efforts comparing the proposed radioactive decay calculation algorithms are presented.

## 초록

This is an Open-Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

## 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 [1-8].

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 make-up 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 [11-14]. 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 graph-based 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.

$x ˙ k = ∑ i = 1 N a k , i ⋅ x i ( t ) ; ​ for k = 1 to N$
(1)

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.

$X ˙ ( t ) = A ⋅ X ( t )$
(2)

where

$X ˙ ( t ) = : [ x ˙ 1 ( t ) x ˙ 2 ( t ) ⋮ x ˙ N ( t ) ] ∈ R N × 1 ; X ( t ) : = [ x 1 ( t ) x 2 ( t ) ⋮ x N ( t ) ] ∈ R N × 1 ; A : = [ a 1 , 1 a 1 , 2 ⋯ a 1 , N a 2 , 1 a 2 , 2 ⋯ a 2 , N ⋮ ⋮ ⋱ ⋮ a N , 1 a N , 2 ⋯ a N , N ] ∈ R ​ N × N$
(3)

The constraint is that the solution X(t) satisfies the initial state X0, i.e., X(0) = X0. Then, the solution is expressed with matrix multiplication of the initial state X0 and the matrix exponential eAt as below.

$X ( t ) = e A t ⋅ X 0$
(4)

In this article, we refer matrix A as the rate matrix and eAt as the transition matrix.

A directed graph G is an ordered pair (V, E) where V is the set of vertices and EV×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,

$V = { 1 , 2 , ⋯ , N } and ( i , j ) ∈ E ⊆ V × V ⇔ a j , i ≠ 0$

The set {Vi} is called the partition of a set V if

$[ V = ∪ i V i ] ∧ [ V i ∩ V j = ∅ i f i ≠ j ]$

A path of length k in a directed graph G = (V, E) is the sequence of vertex {v0, v1, …, vk} such that (vi–1, vi)∈E, for all i = 1, …, k. With this, we call that uV is reachable from vV if there exists a positive integer k such that {v0, v1, …, vk} forms a path of length k in G = (V, E) with v = v0, and u = vk. Let’s define the reverse edge set of E as below.

$E − 1 = : { ( ν , u ) : ( u , ν ) ∈ E }$

We call that a directed graph G = (V, E) is weakly connected, for every pair of vertices u, vV, if u is reachable from v in G' = (V, EE–1). A graph Gs = (Vs, Es) is called a subgraph of G = (V, E) if VsV and EsE. A graph Gs = (Vs, Es) is called a strict subgraph of G = (V, E) if VsV or EsE. A subgraph Ga = (Va, Ea) of a directed graph G = (V, E) is called a weakly connected component of a directed graph G if Ga is weakly connected and there is no subgraph Gb of G such that Gb is weakly connected and Ga is a strict subgraph of Gb. With this, we define the weakly connected components of a directed graph G = (V, E) denoted by {Gc,i := (Vc,i, Ec,i)} such that for all i, Gc,i is a weakly connected component and {Vc,i} and {Ec,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.

$x ˙ 1 = x 2 ; x ˙ 2 = x 3 ; x ˙ 3 = − x 3 ; x ˙ 4 = 0$
(5)

Rate matrix A in (3) can be constructed from the linear dynamic system in (5) as below.

$A = [ 0 1 0 0 0 0 1 0 0 0 − 1 0 0 0 0 0 ]$

Digraph G = (V, E) induced from rate matrix A is given in (6); Fig. 1 gives its graphical representation.

$V = { 1 , 2 , 3 , 4 } a n d E = { ( 2 , 1 ) , ( 3 , 2 ) , ( 3 , 3 ) }$
(6)

There are two weakly connected components: Gc,1 = ({1,2,3}, {(2,1),(3,2),(3,3)}) and Gc,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

$X 0 : = [ x 1 , 0 , x 2 , 0 ⋯ ∈ x N , 0 ] T$
(7)

Consider the following index set collecting the index of the nonzero elements in Equation (7). We call this the initially activated vertices.

$I ^ = { i : x i , 0 ≠ 0 }$

The reach of iÎ in the induced graph G = (V, E) can be defined as below.

$R i : = { ν ∈ V | ν is reachable from i ∈ I ^ }$

The union of Ri gives the union of the reachable vertices from the initially activated vertices, denoted by RÎ. That is,

$R I ^ = ∪ i ∈ I ^ R i$

Then, we can construct the subgraph induced with RÎ, $G R I ^ = ( R I ^ , E R I ^ )$, where

$E R I ^ = { ( i , j ) ∈ E : i , j ∈ R I ^ }$

The subgraph $G R I ^ = ( R I ^ , E R I ^ )$ 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 I ^ c$, is the complement of RÎ. That is, $R I ^ c { 1 , ⋯ , N } \ R I ^$. By grouping the reachable and unreachable vertices, Equation (2) can be rewritten as below.

$[ X ˙ R I ^ ( t ) X ˙ R I ^ c ( t ) ] = [ A R I ^ → R I ^ A R I ^ c → R I ^ A R I ^ → R I ^ c A R I ^ c → R I ^ c ] [ X R I ^ ( t ) X R I ^ c ( t ) ]$

where $X R I ^ ( t )$ is the collection of the states corresponding to the reach of the initially activated vertices. On the other hand, $X R I ^ c ( t )$ represents the remaining vertices that are not in the reach of the initially activated vertices. Note that

$X ˙ R I ^ ( t ) = A R I ^ → R I ^ X R I ^ ( t ) + A R I ^ c → R I ^ X R I ^ c ( t )$
(8)

and

$X ˙ R I ^ c ( t ) = A R I ^ → R I ^ c X R I ^ ( t ) + A R I ^ c → R I ^ c X R I ^ c ( t )$
(9)

As $R I ^ c$ is not reachable from RÎ, we have $A R I ^ → R I ^ c = 0$. Then, Equation (9) becomes

$X ˙ R I ^ c ( t ) = A R I ^ c → R I ^ c X R I ^ c ( t )$

Thus, the unreachable part of the linear dynamic system becomes

$X R I ^ c ( t ) = e A R I ^ c → R I ^ c t X R I ^ c ( 0 ) = e A R I ^ c → R I ^ c t 0 = 0$

With $X R I ^ c ( t ) = 0$, Equation (8) is further simplified to

$X ˙ R I ^ ( t ) = A R I ^ → R I ^ X R I ^ ( t )$
(10)

The solution of Equation (10) is

$X R I ^ ( t ) = e A R I ^ → R I ^ t ⋅ X R I ^ ( 0 )$
(11)

Given the subgraph $G R I ^ = ( R I ^ , E R I ^ )$, building the rate matrix $A R I ^ → R I ^$ takes O(|RÎ|2). The matrix exponential $e A R I ^ → R I ^ t$ can be calculated with O(|RÎ|3) [15]. The multiplication, $e A R I ^ → R I ^ t X R I ^ ( 0 )$, 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 ( | R I ^ | + | E R I ^ | )$ 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 ith weakly connected component of $G R I ^ ( R I ^ , E R I ^ )$ by Gc,i := (Vc,i, Ei) where Vc,iRÎ and Ec,i is the edges defined over Vc,i. That is,

$E c , i : = { ( u , ν ) ∈ E : u , ν ∈ V c , i }$

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 Ac,i is constructed from the connected component Gc,i := (Vc,i, Ei). We note that the value of off-diagonals elements representing the connection between the weakly connected components is zero. That is,

$A R I ^ = [ A c , 1 0 ⋯ 0 0 A c , 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ A c , M ]$

With the above rate matrix, the corresponding linear dynamic system can be described as below.

$[ X ˙ c , 1 X ˙ c , 2 ⋮ X ˙ c , M ] = [ A c , 1 0 ⋯ 0 0 A c , 2 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ A c , M ] [ X c , 1 X c , 2 ⋮ X c , M ]$

Then, for i = 1, …, M,

$X c , i ( t ) = e A c , i t X c , i ( 0 )$
(12)

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 Ac,i takes O(|Gi|2) and building all diagonal block matrices takes $∑ i M O ( | G i | 2 )$. As $| R I ^ | = ∑ i M | G i |$, one can see that

$∑ i M O ( | G i | 2 ) ≤ O ( | R I ^ | 2 )$
(13)

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 I ^ → R I ^$ from $G R I ^ = ( R I ^ , E R I ^ )$ Also, the overall complexity of building matrix exponentials of all subdynamics is $∑ i M O ( | G i | 3 )$, which is not worse than O(|RÎ|3), the complexity of computing matrix exponential $e A R I ^ → R I ^ t$. That is,

$∑ i M O ( | G i | 3 ) ≤ O ( | R I ^ | 3 )$
(14)

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 I ^ = ( R I ^ , E R I ^ )$ is acyclic. The minimum starting vertex set of $G R I ^$ is the collection of vertices that do not have incoming edges. That is,

$S ( G R I ^ ) = { ν ∈ R I ^ : ∀ u ∈ V \ { ν } , ∃ ( u , ν ) ∈ E R I ^ }$

As $G R 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 ( G R I ^ )$ as below.

$R I ^ = ∪ i ∈ S ( G R I ^ ) R i$

For $i ∈ S ( G R I ^ )$, one can construct the disjoint set of the reachable vertices as below.

$R i − = R i \ ∪ j < i R j$

Then, ${ R i − : i ∈ S ( G R I ^ ) }$ becomes a partition of RÎ. That is,

$[ R I ^ = ∪ i ∈ S ( G R I ^ ) R i − ] ∧ [ R i − ∩ R j − = ∅ if i ≠ j ]$

For A⊆{1, …, N}, the projection map PA : RNRN, is defined as:

$P A ( [ x 1 x 2 ⋯ x N ] T ) : = [ z 1 z 2 ⋯ z N ] T$

where

$z i = { x i i f i ∈ A 0 o . w .$

Denote the complement of RÎ by $R I ^ c : = { 1 , ⋯ , N } \ R I ^$. As the non-zero initial condition index Î and $R I ^ c$ are disjoint, we have

$P R I ^ c ( [ x 1 , 0 , x 2 , 0 ⋯ x N , 0 ] T ) = 0$

Note that, with any partition {V1, …, VM} of RÎ, ${ V 1 , ⋯ , V M , R I ^ c }$ 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 ${ V 1 , ⋯ , V M , R I ^ c }$ of RÎ. That is,

$X ( 0 ) = P R I ^ c ( X ( 0 ) ) + ∑ i = 1 M P V i ( X ( 0 ) ) = ∑ i = 1 M P V i ( X ( 0 ) )$

Note that X(t) can be expressed as the sum of respective state evolution with the projected initial state X0 induced from any partition {Vi} of the index set V. That is,

$X ( t ) = e A t ∑ i P V i ( X ( 0 ) ) = ∑ i e A t P V i ( X ( 0 ) )$

As ${ R i − : i ∈ S ( G R I ^ ) }$ is a partition of RÎ, ${ R i − : i ∈ S ( G R I ^ ) }$ and $R I ^ c$ become the partition of the index set V. Then, we have

$X ( 0 ) = P R I ^ c ( X ( 0 ) ) + ∑ i ∈ S ( G R I ^ ) P R i − ( X ( 0 ) )$

As $P R I ^ c ( X ( 0 ) ) = 0$, we have

$X ( t ) = ∑ i ∈ S ( G R I ^ ) e A t P R i − ( X ( 0 ) )$

Let us focus on the following dynamics below.

$X ˙ i ( t ) = A X i ( t ) ; X i ( 0 ) = P R i − ( X ( 0 ) )$
(15)

The solution of the dynamics given in Equation (15) is

$X i ( t ) = e A t P R i − ( X ( 0 ) )$

The set of unreachable vertices of G from $i ∈ S ( G R I ^ )$, denoted by $R I ^ c$, is the complement of RÎ, that is, $R I ^ c { 1 , ⋯ , N } \ R I ^$. Note that $R i − ⊆ R I ^$. Thus, we have $R i − ∩ R I ^ c = ∅$ and all components of X(0) belonging to $R I ^ c$ become zero with the projection map $P R I ^ −$. Grouping the reachable and unreachable vertices and rewriting Equation (15) result in

$[ X ˙ R i ( t ) X ˙ R i c ( t ) ] = [ A R i → R i A R i c → R i A R i → R i c A R i c → R i c ] [ X R i ( t ) X R i c ( t ) ] ; [ X R i ( 0 ) X R i c ( 0 ) ] = [ X R i ( 0 ) 0 ]$

With $A R I ^ → R I ^ c = 0$, one can see that

$X ˙ R i c ( t ) = A R i → R i c X R i ( t ) + A R i c → R i c X R i c ( t ) = A R i c → R i c X R i c ( t )$

Then, the solution of the unreachable part of the dynamics becomes for all t ≥ 0,

$X R i c ( t ) = e A R i c → R i c t X R i c ( 0 ) + e A R i c → R i c t 0 = 0$

The dynamics of the reachable part is

$X ˙ R i ( t ) = A R i → R i X R i ( t ) + A R i c → R i X R i c ( t ) = A R i → R i X R i ( t )$

The solution of the reachable part dynamics becomes

$X R i ( t ) = e A R i → R i t X R i ( 0 )$
(16)

The above derivation implies that the state evolution activated with the initial state X0 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 R1 and R5. The partition of the vertex set V is {{1,2,3,4,6},{5}}. The corresponding initial value partition are $P R ¯ 1 ( X 0 )$ and $P R ¯ 5 ( X 0 )$.

## 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

$x ˙ k = − λ k x k + ∑ m ∈ ( 1 , ⋯ , N ) \ { k } f k , m ⋅ λ m ⋅ x m$
(17)

In Equation (17), xk(t) represents the amount of the kth nuclide at decayed to time t and k is its time derivative. λk is the decay constants for the kth nuclide and fk,m is the moles of the kth nuclide created from the decay of one mole of the mth nuclide, i.e., branching fraction. Then, the rate matrix can be constructed as below:

$A = [ − λ 1 f 1 , 2 λ 2 ⋯ f 1 , N λ N f 2 , 1 λ 1 − λ 2 ⋯ f 2 , N λ N ⋮ ⋮ ⋱ ⋮ f N , 1 λ 1 f N , 2 λ 2 ⋯ − λ N ] ∈ R N × N$

Rate matrix A is constructed with Idaho National Laboratory Mass Tracking System Database built from ENDF/ B-VII 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:

• Processor: Intel® Core™ i7-6600U CPU @2.60GHz

• RAM: 16.0GB

• MATLAB Version 9.3 (R2017b)

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 worst-case 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.

## 6. Acknowledgements

Rate matrix A in Section 4 is adapted from the MTG decay rate matrix implemented by DeeEarl Vaden.

This submitted manuscript was authored by a contractor of the U.S. Government under DOE Contract No. DE-AC07-05ID14517. Accordingly, the U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for U.S. Government purposes.

This information was prepared as an account of work sponsored by an agency of the U.S. Government. Neither the U.S. Government nor any agency thereof, nor any of their employees, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness, of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. References herein to any specific commercial product, process, or service by trade name, trade mark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the U.S. Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the U.S. Government or any agency thereof.

## Figure

Induced digraph G = (V,E) from rate matrix A.

Illustrative example describing a value partition induced from the starting vertices.

Isotope decay chain graph (3631 isotopes).

Number of reachable isotopes from each isotope.

Largest decay chain (36 isotopes counting neutron as isotope).

Computation time comparison with respect to the number of activated initial variables.

## Table

Relative average computational time with varying number of nonzero initial variables

## Reference

1. H. Bateman, “The solution of a system of differential equations occurring in the theory of radio-active transformations”, Proc. Cambridge Phil. Soc., 1908 15, 423- 427 (1908).
2. D.R. Vondy, Development of a General Method of Explicit Solution to the Nuclide Chain Equations for Digital Machine Calculations (Thesis), No. ORNL-TM-361, Oak Ridge National Lab., Tenn., 1962.
3. M.J. Bell, ORIGEN: the ORNL isotope generation and depletion code. No. ORNL--4628. Oak Ridge National Lab., (1973).
4. J. Cetnar. “General solution of Bateman equations for nuclear transmutations”, Annals of Nuclear Energy, 33 (7), 640-645 (2006).
5. M. Pusa and J. Leppänen, “Computing the matrix exponential in burnup calculations”, Nuclear science and engineering, 164(2), 140-150 (2010).
6. A.E. Isotalo and P.A. Aarnio, “Comparison of depletion algorithms for large systems of nuclides”, Annals of Nuclear Energy, 38(2-3), 261-268 (2011).
7. R. Dreher, “Modified Bateman solution for identical eigenvalues”, Annals of Nuclear Energy, 53, 427-438 (2013).
8. D. Vaden and T. Yoo, “Radioactive Decay Computation with Dynamic Source Terms”, Nuclear Science and Engineering, 193(5), 549-553 (2019).
9. J. Serp, M. Allibert, O. Beneš, S. Delpech, O. Feynberg, V. Ghetta, D. Heuer, D. Holcomb, V. Ignatiev, J.L. Kloosterman, and L. Luzzi, “The molten salt reactor (MSR) in generation IV: overview and perspectives”, Progress in Nuclear Energy, 77, 308-319 (2014).
10. T.J. Dolan ed., Molten Salt Reactors and Thorium Energy, Woodhead Publishing (2017).
11. H. Lee, G. Park, K. Kang, J. Hur, J. Kim, D. Ahn, Y. Cho, and E. Kim, “Pyroprocessing technology development at KAERI”, Nuclear engineering and technology, 43(4), 317-328 (2011).
12. T. Inoue, T. Koyama, and Y. Arai, “State of the art of pyroprocessing technology in Japan”, Energy Procedia, 7, 405-413 (2011).
13. T.Y. Karlsson, G.L. Fredrickson, T. Yoo, D. Vaden, M.N. Patterson, and V. Utgikar, “Thermal analysis of projected molten salt compositions during FFTF and EBR-II used nuclear fuel processing”, Journal of Nuclear Materials, 520, 87-95 (2019).
14. T. Yoo and D. Vaden, “A new inventory tracking method for Mark-V electrorefiner”, Annals of Nuclear Energy, 128, 406-413 (2019).
15. C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later”, SIAM review, 45(1), 3-49 (2003).
16. T.H. Cormen, C.E. Leiserson, R.L. Rivest, and C. Stein, Introduction to algorithms, MIT press (2009).
17. D. Coppersmith and S. Winograd, “Matrix multiplication via arithmetic progressions”, Journal of symbolic computation, 9(3), 251-280 (1990).
18. M. B. Chadwick, M. Herman, P. Obložinský, M.E. Dunn, Y. Danon, A.C. Kahler, D.L. Smith, B. Pritychenko, G. Arbanas, R. Arcilla, and R. Brewer, “ENDF/B-VII. 1 nuclear data for science and technology: cross sections, covariances, fission product yields and decay data”, Nuclear data sheets, 112(12), 2887- 2996 (2011).
19. H. Padé, Sur la représentation approchée d'une fonction par des fractions rationnelles, No. 740, Gauthier- Villars et fils, 1892.
20. R.C. Ward, “Numerical computation of the matrix exponential with accuracy estimate”, SIAM Journal on Numerical Analysis, 14(4), 600-610 (1977).
1. SEARCH
2. #### JNFCWT

Online Submission

4. #### Editorial OfficeContact Information

- Tel: +82-42-861-5851, 866-4205
- Fax: +82-42-861-5852
- E-mail: krs@krs.or.kr