Higher-order Chebyshev Rational Approximation Method (CRAM)

The burnup equations can in principle be solved by computing the exponential of the burnup matrix. However, due to the difficult numerical characteristics of burnup matrices, the problem is extremely stiff, and the matrix exponential solution was long considered infeasible for an entire burnup system containing over a thousand nuclides. After discovering that the eigenvalues of burnup matrices are generally confined to a region near the negative real axis, the Chebyshev rational approximation method (CRAM) was introduced as a novel method to solve the burnup equations. It can be characterized as the best rational function on the negative real axis and it has been shown to be capable of simultaneously solving an entire burnup system both accurately and efficiently. The main difficulty in using CRAM for computing the matrix exponential is determining the coefficients of the rational function for a given approximation order. Some polynomial CRAM coefficients have been published in 1984, and based on these literature values, CRAM approximations up to the order 16 have been thus far applied in burnup calculations. The topic of this paper is the computation of CRAM approximations and their application to burnup equations. A Remez-type method utilizing the equioscillation property of best approximations is used to construct the CRAM approximants for approximation orders 1, ..., 50. Numerical results are presented for a large burnup system and for a decay system. It is demonstrated that higher-order CRAM can be used to accurately solve the burnup equations even with time steps of the order of millions of years.