Generalized and Stability Rational Functions for Dynamic Systems of Reactor Kinetics

The base of reactor kinetics dynamic systems is a set of coupled stiff ordinary differential equations known as the point reactor kinetics equations.These equations which express the time dependence of the neutron density and the decay of the delayed neutron precursors within a reactor are first order nonlinear and essentially describe the change in neutron density within the reactor due to a change in reactivity. Outstanding the particular structure of the point kinetic matrix, a semianalytical inversion is performed and generalized for each elementary step resulting eventually in substantial time saving. Also, the factorization techniques based on using temporarily the complex plane with the analytical inversion is applied. The theory is of general validity and involves no approximations. In addition, the stability of rational function approximations is discussed and applied to the solution of the point kinetics equations of nuclear reactor with different types of reactivity. From the results of various benchmark tests with different types of reactivity insertions, the developed generalized Padé approximation (GPA) method shows high accuracy, high efficiency, and stable character of the solution.


Introduction
The multigroup neutron diffusion equations, which are derived from the general mathematical form of neutron transport equation, are represented generally by a system of time-and space-dependent coupled partial differential equations for which approximate solutions are sought from neutron transport equation.The point reactor kinetics equations are a useful simplification of the time-dependant neutronic neutron diffusion equations for most reactors operating near delayed critical.Included in that system are equations which describe the neutron level, reactivity, an arbitrary number of delayed neutron groups, and any other variables that enter into the reactivity equation in the case of reactivity feedback.
One of the important properties in a nuclear reactor is the reactivity, due to the fact that it is directly related to the control of the reactor.The startup process of a nuclear reactor requires that reactivity is varied in the system by lifting the control rods discontinuously.In practice, the control rods are withdrawn at time intervals such that reactivity is introduced in the reactor core linearly, to allow criticality to be reached in a slow and safe manner.For safety analysis and transient behaviour of the reactor dynamics, the neutron density and the delayed neutron precursor concentration are important parameters to be studied.Under simplified conditions, an analytical solution to the point reactor kinetic equations can be obtained, for example, such as there being no extraneous neutron source, assuming a prompt jump approximation [1][2][3][4][5][6], but they are elusive when the coefficients vary with time.Moreover, a time-dependent reactivity inserted into a point reactor is coupled multiplicity with neutron density to form a set of linear equations with time-dependent coefficients.This time dependence makes it difficult to obtain analytical solution, consuming CPU time on Pc's, and numerical integration is usually employed.Furthermore, the number of computations needed for each additional term in the explicit method for the exponential function makes the computing time prohibitive.For this and to avoid the instabilities associated with the computational effort involved in using the explicit methods, a particular class of approximations for the exponential function, namely, rational function approximations, is required.On the other International Journal of Nuclear Energy hand as is well known, a long-standing problem in point reactor kinetics is the stiffness arising from the orders of magnitude difference between prompt and delayed neutron lifetimes, which results in the requirement for very small time step increments in numerical solutions to the kinetics equations.
There has been a great deal of research that focused on eliminating the problem of stiffness.Several methods especially adapted for solving the initial value problems for stiff systems of ordinary differential equations.Among the methods are numerical integration using Simpson's rule [7], finite element method [8], Runge-Kutta procedures [9,10], piecewise polynomial approach [11,12], singular perturbation method [13][14][15], finite difference method [16], and other methods which are based on integral equation formulations with the slowly varying factor in each integrand represented by an assumed functional form [17][18][19][20].Most of these methods are successful in some specific problems but still suffer, more or less, from disadvantages, for example, limitations on the maximum permissible time step for ensuring computational stability.For slow transients in fast systems, even the analytic continuation method [21] should be changed in its basic procedure to achieve efficient computation, and inability to handle point reactor kinetics equations in their full generality, and so forth.
For quantitative computations, a numerical calculation is used.Furthermore, a very small step size has to be adopted for stability due to the problem of stiffness in the differential equations when applying a general numerical method.Many numerical methods have been proposed in the literature, such as stiffness confinement method [22] and power series solution [23][24][25][26].Several other creative schemes have been implemented as well as da Nóbrega [27], Lawrence and Doring [28], and Aboanber and Nahla [29] use a Padé approximation to the solution of the point kinetics equations.Also, the Padé approximation technique is applied to the nonlinear system including temperature feedback via analytical inversion method [30], CORE [31].Recently, a very efficient method for obtaining a highly accurate solution of reactor kinetic equation was proposed in [32][33][34].The basic idea of this approach is to combine classical discretization schemes with extrapolation techniques to reduce numerical error (e.g., discretization error and truncation error).
The present work is organized as follows.In Section 2 a newly developed general expression for a special type of functions, which allows us to approximate the exponential function in an economical manner, is derived.In addition, the factorization methods based on using temporarily the complex plane with the analytical inversion are generalized and applied to obtain the general solution of the point kinetics equations.In order to ensure the validity and stability of this class of approximations, the () test is applied to the rational function approximations (Padé approximations).Numerical examples of applying the method to a variety of primary problems confirm that the time step increment size can be greatly increased and that much computing time (CPU) can be saved, as compared to other conventional methods, in Section 3. Concluding remarks are summarized in Section 4.

Mathematical Background and Point Kinetics Model
2.1.Basic Idea.In most situations encountered in the analysis of different types of reactors, it is sufficient to model the neutronic behaviour of the reactor by a low order approximation to the formally exact neutron transport equation.Practically, analytical solutions of the coupled time-dependent transport and precursor equations for reactor kinetics problems are stiff and difficult; so approximate methods are usually used.The most widely used of these approximations is the multigroup neutron diffusion theory.In many operations and accident transients, large and fast changes in neutron distribution can occur.Therefore, a large number of recalculations of the shape function are needed to obtain an acceptable accuracy of calculations.The basic problem which it attacks is that of finding solutions to the space-time neutron group diffusion equation derived from the neutron Boltzmann transport equations for a reactor.Every approximation and stability study presented in this work is based on the point reactor kinetics equations, whose derivation and particular analytical solutions may be found in the standard reactor physics literature.In particular, Duderstadt and Hamilton [35] provide a rather intuitive approach to the equations while Henry [36] shows how they can be rigorously obtained from the timedependant diffusion equation.The point reactor kinetics equations are a useful simplification of the time-dependant neutronic Boltzmann problem for most reactors operating near delayed critical.The variables and parameters that appear in these equations are the result of integrating the original partial derivatives formulation of the neutron conservation problem with respect to both the spatial and angular coordinates and to all available neutron energies.Nevertheless, these particular details are not of interest for the analysis that follows.The "global" dynamic behavior of the reactor can be described adequately using the point reactor model.For the G-delayed neutron precursor group, the point reactor model is written as follows: with where () is the absolute reactivity defined as = ( − 1)/ =  ex /, which is generally a function of time, () is the reactor power (or neutron density),  is the effective multiplication factor, Λ() is the neutron generation time,   () is the precursor concentrations of a group  and   is the average decay constant,   is the effective delayed neutron fraction of a group, and () is the source.
Equations (1) may be compactly written in vector form.If denotes a  + 1 component column vector and denotes a ( + 1) × ( + 1) matrix, then (1) are simply A wealth of more or less economical methods have been devised to solve the system of (5).If the matrix elements [  ], in particular the reactivity, are constants, then exact solution of (5) can be easily obtained, for example, Porsching [19]; Duderstadt and Hamilton [35] and Henry [36] provide a rather intuitive approach to the equations.However, if as is usually the case () (and therefore ()) is a function of time, then (5) no longer has an exact solution.The main problem that has been recognized for a long time consists in the so-called stiffness of the system of (5).Stiff differential equations frequently arise in physical situation characterized by the existence of greatly differing time constants and are common in particular to such other fields as common engineering applications, circuit analysis, missile guidance system, chemical engineering, and chemical analysis.For a review of the numerical implications of this problem, see, for example, Gear's book [37].
Multiplying both sides of (5) by the integrating factor  −A() , we obtained the exact solution where the solution within each time step is sought starting from a piecewise constant approximation substituted in the interval [  ,  +1 ], with a constant, namely,   , which can be its value at some point in the interval (e.g.,   ≈ (  )    ≈ 1/2((  ) + ( +1 ))) as Integrating both sides of the equation over the th time step from   to  +1 gives Applying the usual constant parameters assumption, where () and () are replaced by their mean values  and , on the interval under consideration with ℎ =  +1 −   .Equation ( 8) reduces after subtracting and evaluating the integral to the equation of the form The numerical scheme proposed in the present investigation is based on this equation.One class of numerical methods for solving the kinetics equations is based upon the Taylor series approximation to the matrix series of ( 9).The number of computations needed for each additional term makes the computing time prohibitive.For this and to avoid the instabilities associated with the computational effort involved in using the explicit methods, (9) requires a particular class of approximations for the exponential function, namely, Padé rational approximations.

Generalities of Analytical Inversion Method.
Let us first state in this section the following important theorem before proceeding with the generalized rational function.The stated theorem for the matrices A and A  , Abaonber and Nahla [30], provides us with a mathematical formula, which permits to approximate some special functions.Let the matrix A have eigenvalues   and eigenvectors   ( = 1, 2, . . ., ℓ).Also, let A  have corresponding eigenvalues   , similar to its eigenvalues   which are the same as those of A. Thus, the following theorem is valid for any () for which (  ) is bounded for all .If A  =     and A    =     with      = 1 that is normalized to unity, so the following expression holds for any F(A) where F(A) represents a special type of functions (e.g., exponential functions, Trigonometric functions, and Hyperbolic functions) where the eigenvalues of both the matrices A and A  form a biorthonormal set when properly normalized.Equation ( 10) is so far a mere mathematical manipulation.It has, however, a form that permits us to approximate some special functions (e.g., the exponential function, exp(ℎ)) in an economical manner.The important note here is that if (  ) is a good approximation for F(  ), the th term in the summation should drop and consequently it will have a very small coefficient, Since exp(ℎ  ) ≈ (ℎ  ), thus, to a high degree of accuracy, we have where the sum ∑ /  is over only those  for which (11) does not hold.
The vectors U  and V  are easily calculated from their defining equation: or in the normalized form: where ]  is the normalization factor, satisfying the normalization condition      = 1, given by In this connection the analytical inversion method is applied to the point kinetic matrix A. For numerical purpose, exp(ℎ) is then approximated by some simple function of (ℎ), typically a real rational function of type (, ): with some normalization (usually  0 or  0 = 1).Whatever the degree, , of the denominator, it can always be factored into real linear and quadratic factors: For the sake of simplicity, suppose that  = 0, and Ψ +1 is approximated by which can be obtained by solving recursively the matrix problems: with Due to the particular structure of , a semianalytical inversion (Analytical Inversion Methods (AIMs) [30]) can be performed for each elementary step resulting eventually in substantial time savings.They developed a general expression for the inverse of [I−A], where for a real  and a point kinetic matrix A() of ( 4) with the reactivity () appropriately averaged, the following expression is introduced: where, generally for  delayed neutron group: Here, A = [a  ] is any point kinetics-like matrix, that is,   = 0 except when  = 1 or  = 1 or  = .Consequently, the previous expression is more general compared with that introduced by Aboanber and Nahla [30].The  ℓ in ( 18) being (generally) real numbers, the quadratic factors   will have complex conjugate roots.The expression ( 21) is of no great advantage by itself since we can solve directly the system of equations implied by the inverse shown with the same computational effort spending essentially in this case.However, the utility of the analytical inversion is evident when   form a complex conjugate pair.In this connection the following pair of factors is considered, with  =   +   , and  =   −   which is a real matrix and thus has a real inverse.This expression can be expressed in the general form after carrying out some involved but straightforward algebra: where with When A is the point kinetics matrix, this general expression reduces, of course, to the formula previously derived by Aboanber and Nahla [30] for the particular case of the Padé 02 approximations where we have   =   = ℎ/2.The remarkable feature of expressions ( 21) and (24) for the point kinetics matrix is that only the coefficient  depends on the reactivity; unless the rational approximation chosen changes from one time step to the other, it is thus possible to predetermine one for all the vectors a and b and the matrix D, while only modifying  at each time step, if at least the reactivity various with time.In general, the time dependence of the macroscopic cross-sections would also affect the   =   /Λ entries of the point kinetics matrix, so that these terms should generally be regarded as time-dependent quantities; in this case, one should also reevaluate vector a, but this should probably be not done at each time step, as the   is likely to change much slower than .

Generalized Padé Approximations.
Numerical methods for differential equations when applied to a linear test problem generate a sequence of approximations which satisfy a difference equation.In the problem of finding an approximation to the solution of such differential equations by numerical methods, many investigators show that rational approximation to the exponential function is of great importance.Among the rational approximations to the exponential, the Padé approximations are the most common.The Padé approximation often gives better approximation of the function than truncating its Taylor series, and it may still work where the Taylor series does not converge.For these reasons Padé approximations are used extensively in computer calculations.Some famous one-step integration methods have since long been recognized as particular Padé approximations, like the forward Euler (FE), backward Euler (BE), and Crank-Nicholson (CN) schemes that are equivalent to the Padé 10, Padé 01 and Padé 11 approximations, respectively.Many classical integration schemes of the Runge-Kutta type have been shown, moreover, to lead to the Padé approximation when applied to the equation Ψ/ = Ψ.
Given a rational function   (), the Padé ( + ) approximation to exp(ℎ) is obtained if the remaining ( +  + 1) coefficients are chosen so that the Maclaurin expansion of   () agrees with the Maclaurin expansion of exp() to ( +  + 1) terms.As such, the Padé approximations of exp() appear as local approximations of exp() at  = 0 and are consequently rather poor approximations for large  leading to restrictions on the time steps that can be taken.
In the present work, we follow the same ideas, but with different techniques; in particular, we derive and apply to the point kinetics equations global approximations of Padé in an attempt to increase the time steps to be used or, equivalently, to minimize the amount of work needed to achieve certain accuracy.For numerical purposes, exp(ℎ) is then approximated by some simple functions of (ℎ), typically a real rational function of type (, ) given by ( 16), with usual renormalization ( 0 or  0 = 1).
In the case of exponential function of ( 16), the explicit expressions for the   's and  ℓ 's are obtained.Dealing with the following representative equations for the exponential functions: the following system of algebraic equations of the unknown coefficients   's and  ℓ 's are obtained [29]: The utility of the previous system of algebraic equations is evident, especially, when one tries to determine any type of Padé approximations of any order.Actually, the Padé (, ) approximation to exp() is a one-point Padé approximation this point being  = 0, and we have  ()  (0) = exp () (0),  = 0, 1, . . .,  + .

Stability Considerations.
The stability property is strongly problem dependent.In order to get rid of this constraint, we pay more attention to convergence and stability problems, by introducing a test problem Ψ = Ψ and Ψ(0) = Ψ 0 , where  is a complex constant.Convergence guarantees us that we can approach the exact solution as closely as we want by diminishing correspondingly the time step, and the order of convergence gives us an idea of how fast this can be achieved when ℎ is reduced; these notions, especially the last, turn out to be of little practical interest when we are interested in approximating correctly the exact solution in the most economical way, that is, in a minimum number of time steps.This is true also when stability is considered, following the concept of time A-stability, a method is said to be A-stable if all numerical approximations tends to zero as the number of time steps tends to infinity when it is applied to the test differential equation Ψ = Ψ with a fixed positive ℎ and any complex constant , with Re() < 0 [37].This concept of stability is clearly much stronger than the usual (traditional) ones and must usually replace them when dealing with stiff systems of ODEs.
Applying the rational approximation methods developed in the previous section for different type of Padé approximations to the differential equation Ψ = Ψ, we obtain where ℎ = Δ and A-stability requires        (ℎ)      ≺ 1, ∀ with Re() ≺ 0, ∀ℎ ≻ 0, where   (ℎ) is any of the rational (, ) approximations considered.If Re() ≻ 0, the exact solution itself is growing exponentially and accuracy becomes the only requirement.A given integration method, (29), will be A-stable if it satisfies (30).This, in turn, will be true if and only if        (ℎ)      ≤ 1 on Re() = 0,   (ℎ) analytic in Re() ≺ 0, as a consequence of the maximum modulus theorem.Varga [38] has shown that all the Padé approximations with  ≥  are absolutely stable on the real negative axis:        (ℎ)      ≤ 1, ∀,  ≥ , ℎ real negative.
These approximations are thus possible if one is sure that the spectrum of A is real.Truly, this is not the case for the spacetime kinetics described in our case.A method is A()-stable if its stability domain is a wedge in the left half of the complex plane where the angle between the edge and the real axis is equal to .Consequently, it is preferable to use an A()-stable scheme for which a rational function   (ℎ) is such that where then the method is A()-stable (or acceptable), where  ∈ (0, /2).The angle  is only one of a number of parameters which have been proposed for measuring the extent of the stability region, but it is probably the best such measure.
As will be shown in the next section, many approximations which are proved to be accurate in most applications are not A-stable.

Numerical Results for Selected Benchmark Test Cases
In order to verify the developed method proposed in this study for the prediction of the neutron density distribution, the Padé approximation method was used for the numerical solution of the point kinetics equations.The nuclear parameters used in the validation of formulation proposed in this study which are to be used for the proposed benchmarks against the formulation proposed by [32][33][34] and other reference literatures are reported in Table 2. Several   authors claim to present benchmark solutions where they only report 2 to 3 accurate digits.Several reactor configurations are considered, including thermal and fast reactors and multiprecursor materials with varying delayed emission fractions.All benchmark problems correspond to transient cases starting from a critical state, with equilibrium between neutron and precursor concentration.The initial conditions  ( = 0) and   ( = 0) should be given to complete the formulation of the problem.In the absence of an external neutron source, the equilibrium condition for a steady-state density  0 is given by   (0) = (  /  Λ) 0 .In the following, the accuracy of the developed Padé method is benchmarked against the extremely accurate EPCA (Ganapol, 2013) and CATS methods [40], which can produce 9-digit or more reference values.The comparison with other results are also reported to highlight the performance of developed Padé technique.

Benchmark Series 1:
Step Reactivity Insertion.In order to compare the different methods a series of benchmarks is considered.The first example simulates a step-reactivity insertion without feedback.The Padé approximation method results are compared to well-established benchmarks.The first benchmark considered is for a step reactivity insertion of $−0.5 and $1 into a thermal reactor.The neutron generation time is 5 × 10 −4 s and the kinetics parameters are given in Table 2. Table 3 presents the results corresponding to this case for both negative and positive reactivity insertion in a thermal reactor configuration, thermal I.The average step size used in the developed FORTRAN code based on the present generalized Padé method (GPA) is 0.1 s.The converged densities agree to all of the digits of the recent reference solution EPCA and CATS.The approximately zero relative percentage error between GPA, EPCA, and GPA indicates that the relative error of the neutron density is much lower than the predefined error tolerance with the advantage of small time step size and of course CPU time.This simple preliminary test shows that the presented methods are consistent with the both EPCA and CATS solutions, when rounded, as well as with the analytical solution to the 9th digit.
Benchmark Series 2: Ramp Reactivity Insertion.The second type of reactivity triggering a transient is a ramp insertion.The transient is studied here in both a thermal (thermal II) and a fast reactor (fast I) configurations, for two values of the ramp.The reactivity driving function is  = 0.0007 for thermal reactor and  = 0.0044 for fast reactor.The neutron density at five time points between 0 and 10 s is computed and shown in Table 4.The reference value was provided by Hermite polynomial method with the fixed step size of 0.0001 s, Li et al. [20].An average step size of 0.001 s was obtained with the developed (GPA) method.By examining the relative error given in Table 4, the developed method gives accurate values consistently.All relative errors are much lower than the predefined error tolerance.-6.839867192 (07) A fast ramp of $1/s produces the neutron density given in Table 4.The stiffness is relatively strong.The neutron density and delayed neutron precursor density exponentially increase after 0.5 sec.The relatively small average step size of 0.001 sec of this case also shows that the step size can be automatically determined based on the stiffness of the problem.The relative error that is lower than the predefined tolerance in Table 4 also proves that the GPA method gives the exact results as the mathematical expectation.
For the ramp variation of reactivity as mentioned by Ganapol et al., 2013, even in the case of linear kinetics the PCA is inaccurate due to the variation of the reactivity within the discretization step, the problem which is not included in the GPA method.Indeed the numerical error due to the constant approximation can be reduced by choosing a time step sufficiently short as compared with reactivity variation and neutron lifetime.In the EPCA, this problem is reduced through coupling with the submesh that ensures an update of the reactivity within each time step.Results in Table 4 compare PCA and EPCA as well as CATS with the GPA benchmark solutions.While PCA gets the first four to at most seven digits, GPA is exact to at worst the 5th digit, which agrees very well with EPCA.Also the reference values reported in Table 4, showing that its accuracy is at best 4 digits of the true benchmark result.In addition, the GPA average step size of 0.01 sec is smaller than the GRK average step size of 0.02 sec [10].By comparing the GPA CPU time with PWS CPU time of 3.02 sec [24], the GPA method is about a hundred times faster than the PWS method.
Benchmark Series 3: Zigzag Reactivity Ramp.A zigzag ramp reactivity is inserted into a thermal reactor taken from [27], where the GPA method is applied to compute the response of a reactor core to a zigzag ramp input of reactivity.Consider The physical parameters in this example described in Table 2, Thermal Reactor I.The results obtained with the adopted GPA method are compared against those obtained by PCA and EPCA methods as well as with Chao and Attard [22] method (SCM).These computations are done with a fixed time step size of 0.01 s.According to all zero relative errors in Table 5, the GPA method provides accurate results.The errors associated with the four methods show the consistency of the GPA method.The neutron density profiles versus time are shown in Figure 1.The zigzag reactivity ramp brings fluctuation to the neutron density profile whose curve shape follows the reactivity driving function, but the relative magnitude of  the response of the delayed neutron precursor density to the zigzag reactivity ramp is relatively small, Figure 2.
Benchmark Series 4: Sinusoidal Reactivity.For further investigation of the developed method, more calculations were made.A sinusoidal reactivity insertion transient is simulated in Figure 3 for a thermal configuration (Thermal III).In this case the reactivity has the following form: () =  0 +  sin(/180), where the notation has its conventional meaning.The equations of reactivity combined with the system of point kinetics will give a crude simulation of a slow selflimiting excursion, where  0 is the constant part of the excess reactivity, and sin(/180) is a given function characterizing the time dependence of the reactivity.The parameter  is a positive number that represents the magnitude of the variable part of the excess reactivity in dollars.This parameter is assumed to be sufficiently small compared to unity.That is a real assumption since, except in the accidental case, the excess reactivity is always less than one dollar.The wave form oscillations of the flux is such that a modified sinusoidal curve is obtained for different values of the  parameter in Figure 3 with reactivity oscillation versus time.
Benchmark Series 5: Feedback Reactivity.The presence of temperature feedback is useful in providing an estimate of the transient behaviour of a reactor power.The previous benchmark series of reactivity inputs discussed the neutron kinetics (i.e., the response of neutron density in a nuclear reactor to an external reactivity input) under the implicit assumption that the level of the neutron density does not affect the properties of the system that determine the neutron kinetics, most notably the reactivity.This is the situation when the neutron density is sufficiently small that the fission heat does not affect the temperature of the system (i.e., zero power).However, in an operating nuclear reactor the neutron density is large enough that any change in fission heating resulting from change in neutron density will produce changes in temperature, which in turn will produce changes in reactivity, or in reactivity feedback.It is assumed that the reactor has a negative temperature coefficient of reactivity  ( > 0) when a small step reactivity  0 ( 0 < ) is inserted.When we consider the temperature feedback, the real reactor reactivity is given by where  is the temperature increment of the reactor.After the reactivity  0 is inserted into the reactor, the power responds quickly and the adiabatic model can be used for the calculation of reactor temperature [41,42].Then the derivative of temperature to time can be given as follows: where   is the reciprocal of thermal capacity of reactor.Combining ( 36) and (37) results in Generally, where () represented the impressed reactivity variation (generally polynomial in ), and () is the "shutdown coefficient" of the reactor system.It is usually taken as a negative constant of magnitude ranging from 10 −13 cm 3 /sec for slow system to 10 −7 cm 3 /sec for fast metal systems.The results of comparison for this case are shown in Table 6.

Conclusions
The point reactor kinetics equations are a useful simplification of the time-dependant neutronic Boltzmann problem that allow studying the stability of a reactor core using dynamical systems theory.From a strictly mathematical point of view, the point kinetics of reactor dynamic systems are of great importance for understanding the time-dependent behaviour of the neutron density.In nuclear reactor the response to either a planned or unplanned change in the reactor conditions is of great importance to the safe and reliable operation of the reactor.It is therefore important to understand the response of the neutron density and how it relates to the speed of lifting control rods (reactivity variations).The generalized Padé rational approximation (GPA) is used to approximate the exponential function in the point kinetics equation solution.General expression for the coefficients is developed for any type of   approximations.In addition, the roots of the inhour formula are used as the eigenvalues of the point kinetics matrix.A semianalytical inversion and factorization methods are developed and generalized for each elementary step resulting eventually in substantial time saving.It has been shown that GPA can exhibit a significant computational advantage by reducing the CPU time over more conventional approximations, and that it can be applied to a wide range of important problems.By comparing the GPA CPU time with PWS CPU time of 3.02 sec, the GPA method is about a hundred times faster than the PWS method as mentioned in benchmark.If larger time steps are considered, as reported in the presented tables, the GPA would still be quite useful over fairly time steps, especially at the extreme values of the inserted reactivity.The methods presented here should, however, retain much of their interest when applied to all cases where, after a small step-like transient, the reactivity reaches new stationary values.The results confirm the theoretical analysis and indicate the range of applicability of the presented methods.Their great advantage is that all we have said earlier remains valid for full space-time kinetics.In addition, the numerical tests demonstrated that high accuracy could be achieved in strong step reactivity insertion.The low step size was required to compute the ramp reactivity to acquire the same accuracy.

Figure 1 :Figure 2 :
Figure 1: Variations of neutron density for zigzag reactivity at different initial values.

Figure 3 :
Figure 3: Variations of neutron density for oscillatory reactivity at different initial values.

Table 1 :
Coefficients and explicit Padé approximations with their stability for the exponential function.exp(ℎA) = ∑

Table 2 :
Delayed-neutron parameters of typical reactors used in the benchmark test cases.

Table 3 :
Neutron density of benchmark 1: step of $−0.5 and $1 into a thermal reactor.

Table 6 :
Neutron density of benchmark 5: compensated response to ramp function reactivity changes with prompt neutron generation time, Λ = 5 * 10 −4 sec, and shutdown coefficients are in the range of 10 −11 cm 3 /sec to 10 −13 cm 3 /sec.