\todo{Anthony will need to review this text, or rewrite it if necessary.}

In addition to numeric tests, it is possible to run a sanity check on the CRAM

method based on physical tests. If $A$ is a transmutation matrix with no

fission, $e^{-At}$ represents {\color{red}\ldots}.

The $i$th column of $e^{-At}=e^{-At}I$ is equal to $e^{-At}\delta_{i}$, where

$\delta_i$ is a column vector with $1$ in the $i$th row and $0$ elsewhere.

That is, the $i$th column of $e^{-At}$ represents the transmutation of a unit

mass of the $i$th nuclide after $t$ seconds with no fission. By the

conservation of mass, these columns should each sum to 1, the mass that was

started. While pure transmutation of an isolated unit mass of each nuclide may

not be physically feasible, this provides a mathematical sanity check for

different methods of computing $e^{-At}$.

Mathematically, $A$ is a matrix where each non-diagonal entry is nonnegative,

and each diagonal entry is the negative of the sum of the non-diagonal entries

in its column, so that the columns of $A$ each sum to $0$. Such matrices are

called \textit{rate matrices}~\cite{glasser1980properties}, and it can be

easily shown that if $A$ is a rate matrix, then $e^{-At}$ is Markov, that is,

its columns sum to 1.\footnote{If $A$ is a rate matrix,

$[1 \cdots 1] A = [0 \cdots 0]$. For $n\geq 1$, we have $[1 \ldots 1] A^n = [0 \cdots 0] A^{n-1} = [0\cdots 0]$. Using

Equation~\ref{eq:matrix-exponential},

[1\ldots 1] + [0 \ldots 0] + [0 \ldots 0] + \cdots = [1\ldots 1]$. That is,

the columns of $e^{-At}$ each sum to $1$.}

Figures~\ref{fig:nofission-pwru50-1-day}, \ref{fig:nofission-pwru50-1-year},

\ref{fig:nofission-pwru50-1000-years},

and~\ref{fig:nofission-pwru50-1-million-years} show the results of this sanity

check for various solvers and time steps. The matrix $A$ was generated from

the ORIGEN data library \texttt{pwru50.lib}, with fission product yields

omitted so that its columns sum to $0$. The figures show histograms of the

column sums of $e^{-At}$ minus 1 for $t$ equal to 1 day, 1 year, 1000 years,

and 1 million years. A perfect calculation would have every column sum minus 1

equal exactly zero.

The leftmost histograms show

\texttt{scipy.\allowbreak{}sparse.\allowbreak{}linalg.\allowbreak{}expm} for

comparison. The remaining show CRAM with degree 14 using

\texttt{sympy.\allowbreak{}lambdify} using

\texttt{scipy.\allowbreak{}sparse.\allowbreak{}linalg} with the UMFPACK and

SuperLU backends, and a C solver generated by transmutagen against the

sparsity pattern of $A$. Note that for CRAM with degree 14, the absolute error

is expected to be on the order of $10^{-14}$. These figures show that CRAM is

highly accurate even for large time steps. The \texttt{expm} implementation in

\texttt{scipy.\allowbreak{}sparse} is inaccurate even for small time steps,

and is extremely inaccurate for large time steps. At

$t = 1 \mathrm{\ million\ years}$, it gave a resulting matrix with NaN values.

Finally, CRAM with the UMFPACK backend has issues with a handful of nuclides

at small timesteps and with most nuclides at 1 million years, whereas SuperLU

and the transmutagen generated C solver backends have nearly identical

precision characteristics, staying accurate for all nuclides even at large

time steps.

\begin{figure}[!ht]

\resizebox{0.9\textwidth}{!}{\input{nofission-pwru50-1-day.pgf}}

\caption{Sanity check for 1 day.}

\label{fig:nofission-pwru50-1-day}

\end{figure}

\begin{figure}[!ht]

\resizebox{0.9\textwidth}{!}{\input{nofission-pwru50-1-year.pgf}}

\caption{Sanity check for 1 year.}

\label{fig:nofission-pwru50-1-year}

\end{figure}

\begin{figure}[!ht]

\resizebox{0.9\textwidth}{!}{\input{nofission-pwru50-1000-years.pgf}}

\caption{Sanity check for 1000 years.}

\label{fig:nofission-pwru50-1000-years}

\end{figure}

\begin{figure}[!ht]

\resizebox{0.9\textwidth}{!}{\input{nofission-pwru50-1-million-years.pgf}}

\caption{Sanity check for 1 million years.}

\label{fig:nofission-pwru50-1-million-years}

\end{figure}