Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics

Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics 2 corresp first-author ahmed.aboanber@science.tanta.edu.eg ahmed.aboanber@science.tanta.edu.eg ahmed.aboanber@science.tanta.edu.eg 1 2 corresp a.nahla@science.tanta.edu.eg a.nahla@science.tanta.edu.eg a.nahla@science.tanta.edu.eg 2 corresp adel.edress@science.tanta.edu.eg adel.edress@science.tanta.edu.eg adel.edress@science.tanta.edu.eg Department of Mathematics, Faculty of Science, Taibah University, Al-Madinah 41411, Saudi Arabia Department of Mathematics, Faculty of Science, Taibah University, Al-Madinah 41411, Saudi Arabia Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt Fractional stochastic kinetics equations have proven to be valuable tools for the point reactor kinetics model, where the nuclear reactions are not fully described by deterministic relations. A fractional stochastic model for the point kinetics system with multi-group of precursors, including the effect of temperature feedback, has been developed and analyzed. A major mathematical and inflexible scheme to the point kinetics model is obtained by merging the fractional and stochastic technique. A novel split-step method including mathematical tools of the Laplace transforms, Mittage-Leffler function, eigenvalues of the coefficient matrix and its corresponding eigenvectors have been used for the fractional stochastic matrix differential equation. The validity of the proposed technique has been demonstrated via calculations of the mean and standard deviation of neutrons and precursor populations for various reactivity: step, ramp, sinusoidal, and temperature reactivity feedback. The results of the proposed method agree well with the conventional one of the deterministic point kinetics equations. Itĉ stochastic point kinetics equations Temperature feedback effects Wiener process Fractional calculus Mittage-Leffler function 1 Introduction 1 The dynamical processes described by the linear or nonlinear point kinetics equations are random processes in nature, that is due to the neutron population and precursor concentrations of delayed neutrons varying randomly with time. At the levels of high power, the stochastic manner is imperceptible, while at low-power levels, for example, at the start up of reactor operation, random fluctuation in the neutron population density and neutron precursor concentrations can be useful. The behavior variations of neutron population and precursor concentrations for nuclear reactors have been described by several innovators through employing stochastic models. Hayes and Allen [ 1 1 ] are the first authors who derived the stochastic model of the linear point reactor kinetics equations. They introduced a simplified stochastic model based on the Itô stochastic differential equations. The numerical results of this model using stochastic piecewise constant approximation (SPCA) have been compared with the Monte Carlo (MC) calculations and the experimental measurements [ 2 2 ]. Ha and Kim [ 3 3 , 4 4 ] have presented the stochastic space-dependent kinetics model (SSKM) to solve the one dimensional monoenergetic space-time reactor kinetics. In 2012, Ray [ 5 5 ] developed Taylor 1.5 strong order methods and the Euler–Maruyama to solve the stochastic point kinetics equations with step reactivity, while in 2013, Ray and Patra [ 6 6 ] presented the same techniques with sinusoidal reactivity. The stochastic partial differential equation and stochastic difference equations have been presented by Ref. [ 7 7 ] for neutron transport equation. Furthermore, the power doubling time for a subcritical reactor is identified through the point kinetics system by Allen as a stochastic first-passage time problem [ 8 8 ]. In 2013, Ray and Patra [ 9 9 ] were the first authors who presented a numerical solution of f ractional stochastic neutron point kinetic equations. Ayyoubzadeh and Vosoughi [ 10 10 ] simplified the system of Itô stochastic differential equations via alternative derivation of the stochastic differential equations. In 2016, Nahla and Edress [ 11 11 ] utilized the analytical exponential model (AEM) for the simplest formula of the stochastic point reactor kinetics system with various reactivity. They also proposed an efficient stochastic model (ESM) for the point kinetics model in Ref. [ 12 12 ]. da Silva et al. [ 13 13 ] presented a solution for the stochastic neutron point kinetics model. In 2017, Nahla [ 14 14 ] developed the analytical exponential technique (AET) to solve a stochastic nonlinear system of the point reactor kinetics equations with Newtonian temperature feedback reactivity. Finally, Singh and Ray [ 15 15 ] presented a comparison of two split-step methods for the numerical simulation of stochastic point kinetics equations in the presence of Newtonian temperature feedback effects. The fundamental objective of this work is speculation of the stochastic point kinetics system to a fractional stochastic point kinetics system including multi-group of the precursor concentration. To overcome the difficulty arising from the merging of the fractional and stochastic techniques, a developed mathematical technique is presented for solving the equation in the matrix form of the proposed fractional stochastic model. This technique is based on split-step technique, Laplace transforms, the Mittage-Leffler function, eigenvalues of the coefficient matrix and its corresponding eigenvectors. The proposed method is applied to the fractional stochastic point kinetics system with various reactivity and different fractional order. The paper is sorted out as follows: the preliminaries of the stochastic model, as well as the definitions of the Wiener process and fractional calculus, are introduced in Sect. 2. The solution of the fractional stochastic point kinetics equations with multi-group of a delayed precursor is derived in Sect. 3. The computational numerical results of the proposed system are discussed and compared with various stochastic techniques in Sect. 4. General conclusions including future work are presented in Sect. 5. 2 Preliminaries 1 In the following subsections, preliminaries of the Itô stochastic model, Wiener process, and fractional calculus are introduced briefly. 2.1 Itô stochastic model 2 Let us consider the following Itô stochastic differential model: [ 16 16 ] <math display="block" id="M1"><mrow><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>g</mi><mo stretchy="false">(</mo><mi>t</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext>d</mtext><mi>t</mi><mo>+</mo><mi>f</mi><mo stretchy="false">(</mo><mi>t</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> where, <math id="M2"><mrow><mi>g</mi><mo>:</mo><msup><mi>ℝ</mi><mo>+</mo></msup><mo>×</mo><msup><mi>ℝ</mi><mi>m</mi></msup><mo>→</mo><msup><mi>ℝ</mi><mi>m</mi></msup></mrow></math> and <math id="M3"><mrow><mi>f</mi><mo>:</mo><msup><mi>ℝ</mi><mo>+</mo></msup><mo>×</mo><msup><mi>ℝ</mi><mi>m</mi></msup><mo>→</mo><msup><mi>ℝ</mi><mrow><mi>m</mi><mo>×</mo><mi>m</mi></mrow></msup></mrow></math> are locally bounded and measurable functions, W( t ), is an m-dimensional Wiener process which is defined as the diffusion term. The general solution of Eq. (1) can be written as: <math display="block" id="M4"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><mi>g</mi><mo stretchy="false">(</mo><mi>u</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext>d</mtext><mi>u</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><mi>f</mi><mo stretchy="false">(</mo><mi>u</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> Here, the integral <math id="M5"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup><mi>g</mi></mrow></mstyle><mo stretchy="false">(</mo><mi>u</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mi>d</mi><mi>u</mi></mrow></math> is typically the Riemann-Lebesgue integral, while the integral <math id="M6"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><mi>f</mi><mo stretchy="false">(</mo><mi>u</mi><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>X</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>u</mi><mo stretchy="false">)</mo></mrow></math> is considered as an Itô integral. Generally, the analytical solution of the multi-dimensional Eq.(2) to a great or significant extent is not possible, and numerical techniques are robustly utilized. 2.2 Wiener process 2 Recall that the standard Wiener process is a continuous-time stochastic process which is also called the standard Brownian motion. The Wiener process W ( t ) over [0, T ] is a random variable, W ( t ), that depends on a continuous time, t , and is characterized by three conditions as follows [ 17 17 ] 1. W ( t )=0 for t =0 2. W ( t ) has independent increments with <math id="M7"><mrow><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><mi>W</mi><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>∼</mo><msqrt><mrow><mi>t</mi><mo>−</mo><mi>s</mi></mrow></msqrt><mi mathvariant="script">N</mi><mo stretchy="false">(</mo><mn>0,1</mn><mo stretchy="false">)</mo></mrow></math> for <math id="M8"><mrow><mi>s</mi><mo>∈</mo><mo stretchy="false">[</mo><mn>0,</mn><mi>t</mi><mo stretchy="false">]</mo></mrow></math> ,where 𝒩 (0,1) is the normal distribution with zero mean and unit variance. 3. For 0 ≤ s 1 < t 1 < s 2 < t 2 ≤ T the increments W ( t 1 )- W ( s 1 ) and W ( t 2 )- W ( s 2 ) are independent random variables. For numerical computational object, it is helpful to assume the discretized Brownian motion, where W ( t ) is specialized at discrete values, t . Accordingly, for some positive integer, N , consider h = T/N and let W i denote W ( t i ) with t i = i h . The first condition states that W 0 =0 with the probability equal to one and the other conditions state that W i = W i -1+d W i , i =1, 2,..., N , where each ∂ W i is an independent random variable of the form <math id="M9"><mrow><msqrt><mi>h</mi></msqrt><mi mathvariant="script">N</mi><mo stretchy="false">(</mo><mn>0,1</mn><mo stretchy="false">)</mo></mrow></math> . 2.3 Fractional calculus 2 The Riemann-Liouville fractional integral and Caputo derivative are defined respectively as follows: Definition 1: [ 18 18 ] Let g ( t ) be a continues function, α ≥ 0 and t є R . The Riemann-Liouville integral can be defined as: <math display="block" id="M10"><mrow><msubsup><mi>I</mi><mi>x</mi><mi>α</mi></msubsup><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>x</mi></msubsup></mrow></mstyle><msup><mrow><mo stretchy="false">(</mo><mi>x</mi><mo>−</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><mi>g</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>t</mi><mn>,</mn></mrow></math> where, Γ( α ) is the gamma function of the fractional order. Definition 2 [ 19 19 ]: Let n be an integer number, where n -1 < α < n and x > 0. The Caputo fractional derivative of order α for a function g ( x ) is defined as <math display="block" id="M11"><mrow><msubsup><mi>D</mi><mi>x</mi><mi>α</mi></msubsup><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>x</mi></msubsup></mrow></mstyle><msup><mrow><mo stretchy="false">(</mo><mi>x</mi><mo>−</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mi>n</mi><mo>−</mo><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><msup><mi>g</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>t</mi></mrow></math> where <math id="M12"><mrow><msup><mi>g</mi><mrow><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msup><mo>∂</mo><mi>n</mi></msup><mi>g</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><msup><mi>t</mi><mi>n</mi></msup></mrow></mfrac></mrow></math> . The relation between the Caputo fractional derivative [ 20 20 ] and the Riemann-Liouville fractional integral is introduced in the following formula: <math display="block" id="M13"><mrow><msubsup><mi>I</mi><mi>x</mi><mi>α</mi></msubsup><msubsup><mi>D</mi><mi>x</mi><mi>α</mi></msubsup><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>=</mo><mi>g</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>r</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover></mstyle><msup><mi>g</mi><mrow><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></msup><msub><mrow><mo stretchy="false">(</mo><mn>0</mn></mrow><mo>+</mo></msub><mo stretchy="false">)</mo><mfrac><mrow><msup><mi>x</mi><mi>r</mi></msup></mrow><mrow><mi>r</mi><mo>!</mo></mrow></mfrac><mtext> </mtext><mtext> </mtext><mtext> </mtext><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>n</mi><mo>−</mo><mn>1</mn><mo><</mo><mi>α</mi><mo>≤</mo><mi>n</mi><mn>.</mn></mrow></math> 3 Developed mathematical technique 1 The fractional stochastic model of the point reactor kinetics equation can be written as: [ 1 1 , 2 2 , 21 21 , 22 22 ] <math display="block" id="M14"><mrow><msubsup><mi>D</mi><mi>t</mi><mi>α</mi></msubsup><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mo>+</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><msub><mi>D</mi><mi>t</mi></msub><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mn>,</mn></mrow></math> where, <math display="block" id="M15"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>c</mi><mi>I</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><mi>ρ</mi><mo>−</mo><mi>β</mi></mrow><mi>Λ</mi></mfrac></mrow></mtd><mtd><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><msub><mi>λ</mi><mn>2</mn></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><msub><mi>λ</mi><mi>I</mi></msub></mrow></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>1</mn></msub></mrow><mi>Λ</mi></mfrac></mrow></mtd><mtd><mrow><mo>−</mo><msub><mi>λ</mi><mn>1</mn></msub></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>2</mn></msub></mrow><mi>Λ</mi></mfrac></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mrow><mo>−</mo><msub><mi>λ</mi><mn>2</mn></msub></mrow></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mi>I</mi></msub></mrow><mi>Λ</mi></mfrac></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd><mtd><mrow><mo>−</mo><msub><mi>λ</mi><mi>I</mi></msub></mrow></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mi>q</mi></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> t is the time, n ( t ) is the neutron population, ρ is the total reactivity, c i ( t ) is the i -group of delayed precursor concentration, λ i is the decay constant of i -group of delayed neutrons, β i is the fraction of i -group delayed neutrons, <math id="M16"><mrow><mi>β</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><msub><mi>β</mi><mi>i</mi></msub></mrow></math> is the total fraction of delayed neutrons, I is the total number of delayed neutron groups, Λ is the prompt neutron generation time, q is the external source of neutrons, http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10329866&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10329865&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10329866&type= and <math id="M17"><mrow><mo>|</mo><mi>Δ</mi><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msqrt><mi>h</mi></msqrt><mo>|</mo><mi>η</mi><mo>〉</mo></mrow></math> such that W 1 ( t ), W 2 ( t ),..., W I +1 ( t ) are the Wiener processes [ 5 5 , 17 17 ]. Notice that the initial condition can be determined as n (0)= n 0 , <math id="M18"><mrow><mfrac><mrow><mtext>d</mtext><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><msub><mo>|</mo><mrow><mi>t</mi><mo>=</mo><mn>0</mn></mrow></msub><mo>=</mo><mn>0</mn></mrow></math> , and <math id="M19"><mrow><mfrac><mrow><mtext>d</mtext><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><msub><mo>|</mo><mrow><mi>t</mi><mo>=</mo><mn>0</mn></mrow></msub><mo>=</mo><mn>0,</mn><mi>i</mi><mo>=</mo><mn>1,2,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi></mrow></math> . This means that <math display="block" id="M20"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>n</mi><mn>0</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>1</mn></msub><msub><mi>n</mi><mn>0</mn></msub></mrow><mrow><mi>Λ</mi><msub><mi>λ</mi><mn>1</mn></msub></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>2</mn></msub><msub><mi>n</mi><mn>0</mn></msub></mrow><mrow><mi>Λ</mi><msub><mi>λ</mi><mn>2</mn></msub></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mi>I</mi></msub><msub><mi>n</mi><mn>0</mn></msub></mrow><mrow><mi>Λ</mi><msub><mi>λ</mi><mi>I</mi></msub></mrow></mfrac></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> Of course, when the variance matrix B=0 and the fractional derivative order α =1, Eq. (6) is reduced to the standard point reactor kinetics model. In what follows, we aim to solve the fractional stochastic differential Eq. (6). The technique of split-step was utilized, that is, Eq. (6) is separated into deterministic and stochastic parts, followed by solving each of them separately. Applying the Laplace transformation on the deterministic part, i.e. consider B=0, of Eq. (6) is as follows: [ 23 23 - 25 25 23-25 ] <math display="block" id="M21"><mrow><msup><mi>s</mi><mi>α</mi></msup><mo>|</mo><mi mathvariant="script">P</mi><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>−</mo><msup><mi>s</mi><mrow><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo>|</mo><mi mathvariant="script">P</mi><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mfrac><mn>1</mn><mi>s</mi></mfrac><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo></mrow></math> where <math id="M22"><mrow><mo>|</mo><mi mathvariant="script">P</mi><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mi mathvariant="script">L</mi><mo stretchy="false">[</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo stretchy="false">]</mo></mrow></math> , <math id="M23"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo></mrow></math> is the solution of the deterministic part for Eq. (6). Consequently, we get: <math display="block" id="M24"><mrow><mo>|</mo><mi mathvariant="script">P</mi><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msup><mi>s</mi><mrow><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><msup><mrow><mo stretchy="false">[</mo><msup><mi>s</mi><mi>α</mi></msup><mstyle mathvariant="bold" mathsize="normal"><mi>I</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">]</mo></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mfrac><mn>1</mn><mi>s</mi></mfrac><msup><mrow><mo stretchy="false">[</mo><msup><mi>s</mi><mi>α</mi></msup><mstyle mathvariant="bold" mathsize="normal"><mi>I</mi></mstyle><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">]</mo></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mn>.</mn></mrow></math> Using the inverse Laplace transformation, we have: <math display="block" id="M25"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msub><mi>E</mi><mrow><mi>α</mi><mn>,1</mn></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>t</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><msup><mrow><mo stretchy="false">(</mo><mi>t</mi><mo>−</mo><mi>τ</mi><mo stretchy="false">)</mo></mrow><mrow><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><msub><mi>E</mi><mrow><mi>α</mi><mn>,</mn><mi>α</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mrow><mo stretchy="false">(</mo><mi>t</mi><mo>−</mo><mi>τ</mi><mo stretchy="false">)</mo></mrow><mi>α</mi></msup><mo stretchy="false">)</mo><mtext>d</mtext><mi>τ</mi><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mn>.</mn></mrow></math> Let us introduce the parameter z = t - τ into Eq.(12) to have: <math display="block" id="M26"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msub><mi>E</mi><mrow><mi>α</mi><mn>,1</mn></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>t</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><msup><mi>z</mi><mrow><mi>α</mi><mo>−</mo><mn>1</mn></mrow></msup><msub><mi>E</mi><mrow><mi>α</mi><mn>,</mn><mi>α</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>z</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mtext>d</mtext><mi>z</mi><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mn>.</mn></mrow></math> Using the integration property of the Mittag-Leffler function which is: <math display="block" id="M27"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><msup><mi>z</mi><mrow><mi>b</mi><mo>−</mo><mn>1</mn></mrow></msup><msub><mi>E</mi><mrow><mi>a</mi><mn>,</mn><mi>b</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>z</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mtext>d</mtext><mi>z</mi><mo>=</mo><msup><mi>t</mi><mi>b</mi></msup><msub><mi>E</mi><mrow><mi>a</mi><mn>,</mn><mi>b</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>t</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mn>,</mn></mrow></math> for b > 0, then Eq. (13) reads as follows: <math display="block" id="M28"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msub><mi>E</mi><mrow><mi>α</mi><mn>,1</mn></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>t</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><msup><mi>t</mi><mi>α</mi></msup><msub><mi>E</mi><mrow><mi>α</mi><mn>,</mn><mi>α</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>t</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mn>.</mn></mrow></math> Equation (15) introduces the general solution of the fractional stochastic point reactor kinetics equations, which depends on the stiff coefficient matrix A. To overcome the stiffness of this matrix, the coefficient matrix A was changed by its eigenvalues, ω j , [ 26 26 - 28 28 26-28 ] and the corresponding eigenvectors, |V i ⟩, of the matrix A [ 29 29 - 32 32 29-32 ]. Furthermore, over a small time interval with step size h , the matrices B and A are considered constant over the specified time interval [ t m , t m +1 ] where t m +1= t m + h and m =0,1,2,..., M -1. As a results of this substitution, the Mittag-Leffler function can be written as: <math display="block" id="M29"><mrow><msub><mi>E</mi><mrow><mi>a</mi><mn>,</mn><mi>b</mi></mrow></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msup><mi>z</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>I</mi></munderover></mstyle><msub><mi>E</mi><mrow><mi>a</mi><mn>,</mn><mi>b</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><msup><mi>z</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><mn>,</mn></mrow></math> where, the coefficient matrix A satisfies the following property in bra-ket space [ 33 33 ] <math display="block" id="M30"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>=</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mi>T</mi></msup><mo>=</mo><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><msub><mi>ω</mi><mi>j</mi></msub><mn>,</mn></mrow></math> and <math display="block" id="M31"><mrow><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>l</mi></msub><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>=</mo><msub><mi>δ</mi><mrow><mi>l</mi><mn>,</mn><mi>j</mi></mrow></msub><mo>=</mo><mrow><mo>{</mo><mtable columnalign="left"><mtr><mtd><mn>1</mn><mo>,</mo><mtext> </mtext><mi>l</mi><mo>=</mo><mi>j</mi></mtd></mtr><mtr><mtd><mn>0</mn><mo>,</mo><mtext> </mtext><mi>l</mi><mo>≠</mo><mtext> </mtext><mi>j</mi></mtd></mtr></mtable></mrow><mn>.</mn></mrow></math> The ket eigenvectors, |V j ⟩, can be written by the following analytical form [ 33 33 ] <math display="block" id="M32"><mrow><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>=</mo><msub><mi>σ</mi><mi>j</mi></msub><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>1</mn></msub></mrow><mrow><mi>Λ</mi><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mn>1</mn></msub><mo stretchy="false">)</mo></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>2</mn></msub></mrow><mrow><mi>Λ</mi><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mn>2</mn></msub><mo stretchy="false">)</mo></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mi>I</mi></msub></mrow><mrow><mi>Λ</mi><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mi>I</mi></msub><mo stretchy="false">)</mo></mrow></mfrac></mrow></mtd></mtr><mtr><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow></mrow></math> and the bra eigenvectors, 〈U j |, is <math display="block" id="M33"><mrow><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><mo>=</mo><msub><mi>σ</mi><mi>j</mi></msub><mrow><mo>(</mo><mrow><mn>1</mn><mtext> </mtext><mfrac><mrow><msub><mi>λ</mi><mn>1</mn></msub></mrow><mrow><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mn>1</mn></msub><mo stretchy="false">)</mo></mrow></mfrac><mfrac><mrow><msub><mi>λ</mi><mn>2</mn></msub></mrow><mrow><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mn>2</mn></msub><mo stretchy="false">)</mo></mrow></mfrac><mo>⋯</mo><mfrac><mrow><msub><mi>λ</mi><mi>I</mi></msub></mrow><mrow><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mi>I</mi></msub><mo stretchy="false">)</mo></mrow></mfrac></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> From the normalization ⟨U j |V j ⟩=1, we can deduce: <math display="block" id="M34"><mrow><msub><mi>σ</mi><mi>j</mi></msub><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub><msub><mi>λ</mi><mi>i</mi></msub></mrow><mrow><mi>Λ</mi><msup><mrow><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><mo>+</mo><msub><mi>λ</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mo>−</mo><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mn>,</mn><mtext> </mtext><mo>∀</mo><mi>j</mi><mo>=</mo><mn>0,1,2,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mn>.</mn></mrow></math> After introducing Eq. (16) into Eq. (15) we get: <math display="block" id="M35"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo></mrow></mtd><mtd columnalign="left"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>I</mi></munderover></mstyle><mo stretchy="false">[</mo><msub><mi>E</mi><mrow><mi>α</mi><mn>,1</mn></mrow></msub><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><msup><mi>h</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>〉</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msup><mi>h</mi><mi>α</mi></msup><msub><mi>E</mi><mrow><mi>α</mi><mn>,</mn><mi>α</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">(</mo><msub><mi>ω</mi><mi>j</mi></msub><msup><mi>h</mi><mi>α</mi></msup><mo stretchy="false">)</mo><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>V</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>|</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>Q</mi></mstyle><mo>〉</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The obtained results, <math id="M36"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>〉</mo></mrow></math> , from Eq. (22), were used to evaluate the variance matrix <math id="M37"><mrow><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>B</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></math> . Furthermore, according to the split-step method, the general solution of Eq. (6) is given by: <math display="block" id="M38"><mrow><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mo>|</mo><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>P</mi><mo>¯</mo></mover></mstyle><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><msqrt><mi>h</mi></msqrt><msup><mstyle mathvariant="bold" mathsize="normal"><mover accent="true"><mi>B</mi><mo>¯</mo></mover></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>|</mo><mi>η</mi><mo>〉</mo><mn>,</mn></mrow></math> where <math id="M39"><mrow><mo>|</mo><mi>Δ</mi><mstyle mathvariant="bold" mathsize="normal"><mi>W</mi></mstyle><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msqrt><mi>h</mi></msqrt><mo>|</mo><mi>η</mi><mo>〉</mo></mrow></math> . Equation (23) represents the general solution of the fractional stochastic point kinetics model with a multi-group of delayed precursor concentration. 4 Computational Results 1 In order to affirm the exactness and validity of the proposed technique over the traditional methods, the developed mathematical technique (DMT) for the fractional stochastic point kinetics model with a multi-group of delayed precursor concentration is tested through Matlab code. The mean and standard deviation of the neutron populations and precursor concentrations in different cases of reactivities, step, ramp, sinusoidal, and also in the presence of temperature feedback, are calculated. The results of the proposed method are compared to five different values of fractional order with the results of the traditional stochastic methods. 4.1 Step reactivity 2 The first computational example simulates two cases of step reactivity ρ =0.007 and ρ =0.003 for an actual reactor with six groups of delayed precursors [ 1 1 ]. The parameters of the two problems are taken from the following references [ 1 1 , 5 5 , 9 9 , 11 11 , 12 12 ] as follows: ν =2.5, Λ=0.00002 s, n 0 =100 (neutrons), λ i =[0.0127, 0.0317, 0.115, 0.311, 1.4, 3.87] s -1 , β i =[0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182], β =0.007, and q =0. The mean and standard deviation at step size h =0.001 s after 500 trails are tabulated in Table 1 Table 1 for reactivity ρ =0.007 and time T =0.001 s. In the calculation of the developed mathematical technique, the partial requests of the fractional order are taken respectively as: α =1, α =0.98, α =0.99, α =1.01, and α =1.02. The Matlab code is tested with the conventional results under the same conditions. The obtained results of the developed mathematical technique (DMT) are compared with the familiar Monte Carlo (MC) [ 1 1 ], stochastic piecewise constant approximation (SPCA) [ 1 1 ], Euler-Maruyama (EM) [ 5 5 ], Taylor 1.5 strong order (T1.5SO) [ 5 5 ], analytical exponential model (AEM) [ 12 12 ], efficient stochastic model (ESM) [ 11 11 ], and the deterministic point kinetics model (DPKM) [ 33 33 ]. 1 1   1 1 Method 1 1 Step size 1 1 α 1 1 E [ n ( t )] 1 1 σ [ n ( t )] 1 1 <math id="M40"><mrow><mi>E</mi><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mn>6</mn></munderover><mrow><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 <math id="M41"><mrow><mi>σ</mi><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mn>6</mn></munderover><mrow><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 2 1 ρ =0.007 t =0.001 1 1   1 1 MC 1 1   1 1 1.0 1 1 135.67 1 1 93.376 1 1 4.464×10 5 1 1 7.8073 1 1   1 1 SPCA 1 1   1 1 1.0 1 1 134.55 1 1 91.242 1 1 4.464×10 5 1 1 19.444 1 1   1 1 EM 1 1   1 1 1.0 1 1 139.57 1 1 92.042 1 1 4.463×10 5 1 1 6.071 1 1   1 1 T1.5SO 1 1   1 1 1.0 1 1 139.57 1 1 92.047 1 1 4.463×10 5 1 1 18.337 1 1   1 1 AEM 1 1 h =0.00001 1 1 1.0 1 1 134.54 1 1 91.234 1 1 4.464×10 5 1 1 19.235 1 1   1 1 ESM 1 1 h =0.00001 1 1 1.0 1 1 134.96 1 1 6.8527 1 1 4.464×10 5 1 1 2.529 1 1   1 1 DPKM 1 1 h =0.000025 1 1 1.0 1 1 135.0 1 1   1 1 4.464×10 5 1 1   1 1 DMT 1 1 h =0.001 1 1 0.98 1 1 140.122221 1 1 7.872980 1 1 446362.78 1 1 3.794234 1 1   1 1 DMT 1 1 h =0.001 1 1 0.99 1 1 137.268228 1 1 7.746689 1 1 446361.57 1 1 3.657256 1 1   1 1 DMT 1 1 h =0.001 1 1 1.0 1 1 134.613600 1 1 7.627336 1 1 446360.52 1 1 3.525117 1 1   1 1 DMT 1 1 h =0.001 1 1 1.01 1 1 132.144579 1 1 7.514619 1 1 446359.63 1 1 3.397625 1 1   1 1 DMT 1 1 h =0.001 1 1 1.02 1 1 129.848341 1 1 7.408245 1 1 446358.87 1 1 3.274607 1 1   1 1 DMT 1 1 h =0.0005 1 1 0.98 1 1 141.287614 1 1 7.427251 1 1 446362.79 1 1 3.262914 1 1   1 1 DMT 1 1 h =0.0005 1 1 0.99 1 1 138.120575 1 1 7.318368 1 1 446361.49 1 1 3.137738 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.0 1 1 135.195181 1 1 7.216352 1 1 446360.39 1 1 3.017724 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.01 1 1 132.493180 1 1 7.120848 1 1 446359.45 1 1 2.902478 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.02 1 1 129.997682 1 1 7.031509 1 1 446358.66 1 1 2.791808 1 1   1 1 DMT 1 1 h =0.0001 1 1 0.98 1 1 142.454816 1 1 7.029539 1 1 446363.37 1 1 2.946374 1 1   1 1 DMT 1 1 h =0.0001 1 1 0.99 1 1 138.562403 1 1 6.927727 1 1 446361.78 1 1 2.796684 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.0 1 1 135.024891 1 1 6.833936 1 1 446360.47 1 1 2.656922 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.01 1 1 131.810129 1 1 6.747624 1 1 446359.39 1 1 2.525245 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.02 1 1 128.888862 1 1 6.668273 1 1 446358.50 1 1 2.401280 2 1 ρ =0.003 t =0.1 1 1   1 1 MC 1 1   1 1 1.0 1 1 183.04 1 1 168.79 1 1 4.478×10 5 1 1 1495.7 1 1   1 1 SPCA 1 1   1 1 1.0 1 1 186.31 1 1 164.16 1 1 4.491×10 5 1 1 1917.20 1 1   1 1 EM 1 1   1 1 1.0 1 1 208.60 1 1 255.95 1 1 4.498×10 5 1 1 1233.38 1 1   1 1 T1.5SO 1 1   1 1 1.0 1 1 199.41 1 1 168.547 1 1 4.497×10 5 1 1 1218.82 1 1   1 1 AEM 1 1 h =0.001 1 1 1.0 1 1 186.30 1 1 164.14 1 1 4.490×10 5 1 1 1911.91 1 1   1 1 ESM 1 1 h =0.001 1 1 1.0 1 1 179.93 1 1 10.555 1 1 4.489×10 5 1 1 94.75 1 1   1 1 DPKM 1 1 h =0.0025 1 1 1.0 1 1 179.95 1 1   1 1 4.489×10 5 1 1   1 1 DMT 1 1 h =0.001 1 1 0.98 1 1 180.922922 1 1 12.847182 1 1 449304.10 1 1 94.953085 1 1   1 1 DMT 1 1 h =0.001 1 1 0.99 1 1 180.466893 1 1 13.190199 1 1 449083.17 1 1 94.453581 1 1   1 1 DMT 1 1 h =0.001 1 1 1.0 1 1 180.036327 1 1 13.547699 1 1 448878.15 1 1 93.938513 1 1   1 1 DMT 1 1 h =0.001 1 1 1.01 1 1 179.629458 1 1 13.919605 1 1 448687.88 1 1 93.403821 1 1   1 1 DMT 1 1 h =0.001 1 1 1.02 1 1 179.244639 1 1 14.305915 1 1 448511.28 1 1 92.845214 1 1   1 1 DMT 1 1 h =0.0005 1 1 0.98 1 1 181.064354 1 1 11.889735 1 1 449347.96 1 1 95.974066 1 1   1 1 DMT 1 1 h =0.0005 1 1 0.99 1 1 180.562722 1 1 12.307332 1 1 449103.51 1 1 95.397178 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.0 1 1 180.093120 1 1 12.742279 1 1 448878.30 1 1 94.806707 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.01 1 1 179.653284 1 1 13.194207 1 1 448670.78 1 1 94.198372 1 1   1 1 DMT 1 1 h =0.0005 1 1 1.02 1 1 179.241129 1 1 13.662739 1 1 448479.56 1 1 93.567468 1 1   1 1 DMT 1 1 h =0.0001 1 1 0.98 1 1 180.596300 1 1 10.464344 1 1 449450.09 1 1 88.949000 1 1   1 1 DMT 1 1 h =0.0001 1 1 0.99 1 1 179.988863 1 1 10.964650 1 1 449149.09 1 1 88.353100 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.0 1 1 179.435711 1 1 11.490519 1 1 448876.38 1 1 87.735258 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.01 1 1 178.932909 1 1 12.042888 1 1 448629.26 1 1 87.089442 1 1   1 1 DMT 1 1 h =0.0001 1 1 1.02 1 1 178.476566 1 1 12.622634 1 1 448405.28 1 1 86.409072 In a similar manner, the results of the second case ρ =0.003 are given in Table 1 Table 1 . The numerical results, in this case, are calculated at time t =0.1 s. Moreover, the pattern of the two individual neutron sample paths and the mean neutron population are shown in Fig. 1 Fig. 1 at different values of the fractional order α =0.98, α =1, and α =1.02, symbolized a , b , and c respectively. The red solid curve is the mean neutron population, while the blue dashed curve and black dot-dashed curve represent the two individual neutron sample paths, which show approximately a real behavior for neutron density into the actual reactors. 4.2 Ramp reactivity 2 The second example simulates a ramp reactivity, where the same parameters are taken from the previous example. The function of reactivity is taken as: ρ =0.1 β t and T =1 s. The result of a comparison of the proposed method (DMT), efficient stochastic model (ESM) [ 11 11 ], and the deterministic point kinetics model (DPKM) [ 33 33 ] is given in Table 2 Table 2 . To show the effect of the time interval, the code for the developed method was run at different time step as well as at different values of the parameter α . The relative percentage errors for E [ n ( t )] with the DMT is recorded as 0.025, 0.026, and 0.159 in compared with the conventional methods SPCA, AEM, and ESM respectively under the same conditions at time step h =0.01 s, Table 2 Table 2 for ramp reactivity, while the relative percentage error is 0.05 for the DPKM at time step h =0.001 s. Furthermore, in the same table, the relative percentage errors for the DMT method at different step time interval ( h =0.01 and h =0.001) are recorded as 0.126, while at h =0.01 and h =0.0001 it is 0.322. In Table 1 Table 1 for step reactivity the effect of the time step interval appears from the relative percentage error at different time steps, e.g. at h =0.001 and h =0.0005 the RPE is -0.432, while at h =0.001 and h =0.0001 it is -0.306. The previous analysis confirms the stability of the DMT method and the effect of time step is acceptable. The behavior of the results show that an increase with the decreasing value of α and vise versa. Finally, we conclude that the validity of the proposed method shows a high agreement with the deterministic method DPKM as well as with the conventional methods (e.g. MC, SPCA, EM, T1.5SO, AEM, ESM, and DMT). 1 1 Method 1 1 Step size 1 1 α 1 1 E [ n (1)] 1 1 σ [ n (1)] 1 1 <math id="M42"><mrow><mi>E</mi><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mn>6</mn></munderover><mrow><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 <math id="M43"><mrow><mi>σ</mi><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mn>6</mn></munderover><mrow><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 SPCA 1 1 h =0.01 1 1 1.0 1 1 113.268077 1 1 13.330142 1 1 448239.846 1 1 3009.93141 1 1 AEM 1 1 h =0.01 1 1 1.0 1 1 113.267707 1 1 13.327291 1 1 448239.798 1 1 3002.68282 1 1 ESM 1 1 h =0.01 1 1 1.0 1 1 113.116433 1 1 4.111150 1 1 448253.780 1 1 47.203115 1 1 DPKM 1 1 h =0.001 1 1 1.0 1 1 113.091124 1 1   1 1 448236.26 1 1 DMT 1 1 h =0.01 1 1 0.98 1 1 113.474956 1 1 9.476546 1 1 448434.483 1 1 46.871815 1 1 DMT 1 1 h =0.01 1 1 0.99 1 1 113.384203 1 1 9.473676 1 1 448331.161 1 1 46.840134 1 1 DMT 1 1 h =0.01 1 1 1.0 1 1 113.296719 1 1 9.470503 1 1 448232.846 1 1 46.811703 1 1 DMT 1 1 h =0.01 1 1 1.01 1 1 113.212358 1 1 9.468112 1 1 448139.299 1 1 46.782352 1 1 DMT 1 1 h =0.01 1 1 1.02 1 1 113.131034 1 1 9.466234 1 1 448050.294 1 1 46.753327 1 1 DMT 1 1 h =0.001 1 1 0.98 1 1 113.415094 1 1 3.971593 1 1 448535.015 1 1 43.751183 1 1 DMT 1 1 h =0.001 1 1 0.99 1 1 113.280768 1 1 4.058934 1 1 448378.858 1 1 43.695078 1 1 DMT 1 1 h =0.001 1 1 1.0 1 1 113.153337 1 1 4.150857 1 1 448233.774 1 1 43.632678 1 1 DMT 1 1 h =0.001 1 1 1.01 1 1 113.032287 1 1 4.247468 1 1 448098.976 1 1 43.555908 1 1 DMT 1 1 h =0.001 1 1 1.02 1 1 112.917106 1 1 4.348945 1 1 447973.736 1 1 43.348216 1 1 DMT 1 1 h =0.0001 1 1 0.98 1 1 113.305430 1 1 3.263194 1 1 448638.295 1 1 47.552005 1 1 DMT 1 1 h =0.0001 1 1 0.99 1 1 113.112671 1 1 3.411900 1 1 448425.709 1 1 47.488998 1 1 DMT 1 1 h =0.0001 1 1 1.0 1 1 112.932052 1 1 3.567299 1 1 448232.798 1 1 47.438442 1 1 DMT 1 1 h =0.0001 1 1 1.01 1 1 112.762814 1 1 3.729720 1 1 448057.756 1 1 47.366197 1 1 DMT 1 1 h =0.0001 1 1 1.02 1 1 112.604169 1 1 3.899621 1 1 447898.902 1 1 47.283241 For various values of fractional orders ( α =0.98, α =1, α =1.02), the pattern of the two individual neutron sample paths and the mean neutron are shown in Figs. 2 Figs. 2 , 3 3 , 4 4 . Furthermore, the intensity of fluctuations for neutron sample paths increases with the mean neutron population. This phenomenon arises from the fact that the variance matrix including are stochastic part is dependent on the mean neutron population and precursors population. 4.3 Sinusoidal reactivity 2 In the third example, the variation of reactivity insertion in the form of a sinusoidal change as <math id="M44"><mrow><mi>ρ</mi><mo>=</mo><msub><mi>ρ</mi><mn>0</mn></msub><mi>sin</mi><mrow><mo>(</mo><mrow><mfrac><mrow><mi>π</mi><mi>t</mi></mrow><mi>T</mi></mfrac></mrow><mo>)</mo></mrow></mrow></math> [ 6 6 ] is considered. The parameters of the nuclear reactor with one group of delayed neutrons are as follows: ρ 0 =0.005333 (0.68 $), β 1 = β =0.0079, λ 1 =0.077, Λ=10 -3 , q =0, n 0 =1 (neutrons) with period time of T =100 s. The two individual neutron sample paths and the mean neutron are shown in Figs. 5 Figs. 5 , 6 6 , and 7 7 using different values of fractional orders ( α =0.98, α =1, α =1.02) respectively, where the time step interval is h =0.01 s after 500 trails. We compared the stochastic and deterministic solutions to deduce the fact that the deterministic solution represents the average of the stochastic approach. Moreover, the individual sample paths oscillates around the deterministic curve. The intensity of fluctuations in these classes of figures arise from the variations of the mean neutron population. 4.4 Temperature feedback reactivity 2 In the most nuclear literature, there are two cases for the external reactivity, step, and ramp external reactivities. The means of the neutron population are calculated for a U 235 nuclear reactor with step and ramp external reactivities. In what follows, the effect of Newtonian temperature feedback introduced into the reactivity is analyzed. The new reactivity form in the presence of temperature feedback is given by: <math display="block" id="M48"><mrow><mi>ρ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><mi>σ</mi><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mn>0</mn><mi>t</mi></msubsup></mrow></mstyle><mi>N</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>τ</mi></mrow></math> where, ρ ex ( t ) represents the external reactivity and σ is the nonlinear coefficient part which represents the product of the reciprocal of the thermal capacity and the temperature coefficient. 4.4.1 Step external reactivity 3 To check the developed mathematical technique for the nonlinear fractional stochastic model in the presence of temperature feedback and step external reactivity, let us take the parameters of the U 235 nuclear reactor as follows [ 14 14 , 15 15 , 34 34 ]: λ i =[0.0124, 0.0305, 0.111, 0.301, 1.13, 3.0] s -1 , β i =[0.00021, 0.00141, 0.00127, 0.00255, 0.00074, 0.00027, 0.00645], β =0.00645, Λ=5.0× 10 -5 s, σ =2.5× 10 -6 MWs) -1 , and N (0)=1 (neutrons). In Table 3 Table 3 , the peak of the mean neutron population with the corresponding time at various step external reactivity, ρ ex =1.0$, ρ ex =1.5$, and ρ ex =2$, are tabulated for different values of fractional order 0.98, 0.99, 1.0, 1.01, and 1.02. Using time step h =0.001 s and after 500 trails, the peak of the mean neutron population is compared with the peak of the mean neutron population using the analytical exponential technique (AET) [ 14 14 ]. In addition, the mean neutron population and the two individual neutron sample paths are plotted in Figs. 8 Figs. 8 , 9 9 , and 10 10 for the step external reactivity ( ρ ex =0.5$) and Figs. 11 Figs. 11 , 12 12 , and 13 13 for the step external reactivity ( ρ ex =0.75$) using different values of fractional derivative order 0.98,1.0, 1.02, respectively. The most important notice is that the mean neutron population increases with time until it reaches the maximum value due to the positive external reactivity. Thereafter, the mean neutron tends to zero with the increasing time due to the effect of temperature reactivity feedback. Again, these figures confirm that the amplitude of fluctuations for neutron sample paths is affected by the direct variations of the mean neutron population. 1 1 Method 1 1 α 1 1 ρ ex ($) 1 1 0.5 1 1 0.75 1 1 1.0 1 1 1.5 1 1 2.0 1 1 DMT 1 1 0.98 1 1 Peak 1 1 53.0279 1 1 189.727 1 1 895.547 1 1 41083.125 1 1 158278.694 1 1   1 1   1 1 Time 1 1 25.378 1 1 7.732 1 1 0.942 1 1 0.151 1 1 0.088 1 1   1 1 0.99 1 1 Peak 1 1 49.565 1 1 176.089 1 1 826.687 1 1 36557.357 1 1 140582.052 1 1   1 1   1 1 Time 1 1 26.284 1 1 8.181 1 1 1.025 1 1 0.162 1 1 0.094 1 1   1 1 1.0 1 1 Peak 1 1 46.289 1 1 163.986 1 1 768.739 1 1 34289.333 1 1 132756.323 1 1   1 1   1 1 Time 1 1 28.142 1 1 8.762 1 1 1.09 1 1 0.173 1 1 0.101 1 1   1 1 1.01 1 1 Peak 1 1 43.137 1 1 152.776 1 1 712.311 1 1 30306.984 1 1 117251.348 1 1   1 1   1 1 Time 1 1 30.303 1 1 9.429 1 1 1.211 1 1 0.186 1 1 0.108 1 1   1 1 1.02 1 1 Peak 1 1 40.239 1 1 142.525 1 1 662.795 1 1 28353.797 1 1 111465.515 1 1   1 1   1 1 Time 1 1 30.708 1 1 9.711 1 1 1.245 1 1 0.199 1 1 0.116 1 1 AET 1 1 1.0 1 1 Peak 1 1 46.235334 1 1 164.209516 1 1 770.5188 1 1 33119.58 1 1 128083.1 1 1   1 1   1 1 Time 1 1 28.142 1 1 8.867 1 1 1.040 1 1 0.174 1 1 0.101 4.4.2 Ramp external reactivity 3 In the test example for ramp external reactivity, the parameters for the U 235 nuclear reactor is taken as the same values from the nonlinear case, where Λ=10×10 -5 s, σ =10 -11 , or 10 -13 (MW s) -1 and the external reactivity is ramp ( ρ ex = 0.01 t , 0.1 t ). The peak of the mean neutron population and its time with various values of fractional order (0.98,1.0, 1.02) are compared with the results of the analytical exponential technique (AET) [ 14 14 ] in Table 4 Table 4 . The results are calculated with time step h =0.001 s and the number of trails is 500. In addition, two individual neutron sample paths and the mean neutron population are drawn for external reactivity ρ ex ( t )=0.003 t and the nonlinear coefficient σ =10 -11 in Figs. 14 Figs. 14 , 15 15 , and 16 16 and σ =10 -13 in Figs. 17 Figs. 17 , 18 18 , and 19 19 . In Figs. 14 Figs. 14 , 15 15 , 16 16 , 17 17 , 18 18 and 19 19 , the mean neutron increases with time until it reaches the peak due to the external reactivity increasing. After that the mean neutron decreases due to the effect of temperature feedback. Therefore, the mean neutron population is almost stable due to the effect of external reactivity, which is equivalent with the effect of temperature feedback. Moreover, the fluctuations of neutron sample paths disappear approximately around the sharp peak due to the logarithmic scale and a large increase in the neutron population at a very small time. Figure 14 Figure 14 shows a pattern of the effect of the neutron population and fluctuations at different sections with time. A sharp increase of fluctuation with a slight increase of the neutron population in the ramp section is followed by a slight increase of the intensity of fluctuations compared with a sharp increase of the neutron population. Moreover, the same remark is observed at the peak and the remainder section of the figures. 2 1 α 1 1   2 1 a =0.003 2 1 a =0.01 2 1 a =0.1 1 1   1 1   1 1   1 1 σ =10 -11 1 1 σ =10 -13 1 1 σ =10 -11 1 1 σ =10 -13 1 1 σ =10 -11 1 1 σ =10 -13 2 1 DMT 1 1   1 1 0.98 1 1 Peak 1 1 4.409101×10 9 1 1 5.771654×10 11 1 1 1.701579×10 10 1 1 2.144910×10 12 1 1 1.876758×10 11 1 1 2.242436×10 13 1 1   1 1   1 1 Time 1 1 2.446 1 1 2.488 1 1 0.831 1 1 0.852 1 1 0.130 1 1 0.136 1 1   1 1 0.99 1 1 Peak 1 1 4.452107×10 9 1 1 5.813629×10 11 1 1 1.704601×10 10 1 1 2.124988×10 12 1 1 1.918749×10 11 1 1 2.290324×10 13 1 1   1 1   1 1 Time 1 1 2.461 1 1 2.504 1 1 0.839 1 1 0.861 1 1 0.132 1 1 0.138 1 1   1 1 1.0 1 1 Peak 1 1 4.477012×10 9 1 1 5.828427×10 11 1 1 1.689973×10 10 1 1 2.102489×10 12 1 1 1.849657×10 11 1 1 2.146883×10 13 1 1   1 1   1 1 Time 1 1 2.475 1 1 2.520 1 1 0.847 1 1 0.869 1 1 0.135 1 1 0.142 1 1   1 1 1.01 1 1 Peak 1 1 4.505250×10 9 1 1 5.847340×10 11 1 1 1.673285×10 10 1 1 2.081893×10 12 1 1 1.860123×10 11 1 1 2.233038×10 13 1 1   1 1   1 1 Time 1 1 2.491 1 1 2.537 1 1 0.856 1 1 0.879 1 1 0.138 1 1 0.144 1 1   1 1 1.02 1 1 Peak 1 1 4.514868×10 9 1 1 5.844036×10 11 1 1 1.656222×10 10 1 1 2.059339×10 12 1 1 1.786014×10 11 1 1 2.131587×10 13 1 1   1 1   1 1 Time 1 1 2.507 1 1 2.554 1 1 0.864 1 1 0.888 1 1 0.141 1 1 0.148 2 1 AET 1 1   1 1 1.0 1 1 Peak 1 1 4.482853×10 9 1 1 5.833649×10 11 1 1 1.687565×10 10 1 1 2.099857×10 12 1 1 1.853996×10 11 1 1 2.217193×10 13 1 1   1 1   1 1 Time 1 1 2.476 1 1 2.521 1 1 0.847 1 1 0.870 1 1 0.135 1 1 0.141 General description and analysis for the results of the developed method can be summarized as follows. Figures 1 Figures 1 , 2 2 , 3 3 , 4 4 , 5 5 , 6 6 , and 7 7 show a large class of changes consists of small variations of the cross sections, positive or negative changes of reactivity attributed to the neutron population, around an expected value E [ n ( t )], corresponding to a critical case in which the perturbation can be induced flux fluctuations, where: δ n ( t )= n ( t )- E [ n ( t )] and δρ ( t )= ρ ( t )- E [ ρ ( t )]. Strong-frequency power fluctuations caused by increasing the reactivity (the fission rate is sufficiently large for the average value equations) of the reactor, as shown in Figs. 1 Figs. 1 , 2 2 , 3 3 , and 4 4 . In the opposite, low-frequency power fluctuations can arise from decreasing the reactivity (i.e. the fission rate) as shown in Figs. 5 Figs. 5 , 6 6 , 7 7 , 8 8 , 9 9 , 10 10 , 11 11 , 12 12 , and 13 13 . Furthermore, Figs. 14 Figs. 14 , 15 15 , 16 16 , 17 17 , 18 18 , and 19 19 deal with the variation of the neutron population, which obtained by studying the fluctuations arise from the temperature feedback. These figures conclude that the randomness of the input is communicated to the output via the response characteristics of the system, where the fluctuation is above or below the mean value at instant time. 5 Conclusion 1 The developed mathematical technique was presented for a linear/nonlinear fractional stochastic model of the point reactor kinetics system with a multi-group of delayed neutron precursors. This system is characterized by its stochastic behavior and can offer the only average or mean values of the modeled populations. However, the neutron population and the delayed neutron precursor concentrations vary randomly with time, meaning the real dynamical process is stochastic. This system was numerically implemented using a stochastic piecewise constant approximation (SPCA) due to the stiffness of these equations. In this paper, the matrix formula for this fractional stochastic model is solved through a developed mathematical technique, which is based on the split-step method, Laplace transforms, Mittage-Leffler function, eigenvalues, and eigenvectors. The mean and standard deviation of the neutron population and the sum of the precursor population were calculated for step, ramp, sinusoidal, and the temperature feedback reactivity insertion which represents the nonlinear fractional stochastic model. Moreover, this fractional differential system was calculated with different values of the fractional derivative order. In order to validate the proposed method (DMT), we present a comparison with the conventional results in the literature of the stochastic model and the deterministic point kinetics model, showing that the method is in agreement with those already established. The future work will be included the derivation and the study of a fractional stochastic model for the time-space kinetics equations. References J.G. Hayes , E.J. Allen , Stochastic point kinetics equations in nuclear reactor dynamics , Ann. Nucl. Energ . 32 , 572 - 587 ( 2005 ). doi: 10.1016/j.anucene.2004.11.009 http://doi.org/10.1016/j.anucene.2004.11.009 Stochastic point kinetics equations in nuclear reactor dynamics J.G. Hayes , Development of stochastic point kinetics equations in nuclear reactor dynamics, thesis , Texas Tech University ( 2005 ). https://ttu-ir.tdl.org/ttu-ir/bitstream/handle/2346/22276/HayesThesis.pdf https://ttu-ir.tdl.org/ttu-ir/bitstream/handle/2346/22276/HayesThesis.pdf Development of stochastic point kinetics equations in nuclear reactor dynamics, thesis P.N. Ha , J.K. Kim , A Stochastic Approach to Monoenergetic Space-Time Nuclear Reactor Kinetics , J. Nucl. Sci. Technol . 47 , 705 - 711 ( 2010 ). doi: 10.1080/18811248.2010.9711646 http://doi.org/10.1080/18811248.2010.9711646 A Stochastic Approach to Monoenergetic Space-Time Nuclear Reactor Kinetics P.N. Ha , J.K. Kim , Further evaluation of a stochatic model applied to monoenergetic space-time nuclear reactor kinetics , Nucl. Eng. Technol . 43 , 523 - 530 ( 2011 ). doi: 10.5516/NET.2011.43.6.523 http://doi.org/10.5516/NET.2011.43.6.523 Further evaluation of a stochatic model applied to monoenergetic space-time nuclear reactor kinetics S.S. Ray , Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor , Ann. Nucl. Energ . 49 , 154 - 159 ( 2012 ). doi: 10.1016/j.anucene.2012.05.022 http://doi.org/10.1016/j.anucene.2012.05.022 Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor S.S. Ray , A. Patra , Numerical solution for stochatic point kinetics equations with sinusoidal reactivity in dynamical system of nuclear reactor, Int . J. Nucl. Energ. Sci. Technol . 7 , 231 - 242 ( 2013 ). doi: 10.1504/IJNEST.2013.052165 http://doi.org/10.1504/IJNEST.2013.052165 Numerical solution for stochatic point kinetics equations with sinusoidal reactivity in dynamical system of nuclear reactor, Int E.J. Allen , Stochastic difference equations and a stochastic partial differential equation for neutron transport , J. Differ. Equ. Appl . 18 , 1267 - 1285 ( 2012 ). doi: 10.1080/10236198.2010.488229 http://doi.org/10.1080/10236198.2010.488229 Stochastic difference equations and a stochastic partial differential equation for neutron transport E.J. Allen , A stochastic analysis of power doubling time for a subcritical system, Stoch . Anal. Appl . 31 , 528 - 537 ( 2013 ). doi: 10.1080/07362994.2013.777287 http://doi.org/10.1080/07362994.2013.777287 A stochastic analysis of power doubling time for a subcritical system, Stoch S.S. Ray , A. Patra , Numerical solution of fractional stochastic neutron point kinetic equation for nuclear reactor dynamics , Ann. Nucl. Energ . 54 , 154 - 161 ( 2013 ). doi: 10.1016/j.anucene.2012.11.007 http://doi.org/10.1016/j.anucene.2012.11.007 Numerical solution of fractional stochastic neutron point kinetic equation for nuclear reactor dynamics S.M. Ayyoubzadeh , N. Vosoughi , An alternative stochastic formulation for the point kinetics , Ann. Nucl. Energ . 63 , 691 - 695 ( 2014 ). doi: 10.1016/j.anucene.2013.09.013 http://doi.org/10.1016/j.anucene.2013.09.013 An alternative stochastic formulation for the point kinetics A.A. Nahla , A.M. Edress , Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors , Nucl. Sci. Tech . 27 , 20 : 1 - 8 ( 2016 ). doi: 10.1007/s41365-016-0025-6 http://doi.org/10.1007/s41365-016-0025-6 Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors A.A. Nahla , A.M. Edress , Efficient stochastic model for the point kinetics equations, Stoch . Anal. Appl . 34 , 598 - 609 ( 2016 ). doi: 10.1080/07362994.2016.1159519 http://doi.org/10.1080/07362994.2016.1159519 Efficient stochastic model for the point kinetics equations, Stoch M.W. da Silva , R. Vasques , B.E.G. Bodmann et al ., A nonstiff solution for the stochastic neutron point kinetics equations , Ann. Nucl. Energ . 97 , 47 - 52 ( 2016 ). doi: 10.1016/j.anucene.2016.06.026 http://doi.org/10.1016/j.anucene.2016.06.026 A nonstiff solution for the stochastic neutron point kinetics equations A.A. Nahla , Stochastic model for the nonlinear point reactor kinetics equations in the presence Newtonian temperature feedback effects , J. Differ. Equ. Appl . 23 , 1003 - 1016 ( 2017 ). doi: 10.1080/10236198.2017.1308507 http://doi.org/10.1080/10236198.2017.1308507 Stochastic model for the nonlinear point reactor kinetics equations in the presence Newtonian temperature feedback effects S. Singh , S.S. Ray , On the comparison of two split-step methods for the numerical simulation of stochastic point kinetics equations in presence of Newtonian temperature feedback effects , Ann. Nucl. Energ ., 110 , 865 - 873 ( 2017 ). doi: 10.1016/j.anucene.2017.08.001 http://doi.org/10.1016/j.anucene.2017.08.001 On the comparison of two split-step methods for the numerical simulation of stochastic point kinetics equations in presence of Newtonian temperature feedback effects E.J. Allen , Modeling with Itô Stochastic Differential Equations , ( Springer , Dordrecht, The Netherlands 2007 ). Modeling with Itô Stochastic Differential Equations D.J Higham , An algorithmic introduction to numerical simulation of stochastic differential equations , SIAM Review 43 , 525 - 546 ( 2001 ). doi: 10.1137/S0036144500378302 http://doi.org/10.1137/S0036144500378302 An algorithmic introduction to numerical simulation of stochastic differential equations M.A. Akinlar , A. Secer , M. Bayram , Numerical Solution of Fractional Benney Equation , Appl. Math. Information Sci . 8 , 1633 - 1637 ( 2014 ). doi: 10.12785/amis/080418 http://doi.org/10.12785/amis/080418 Numerical Solution of Fractional Benney Equation M. Caputo , Elasticità Dissipazione , ( Zanichelli, Bologna, Italy , 1969 ). Elasticità Dissipazione I. Podlubny , Fractional Differential Equations , ( Academic Press , San Diego, CA 1999 ). Fractional Differential Equations S.S. Ray , Fractional Calculus with applications for Nuclear Reactor Dynamics , ( CRC Press , Taylor and Francis group , Boca Raton, New York, USA , 2016 ). Fractional Calculus with applications for Nuclear Reactor Dynamics M.M.R. Williams , Random Processes in Nuclear Reactors , ( Pergamon Press 1974 ). S. Kazem , Exact solution of some linear fractional differential equations by Laplace transform, Int . J. Nonlinear Sci . 16 , 3 - 11 ( 2013 ). Exact solution of some linear fractional differential equations by Laplace transform, Int A.A. Nahla , Analytical solution of the fractional point kinetics equations with multigroup of delayed neutrons during start-up of a nuclear reactor , Ann. Nucl. Energ . 99 , 247 - 252 ( 2017 ). doi: 10.1016/j.anucene.2016.08.030 http://doi.org/10.1016/j.anucene.2016.08.030 Analytical solution of the fractional point kinetics equations with multigroup of delayed neutrons during start-up of a nuclear reactor A.A. Nahla , A. Hemeda , Picard iteration and Padé approximations for stifffractional point kinetics equations , Appl. Math. Comput . 293 , 72 - 80 ( 2017 ). doi: 10.1016/j.amc.2016.08.008 http://doi.org/10.1016/j.amc.2016.08.008 Picard iteration and Padé approximations for stifffractional point kinetics equations A.E. Aboanber , Stiffness treatment of differential equations for the point reactor dynamic systems , Prog. Nucl. Energ . 71 , 248 - 257 ( 2014 ). doi: 10.1016/j.pnucene.2013.12.004 http://doi.org/10.1016/j.pnucene.2013.12.004 Stiffness treatment of differential equations for the point reactor dynamic systems A.A. Nahla , Analytical solution to solve the point reactor kinetics equations , Nucl. Eng. Des . 240 , 1622 - 1629 ( 2010 ). doi: 10.1016/j.nucengdes.2010.03.003 http://doi.org/10.1016/j.nucengdes.2010.03.003 Analytical solution to solve the point reactor kinetics equations J.A.M. Nobrega , A new solution of the point kinetics equations , Nucl. Sci. Eng . 46 , 366 - 375 ( 1971 ). A new solution of the point kinetics equations A.E. Aboanber , A.A. Nahla , Generalization of the analytical inversion method for the solution of the point kinetics equations , J. Phys. A: Math. Gen . 35 , 3245 - 3263 ( 2002 ). doi: 10.1088/0305-4470/35/14/307 http://doi.org/10.1088/0305-4470/35/14/307 Generalization of the analytical inversion method for the solution of the point kinetics equations A.E. Aboanber , A.A. Nahla , Solution of the point kinetics equations in the presence of Newtonian temperature feedback by Padé approximations via the analytical inversion method , J. Phys. A: Math. Gen . 35 , 9609 - 9627 ( 2002 ). doi: 10.1088/0305-4470/35/45/309 http://doi.org/10.1088/0305-4470/35/45/309 Solution of the point kinetics equations in the presence of Newtonian temperature feedback by Padé approximations via the analytical inversion method A.A. Nahla , Generalization of the analytical exponential model to solve the point kinetics equations of Be- and D2O-moderated reactors , Nucl. Eng. Des . 238 , 2648 - 2653 ( 2008 ). doi: 10.1016/j.nucengdes.2008.04.002 http://doi.org/10.1016/j.nucengdes.2008.04.002 Generalization of the analytical exponential model to solve the point kinetics equations of Be- and D2O-moderated reactors A.A. Nahla , An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects , Ann. Nucl. Energ . 38 , 2810 - 2817 ( 2011 ). doi: 10.1016/j.anucene.2011.08.021 http://doi.org/10.1016/j.anucene.2011.08.021 An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects A.A. Nahla , Numerical treatment for the point reactor kinetics equations using theta method, eigenvalues and eigenvectors , Prog. Nucl. Energ . 85 , 756 - 763 ( 2015 ). doi: 10.1016/j.pnucene.2015.09.008 http://doi.org/10.1016/j.pnucene.2015.09.008 Numerical treatment for the point reactor kinetics equations using theta method, eigenvalues and eigenvectors A. Patra , S.S. Ray , On the solution of the nonlinear fractional neutron point-kinetics equation with Newtonian temperature feedback reactivity , Nucl. Technol ., 189 ( 1 ), 103 - 109 ( 2015 ). doi: 10.13182/NT13-148 http://doi.org/10.13182/NT13-148 On the solution of the nonlinear fractional neutron point-kinetics equation with Newtonian temperature feedback reactivity 10.1007/s41365-018-0467-0 2017-09-16 2018-03-10 2018-03-29 2018-07-27 2018-09-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2018 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-018-0467-0 http://dx.doi.org/10.1007/s41365-018-0467-0 . 2018 Nuclear Science and Techniques 09 29