1. Introduction

Stochastic point kinetics equations (SPKEs) are a system of coupled nonlinear stochastic differential equations (SDEs) [1] and are important in nuclear engineering. The SPKEs model a system of Itô SDEs, specifically, neutron population and delayed neutron precursors. The physical dynamical system has been established to be a population process, and techniques have been employed by Hayes and Allen to transform deterministic point kinetics equations into a system of SDEs [2]. The fractional diffusion model is normally applied to large variations of neutron cross sections, which preclude the use of classical neutron diffusion equations [3-6]. Various approximation methods have been developed to improve the control of processes in the nuclear reactor.

In recent developments dynamic, multiphysics phenomena face many challenges regarding accurate numerical schemes, resulting in severe computational requirements. An approach to reduce the severe computational requirements of standard low-order simulations is to employ higher order formulations. In the hierarchy of high-order methods, compact schemes represent an attractive choice for reducing dispersion and anisotropy errors.

Nowak et al. presented the numerical solutions of a fractional neutron point kinetics model for a nuclear reactor [7]. Numerical solutions of SPKEs by implementing stochastic piecewise continuous approximation (PCA) have been obtained by Hayes and Allen [2], which provided a very succinct idea about the randomness of neutron density and precursor concentrations. Saha Ray showed that Euler–Maruyama (EM) and Taylor 1.5 strong order numerical schemes are well-founded estimators compared with stochastic PCA [8]. Saha Ray and Patra applied the Grunwald–Letnikov definition of fractional derivative for solving SPKEs [9]. Nahla and Edress [10] showed the efficiency of analytical exponential model (AEM) for obtaining solutions of SPKEs.

In our present work, we demonstrate the fractional numerical method for obtaining the mean neutron population for various reactivities. In Sect. 2, we introduce the fractional stochastic nonlinear point reactor kinetics equations (SNPKE). The fractional Itô SDE for NPKEs is obtained using the central limit theorem. In Sect. 3, we discuss the higher order approximation method for Caputo derivative, and in Sect. 4, we discuss its application. In Sect. 5, we discuss the obtained numerical solutions for various reactivities.

2. Framework of fractional SNPKE

We now introduce an elementary definition before stating the central limit theorem.

Definition 2.1

Let ${\overline{X}}_M$ be the sample mean of $M$ independent samples $X_1,\mathrm{...},X_M$ of a random variable X such that

where $\mu$ is the mean and ${\sigma }^2$ is the variance for each independent and identically distributed real-valued random variables $X_j$.

Theorem 2.1 (Central Limit Theorem)

If each random variable $X_j$ has a finite second moment with $Var{X}_j={\sigma }^2$, the distribution of ${\overline{X}}_M$ converges to that of a Gaussian random variable with mean $\mu$ and variance $\frac{{\sigma }^2}{M}$. Thus, $X_M^{*}$ converges in distribution to $Z~N(0,{\sigma }^2)$.

The fractional Itô SDE for NPKEs with temperature feedback effects [11-15] obtained from the central limit theorem can be written as follows [16]:

where $\alpha$ is the fractional order of the derivative and $0<\alpha \le 1$, $\Psi (t)=\left(\begin{array}{c}N(t)\\ C_1(t)\\ C_2(t)\\ ⋮\\ C_m(t)\end{array}\right)$, $W(t)=\left(\begin{array}{c}{W}_0(t)\\ W_1(t)\\ W_2(t)\\ ⋮\\ W_m(t)\end{array}\right)$, and $Q=\left(\begin{array}{c}q\\ 0\\ 0\\ ⋮\\ 0\end{array}\right)$. Here, $m$ is the total number of delayed neutron groups and $W_0(t),W_1(t),W_2(t),\dots ,W_m(t)$ are standard Wiener processes as defined in [8], $N(t)$ is the neutron population, and $C_i(t)$ is the precursor concentration of $i$-group of delayed neutrons.

The coefficient matrix $A(t)$ is represented in the following form $=\left(\begin{array}{c}{W}_0(t)\\ W_1(t)\\ W_2(t)\\ ⋮\\ W_m(t)\end{array}\right)$

where $\rho$ is the total reactivity, $\beta ={\displaystyle \sum _{i=1}^m{\beta }_i}$ is the total fraction of delayed neutrons, ${\beta }_i$ is the fraction and ${\lambda }_i$ is the decay constant of $i$-group of delayed neutrons, and $l$ is the prompt neutron generation time.

The covariance matrix $B(t)$, which is evaluated in reference [16], is represented as follows:

where ${\mu }_0(t)=\left(\frac{\rho +\beta }{l}\right)N(t)-{\displaystyle \sum _{i=1}^m{\lambda }_i{C}_i(t)}$, ${\mu }_i(t)=\frac{{\beta }_i}{l}N(t)-{\lambda }_i{C}_i(t)$, and $i=1,2,3,\dots ,m$.

3. Higher order approximation scheme

3.1 Fractional calculus

The Caputo derivative operator for $\alpha \in (0,1)$ is defined as follows:

${}_0{}^{C}D{}_t^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}{\displaystyle \underset{0}{\overset{t}{\int }}{(t-s)}^{-\alpha }{f}^{\prime }(s)ds},$

in which $\Gamma (\cdot )$ is the Euler gamma function.

In this paper, a ${(3-\alpha )}^{th}$ order scheme for Caputo derivative ${}_0{}^{C}D{}_t^{\alpha }\Psi (t)$ with $\alpha \in (0,1)$ has been discussed for obtaining the solutions of fractional SNPKEs. Let $0=t_0<t_1<\mathrm{...}<t_M=T$ and $t_n=t_0+n\tau$ be the equidistant discretized times with $\tau ={\Delta }_n=\frac{(T-t_0)}{M}$ for $M\in Z$ such that $\tau \in (0,1)$. Now using the Taylor expansion to ${\Psi }^{\prime }(s)$, $\Psi (t_{i-1})$, and $\Psi (t_{i+1})$ at point $t=t_i$$(0\le i<n)$, we get

${\Psi }^{\prime }(s)={\Psi }^{\prime }(t_i)+{\Psi }^{″}(t_i)(s-t_i)+\frac{{\Psi }^{‴}(t_i)}{2!}{(s-t_i)}^2+O({(s-t_i)}^3)$, $s\in (t_i,t_{i+1}),$

${\Psi }^{\prime }(t_i)=\frac{\Psi (t_{i+1})-\Psi (t_{i-1})}{2\tau }-\frac{{\Psi }^{‴}(t_i)}{3!}{\tau }^2+O({\tau }^4)$, and

${\Psi }^{″}(t_i)=\frac{\Psi (t_{i+1})-2\Psi (t_i)+\Psi (t_{i-1})}{{\tau }^2}-\frac{{\Psi }^{(4)}(t_i)}{12}{\tau }^2+O({\tau }^4).$

Hence, we can obtain the following:

${\Psi }^{\prime }(s)=\frac{\Psi (t_{i+1})-\Psi (t_{i-1})}{2\tau }+\frac{\Psi (t_{i+1})-2\Psi (t_i)+\Psi (t_{i-1})}{{\tau }^2}(s-t_i)-\frac{{\Psi }^{‴}(t_i)}{3!}{\tau }^2+\frac{{\Psi }^{‴}(t_i)}{2!}{(s-t_i)}^2+O({(s-t_i)}^3)$, $0<s-t_i<\tau$.

Therefore, the Caputo derivative can be discretized as

$\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }\Psi (t){|}_{t=t_n}=\frac{1}{\Gamma (1-\alpha )}{\displaystyle \underset{0}{\overset{t_n}{\int }}{(t_n-s)}^{-\alpha }{\Psi }^{\prime }(s)ds}\\ =\frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}{\displaystyle \underset{t_i}{\overset{t_{i+1}}{\int }}{(t_n-s)}^{-\alpha }{\Psi }^{\prime }(s)ds}}\\ =\frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}{\displaystyle \underset{t_i}{\overset{t_{i+1}}{\int }}{(t_n-s)}^{-\alpha }[\frac{\Psi (t_{i+1})-\Psi (t_{i-1})}{2\tau }+\frac{\Psi (t_{i+1})-2\Psi (t_i)+\Psi (t_{i-1})}{{\tau }^2}(s-t_i)}}\\ -\frac{{\Psi }^{‴}(t_i)}{3!}{\tau }^2+\frac{{\Psi }^{‴}(t_i)}{2!}{(s-t_i)}^2]ds+O({\tau }^3)\\ =\frac{{\tau }^{-\alpha }}{\Gamma (3-\alpha )}{\displaystyle \sum _{i=0}^{n-1}[w_{1,n-i}({\Psi }_{i+1}-{\Psi }_{i-1})+w_{2,n-i}({\Psi }_{i+1}-2{\Psi }_i+{\Psi }_{i-1})]}+r^n.\end{array}$

where $n=1,2,\mathrm{...},M$ and ${\tau }^n$ is the truncation error.

Therefore, the Caputo derivative has the following numerical approximation [17]:

${}_0{}^{C}D{}_{t_n}^{\alpha }\Psi (t_n)=\frac{{\tau }^{-\alpha }}{\Gamma (3-\alpha )}{\displaystyle \sum _{i=0}^{n-1}[w_{1,n-i}({\Psi }_{i+1}-{\Psi }_{i-1})+w_{2,n-i}({\Psi }_{i+1}-2{\Psi }_i+{\Psi }_{i-1})]}+O({\tau }^{3-\alpha }),$ (3.1)

where $0<\alpha <1$, $w_{1,n-i}=\frac{2-\alpha }{2}[{(n-i)}^{1-\alpha }-{(n-i-1)}^{1-\alpha }]$, $w_{2,n-i}={(n-i)}^{2-\alpha }-{(n-i-1)}^{2-\alpha }-(2-\alpha ){(n-i-1)}^{1-\alpha }$, and $r^n$ is the truncation error in the following form:

$r^n=\frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}{\displaystyle \underset{t_i}{\overset{t_{i+1}}{\int }}{(t_n-s)}^{-\alpha }[-C_{\Psi }{\tau }^2+3C_{\Psi }{(s-t_i)}^2]ds+O({\tau }^3)}},$ (3.2)

where $C_{\Psi }=\frac{{\Psi }^{‴}(t_i)}{3!}$ is a constant.

The right-hand side of eq. (3.2) can be expressed as follows:

$\frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}{\displaystyle \underset{t_i}{\overset{t_{i+1}}{\int }}{(t_n-s)}^{-\alpha }[-C_{\Psi }{\tau }^2+3C_{\Psi }{(s-t_i)}^2]ds+O({\tau }^3)}}=\frac{C_{\Psi }}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}(I_1+.3I_2)},$

where

Therefore,

$\begin{array}{l}\frac{C_{\Psi }}{\Gamma (1-\alpha )}{\displaystyle \sum _{i=0}^{n-1}(I_1+.3I_2)}=\frac{C_{\Psi }{\tau }^{3-\alpha }}{\Gamma (2-\alpha )}\{-n^{1-\alpha }-3[{(n-1)}^{1-\alpha }+\mathrm{...}+2^{1-\alpha }+1]-\frac{6}{2-\alpha }[{(n-1)}^{3-\alpha }+\mathrm{...}+2^{3-\alpha }+1]\\ +\frac{6}{(2-\alpha )(3-\alpha )}{n}^{3-\alpha }\}.\end{array}$

Let

$\begin{array}{l}S(n)=-3{\displaystyle \sum _{i=1}^{n-1}{i}^{1-\alpha }}-\frac{6}{2-\alpha }{\displaystyle \sum _{i=1}^{n-1}{i}^{2-\alpha }}+\frac{6}{(2-\alpha )(3-\alpha )}{n}^{3-\alpha }-n^{1-\alpha }\\ ={\displaystyle \sum _{i=1}^{n-1}{a}_i},n\ge 1.\end{array}$

If $n=1$, define $a_0=s(1)=\frac{6}{(2-\alpha )(3-\alpha )}-1$. Then $a_i$($i\le 1$) can be defined as follows:

$a_i=S(i+1)-S(i)=\frac{6}{(2-\alpha )(3-\alpha )}[{(i+1)}^{3-\alpha }-i^{3-\alpha }]-\frac{6}{2-\alpha }{i}^{2-\alpha }-{(i+1)}^{1-\alpha }-2i^{1-\alpha }.$

It can be proven that $|S(n)|$ is bounded for $n\ge 1$ [17-19]. This proves that the series ${\displaystyle \sum _{i=0}^{\infty }{a}_i}$ converges. On further simplification we get,

$\begin{array}{l}{}_0{}^{C}D{}_{t_n}^{\alpha }\Psi (t_n)=\frac{{\tau }^{-\alpha }}{\Gamma (3-\alpha )}[{\Psi }_n(w_{1,1}+w_{2,1})+{\displaystyle \sum _{i=0}^{n-2}{w}_{1,n-i}{\Psi }_{i+1}}-{\displaystyle \sum _{i=0}^{n-1}{w}_{1,n-i}{\Psi }_{i-1}}+{\displaystyle \sum _{i=0}^{n-2}{w}_{2,n-i}{\Psi }_{i+1}}\\ +{\displaystyle \sum _{i=0}^{n-1}{w}_{2,n-i}(-2{\Psi }_i+{\Psi }_{i-1})}]+O({\tau }^{3-\alpha })\end{array}$ (3.3)

In Eq. (3.2), if $i=0$, then ${\Psi }_{i-1}={\Psi }_{-1}$, which lies outside of $[0,T]$. Various options have been provided to approach ${\Psi }_{-1}$. In numerical calculations, the neighboring function values have been used to approximate ${\Psi }_{-1}$, that is, ${\Psi }_{-1}=\Psi (0)-\tau {\Psi }^{\prime }(0)+\frac{{\tau }^2}{2}{\Psi }^{″}(0)+O({\tau }^3)$.

1) When ${\Psi }^{\prime }\text{(}0)={\Psi }^{″}\text{(}0)=0$, then ${\Psi }_{-1}={\Psi }_0+O({\tau }^3)$; the convergence order is $O({\tau }^{3-\alpha })$.

2) When ${\Psi }^{\prime }\text{(}0)=0$, ${\Psi }^{″}\text{(}0)\ne 0$, then ${\Psi }_{-1}={\Psi }_0+\frac{{\tau }^2}{2}{\Psi }^{″}(0)+O({\tau }^3)$; the convergence order of Eq. (3.1) is $O({\tau }^2)$.

3) When ${\Psi }^{\prime }\text{(}0)\ne 0$, then the convergence order is $O(\tau )$.

4. Solution of SPKE by higher order approximation method

In this section, higher order approximation to Caputo derivative has been applied to Eq. (2.1) as follow:

Therefore, the above expression can be simplified into an explicit numerical scheme as follows:

where $\Delta W(t_n)=\sqrt{h}{S}_n$ and $n=1,2,\mathrm{...},M$ with initial condition

where $S_{0_n},S_{1_n},S_{2_n},\dots ,S_{m_n}$ are random variables from $N(0,1)$

with mean equal to 0 and variance equal to 1.

Theorem 4.1

The local truncation error of the scheme is $O({\tau }^{3-\alpha })$.

Proof: The local truncation error of the higher order scheme for Eq. (4.2) can be derived as follows:

5. Numerical results and discussions

Here the solutions of SPKE ($i=6$) have been obtained using higher order approximation scheme.

5.1 Step reactivity

In this section, the numerical solutions of the fractional stochastic point kinetic model [10] using the following parameters have been obtained: ${\lambda }_1=0.0127{\text{s}}^{-1}$, ${\lambda }_2=0.0317{\text{s}}^{-1}$, ${\lambda }_3=0.115{\text{s}}^{-1}$, ${\lambda }_4=0.311{\text{s}}^{-1}$, ${\lambda }_5=1.4{\text{s}}^{-1}$, ${\lambda }_6=3.87{\text{s}}^{-1}$, ${\beta }_1=0.000266$, ${\beta }_2=0.001491$, ${\beta }_3=0.001316$, ${\beta }_4=0.002849$, ${\beta }_5=0.000896$, ${\beta }_6=0.000182$, $\beta =0.007$, $l=2.0\times {10}^{-5}\text{s}$, $v=2.5$, and $q=0$ with $N(0)=N_0=100$ (neutron) and $C_i(0)=\frac{{\beta }_{i}N(0)}{l{\lambda }_i}$.

For different step reactivities, i.e., $\rho =0.003$ and $\rho =0.007$, the mean peak values of $N(t)$ with respect to time for fractional orders $\alpha =0.96$, 0.98, and 0.99 at step size $h=0.001\text{s}$ and for 500 trials are presented in Table 1. In addition, the results have been compared with those obtained by other methods, namely EM [8], Taylor 1.5 strong order [8], AEM [20], and efficient stochastic model (ESM) [10] with graphical representation in Fig. 1. The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher order approximation method.

Table 1:Full-screen View

Mean peak values of $N(t)$ for different step reactivities.

ρ	α = 0.96		α = 0.98		α = 0.99		EM α = 1[8]	Taylor 1.5 strong order α = 1 [8]	AEM α = 1[20]	ESM α = 1[10]
	Peak	Time (s)	Peak	Time (s)	Peak	Time (s)	Peak	Peak	Peak	Peak
0.003	180.796	0.1	180.033	0.095	180.819	0.083	208.6	199.408	186.30	179.93
0.007	128.655	0.001	128.248	0.001	124.113	0.001	139.568	139.569	134.54	134.96

Show more

Fig. 1.Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) step reactivity $\rho =0.003$ and $\alpha =0.96$, (b) step reactivity $\rho =0.003$ and $\alpha =0.98$, (c) step reactivity $\rho =0.003$ and $\alpha =0.99$.

5.2 Ramp reactivity

The numerical solutions of the fractional SNPKE have been obtained using the same parameters as those in Sect. 5.1. Here, the reactivity can be represented as $\rho =0.1\beta t$.

For ramp external reactivity, i.e., $\rho =0.1\beta t$, the mean peak values of $N(t)$ with respect to time for fractional orders $\alpha =0.96$, 0.98, and 0.99 at step size $h=0.001(s)$ and for 500 trials are listed in Table 2. In addition, the results have been compared with those obtained using other methods, namely AEM [20] and ESM [10] with graphical representation in Fig. 2. The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher order approximation method.

Table 2:Full-screen View

Mean peak values of $N(t)$ for ramp reactivity $\rho =0.1\beta t$ and different values of fractional order $\alpha$.

α	σ	α = 0.96		α = 0.98		α = 0.99		AEM α = 1 [20]	ESM α = 1 [10]
		Peak	Time (s)	Peak	Peak	Peak	Time (s)	Peak	Peak
0.1βt	10-11	113.563	0.998	113.275	186.30	113.045	0.1	113.267707	113.116433

Show more

Fig. 2.Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) ramp reactivity $\rho =0.1\beta t$ and $\alpha =0.96$, (b) ramp reactivity $\rho =0.1\beta t$ and $\alpha =0.98$, (c) ramp reactivity $\rho =0.1\beta t$ and $\alpha =0.99$.

5.3 Temperature feedback reactivity

In this section, solutions of fractional stochastic point kinetic model with $i$-group of delayed neutrons ($i=6$) in the presence of Newtonian temperature feedback have been obtained using the higher order approximation method.

The total reactivity of the reactor in the presence of temperature feedback [21] is represented in the following form:

$\rho (t)={\rho }_{ex}(t)-{\rho }_f(t), {\rho }_f(t)=\sigma {\displaystyle \underset{0}{\overset{t}{\int }}N(\tau )d\tau }$

where $\sigma =\alpha K_c$, ${\rho }_{ex}(t)$ is the external reactivity, $T(t)$ is the temperature, $T_0$ is the initial temperature, $\alpha$ is the temperature coefficient, and $K_{\text{c}}$ is the reciprocal of the thermal capacity.

In the interval $[t_k,t_{k+1}]$, the total reactivity can be expressed as follows [15]:

$\rho (t_k)\approx {\rho }_{\text{ex}}(t_k)-h\sigma {\displaystyle \sum _{j=0}^{k}N(t_j)}.$

5.3.1 Step external reactivity

In this section, the numerical solutions of the SPKE of $U^{235}$ nuclear reactor [3] using the following parameters have been obtained: ${\lambda }_1=0.0124{\text{ s}}^{-1}$, ${\lambda }_2=0.0305{\text{ s}}^{-1}$, ${\lambda }_3=0.111{\text{ s}}^{-1}$, ${\lambda }_4=0.301{\text{ s}}^{-1}$, ${\lambda }_5=1.13{\text{ s}}^{-1}$, ${\lambda }_6=3.0{\text{ s}}^{-1}$, ${\beta }_1=0.00021$, ${\beta }_2=0.00141$, ${\beta }_3=0.00127$, ${\beta }_4=0.00255$, ${\beta }_5=0.00074$, ${\beta }_6=0.00027$, $\beta =0.00645$, $l=5.0\times {10}^{-5}\text{s}$, $\alpha =5.0\times {10}^{-5}{\text{K}}^{-1}$, and $K_{\text{c}}=0.05{\text{ K/MW}}_{\text{s}}$ with $N(0)=N_0=1$ (neutron) and $C_i(0)=\frac{{\beta }_{i}N(0)}{l{\lambda }_i}$.

For different step external reactivities, i.e., ${\rho }_{\text{ex}}=0.5\beta$, ${\rho }_{\text{ex}}=0.75\beta$, and ${\rho }_{\text{ex}}=\beta$, the mean peak values of $N(t)$ for fractional orders $\alpha =0.96$, 0.98, and 0.99 and for different step external reactivities $0.5\beta$, $0.75\beta$, and $\beta$ using 500 trials are presented in Table 3(a). These results have been compared with previously obtained results of split-step forward EM method (SSFEMM) and derivative-free Milstein method (DFMM) [15] in Table 3(b) with graphical representation in Figs. 3, 4, and 5. The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher order approximation scheme.

Table 3Full-screen View

Mean peak values of N(t) for ρ_ex=0.5β, ρ_ex=0.75β, and ρ_ex=β, b comparison of mean peak values of N(t) for ρ_ex=0.5β, ρ_ex=0.75β, and ρ_ex=β for α=1, and α=0.98

ρex	α = 0.96		α = 0.98		α = 0.99		ρex	SSFEMM [15] (α = 1)		DFMM [15] (α = 1)		α = 0.98
	Peak	Time (s)	Peak	Time (s)	Peak	Time (s)		Peak	Time (s)	Peak	Time (s)	Peak	Time (s)
(a)							(b)
0.5β	42.6182	28.65	44.789	30.29	45.9708	29.25	0.5β	46.4939	28.34	46.2606	27.84	44.789	30.29
0.75β	159.21	8.875	160.124	8.895	162.99	9.305	0.75β	163.707	8.795	164.22	8.95	160.124	8.895
β	801.166	0.985	795.268	1.03	772.893	1.0625	β	760.589	1.065	769.238	1.0575	795.268	1.03

Show more

Fig. 3Full-screen View

(Color online) Mean N(t) and two arbitrary sample paths for (a) step reactivity ${\rho }_{\text{ex}}=0.5\beta$ and $\alpha =0.96$, (b) step reactivity ${\rho }_{\text{ex}}=0.75\beta$ and $\alpha =0.96$, (c) step reactivity ${\rho }_{\text{ex}}=\beta$ and $\alpha =0.96$.

Fig. 4Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) step reactivity ${\rho }_{ex}=0.5\beta$ and $\alpha =0.98$, (b) step reactivity ${\rho }_{ex}=0.75\beta$ and $\alpha =0.98$, (c) step reactivity ${\rho }_{ex}=\beta$ and $\alpha =0.98$.

Fig. 5Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) step reactivity ${\rho }_{\text{ex}}=0.5\beta$ and $\alpha =0.99$, (b) step reactivity ${\rho }_{\text{ex}}=0.75\beta$ and $\alpha =0.99$, (c) step reactivity ${\rho }_{\text{ex}}=\beta$ and $\alpha =0.99$.

For $\alpha =0.96$, $0.98$, and $0.99$, the mean $N(t)$ with two arbitrary sample paths have been shown for ${\rho }_{\text{ex}}=0.5\beta$, ${\rho }_{\text{ex}}=0.75\beta$, and ${\rho }_{\text{ex}}=\beta$.

5.3.2 Ramp external reactivity

The numerical solutions of the fractional SNPKE of $U^{235}$ nuclear reactor have been obtained using parameters similar to those in subsection 5.3.1. Here, the external reactivity is represented as ${\rho }_{\text{ex}}=0.1t$ and ${\rho }_{\text{ex}}=0.01t$, the nonlinear coefficient $\sigma$ takes ${10}^{-11}$ or ${10}^{-13}$.

For different ramp external reactivities, i.e., ${\rho }_{ex}=0.1t$ and $0.01t$, the mean peak values of $N(t)$ with respect to time for fractional orders $\alpha =0.96$, 0.98, and 0.99 at step size $h=0.001(s)$ and for 500 trials are listed in Table 4(a). These results have been compared with previously obtained results [10, 15] in Table 4(b) with graphical representation in Figs. 6(a)–(c), 7(a)–(c), and 8(a)–(c). The solutions of SPKE obtained using the above-discussed fractional scheme have been tabulated to establish the efficiency of the higher order approximation method.

Table 4Full-screen View

Mean peak values of N(t) for ${\rho }_{ex}=0.1t$ and $0.01t$, $\rho (t)=at-\sigma {\displaystyle \underset{0}{\overset{t}{\int }}N(\tau )d\tau }$, b comparison of mean peak values of $N(t)$ for ${\rho }_{ex}=0.1t$ and $0.01t$, $\rho (t)=at-\sigma {\displaystyle \underset{0}{\overset{t}{\int }}N(\tau )d\tau }$ for $\alpha =1$ and $\alpha =0.98$

α	σ	α = 0.96		α = 0.98		α = 0.99		α	σ	SSFEMM [15] (α = 1)		DFMM [15] (α = 1)		AEM [10] (α = 1)		α = 0.98
		Peak	Time (s)	Peak	Time (s)	Peak	Time (s)			Peak	Time (s)	Peak	Time (s)	Peak	Time (s)	Peak	Time (s)
(a)								(b)
0.01	10-11	1.69869E + 10	1.084	1.76031E + 10	1.103	1.73211E + 10	1.112	0.01	10-11	1.68604E + 10	1.118	1.69492E + 10	1.118	1.673436E + 10	0.854	1.76031E + 10	1.103
	10-13	2.0663E + 12	1.132	2.18262E + 12	1.151	2.13422E + 12	1.161		10-13	2.12034E + 12	1.169	2.12802E + 12	1.168	2.082531E + 12	0.877	2.18262E + 12	1.151
0.1	10-11	1.78236E + 11	0.223	1.90873E + 11	0.232	2.01965E + 11	0.235	0.1	10-11	1.88849E + 11	0.235	1.89642E + 11	0.235	1.790577E + 11	0.142	1.90873E + 11	0.232
	10-13	2.34211E + 13	0.238	2.37627E + 13	0.248	2.41795E + 13	0.252		10-13	2.24448E + 13	0.251	2.26026E + 13	0.251	2.143778E + 13	0.150	2.37627E + 13	0.248

Show more

Fig. 6Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-11}$, and $\alpha =0.96$, (b) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-13}$, and $\alpha =0.96$, (c) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-11}$, and $\alpha =0.96$, (d) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-13}$, and $\alpha =0.96$.

Fig. 7Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-11}$, and $\alpha =0.98$, (b) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-13}$, and $\alpha =0.98$, (c) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-11}$, and $\alpha =0.98$, (d) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-13}$, and $\alpha =0.98$.

Fig. 8Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for (a) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-11}$, and $\alpha =0.99$, (b) ramp reactivity ${\rho }_{ex}=0.01t$, $\sigma ={10}^{-13}$, and $\alpha =0.99$, (c) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-11}$, and $\alpha =0.99$, (d) ramp reactivity ${\rho }_{ex}=0.1t$, $\sigma ={10}^{-13}$, and $\alpha =0.99$.

5.4 Sinusoidal reactivity

In this section, the reactivity is in the form of sinusoidal change, i.e., $\rho ={\rho }_0\sin \left(\frac{\pi t}{T}\right)$. The numerical solution for this reactivity has been obtained using the following parameters: ${\rho }_0=0.005333$, ${\beta }_1=\beta =0.0079$, ${\lambda }_1=0.077$, $\Lambda =0.001$, $q=0$, $N(0)=N_0=1$, and time period $T=100(s)$. The mean peak values of $N(t)$ with respect to time for fractional orders $\alpha =0.96$, 0.98, and 0.99 are presented in Table 5 with graphical representation in Figs. 9.

Table 9:Full-screen View

Mean peak values of $N(t)$ for sinusoidal reactivity $\rho =0.005333\sin \left(\frac{\pi t}{T}\right)$ for different values of fractional order $\alpha$.

ρ	α = 0.96		α = 0.98		α = 0.99
	Peak	Time (s)	Peak	Time (s)	Peak	Time (s)
$0.005333\sin \left(\frac{\pi t}{T}\right)$	38.0005	38.28	45.8029	38.18	49.1345	38.99

Show more

Fig. 9Full-screen View

(Color online) Mean $N(t)$ and two arbitrary sample paths for sinusoidal reactivity $\rho =0.005333\sin \left(\frac{\pi t}{T}\right)$ with (a) fractional order $\alpha =0.96$, (b) fractional order $\alpha =0.98$, and (c) fractional order $\alpha =0.99$.

6. Conclusion

In this study, fractional SPKEs have been solved using the higher order approximation scheme with different fractional orders $\alpha$. The obtained numerical solutions for mean $N(t)$ have been presented in the tables and graphically demonstrated to justify the efficiency of the proposed higher order approximation method. In addition, the results have been compared to some previous works such as [10, 15, 17]. The results obtained by the implemented fractional model is in good agreement with previous results, which further establishes the efficiency of our proposed scheme. The graphical representation for different reactivities shows the behavior of the mean neutron population. The random fluctuations at low power levels and achievement of equilibrium state after reaching its peak value provide us with a succinct idea about the behavior of $N(t)$ for different reactivities.