Higher order approximate solutions of fractional stochastic point kinetics equations in nuclear reactor dynamics

Higher order approximate solutions of fractional stochastic point kinetics equations in nuclear reactor dynamics 1 first-author 1 corresp Email: santanusaharay@yahoo.com Email: santanusaharay@yahoo.com santanusaharay@yahoo.com National Institute of Technology Rourkela Department of Mathematics Rourkela 769008, India National Institute of Technology Rourkela Department of Mathematics Rourkela 769008, India Stochastic point kinetics equations (SPKEs) are a system of Itô stochastic differential equations whose solution have been obtained by higher order approximation. In this study, a fractional model of SPKEs has been analyzed. The efficiency of the proposed higher order approximation scheme has been discussed in the results section. The solutions of SPKEs in the presence of Newtonian temperature feedback have also been provided to further discuss the physical behavior of the fractional model. Fractional stochastic point reactor kinetics equations Fractional calculus Higher order approximation Caputo derivative Neutron population. 1. Introduction 1 Stochastic point kinetics equations (SPKEs) are a system of coupled nonlinear stochastic differential equations (SDEs) [ 1 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 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 3 - 6 6 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 7 ]. Numerical solutions of SPKEs by implementing stochastic piecewise continuous approximation (PCA) have been obtained by Hayes and Allen [ 2 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 8 ]. Saha Ray and Patra applied the Grunwald–Letnikov definition of fractional derivative for solving SPKEs [ 9 9 ]. Nahla and Edress [ 10 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 1 We now introduce an elementary definition before stating the central limit theorem. Definition 2.1 Let <math id="M1"><mrow><msub><mover accent="true"><mi>X</mi><mo>¯</mo></mover><mi>M</mi></msub></mrow></math> be the sample mean of <math id="M2"><mi>M</mi></math> independent samples <math id="M3"><mrow><msub><mi>X</mi><mn>1</mn></msub><mo>,</mo><mn>...</mn><mo>,</mo><msub><mi>X</mi><mi>M</mi></msub></mrow></math> of a random variable X such that <math display="block" id="M4"><mrow><msub><mover accent="true"><mi>X</mi><mo>¯</mo></mover><mi>M</mi></msub><mo>=</mo><mfrac><mn>1</mn><mi>M</mi></mfrac><mo stretchy="false">(</mo><msub><mi>X</mi><mn>1</mn></msub><mo>+</mo><mn>...</mn><mo>+</mo><msub><mi>X</mi><mi>M</mi></msub><mo stretchy="false">)</mo><mo>,</mo><mtext> </mtext><msubsup><mi>X</mi><mi>M</mi><mo>*</mo></msubsup><mo>=</mo><msqrt><mi>M</mi></msqrt><mo stretchy="false">(</mo><msub><mover accent="true"><mi>X</mi><mo>¯</mo></mover><mi>M</mi></msub><mo>−</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where <math id="M5"><mi>μ</mi></math> is the mean and <math id="M6"><mrow><msup><mi>σ</mi><mn>2</mn></msup></mrow></math> is the variance for each independent and identically distributed real-valued random variables <math id="M7"><mrow><msub><mi>X</mi><mi>j</mi></msub></mrow></math> . Theorem 2.1 (Central Limit Theorem) If each random variable <math id="M8"><mrow><msub><mi>X</mi><mi>j</mi></msub></mrow></math> has a finite second moment with <math id="M9"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><msub><mi>X</mi><mi>j</mi></msub><mo>=</mo><msup><mi>σ</mi><mn>2</mn></msup></mrow></math> , the distribution of <math id="M10"><mrow><msub><mover accent="true"><mi>X</mi><mo>¯</mo></mover><mi>M</mi></msub></mrow></math> converges to that of a Gaussian random variable with mean <math id="M11"><mi>μ</mi></math> and variance <math id="M12"><mrow><mfrac><mrow><msup><mi>σ</mi><mn>2</mn></msup></mrow><mi>M</mi></mfrac></mrow></math> . Thus, <math id="M13"><mrow><msubsup><mi>X</mi><mi>M</mi><mo>*</mo></msubsup></mrow></math> converges in distribution to <math id="M14"><mrow><mi>Z</mi><mo>~</mo><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><msup><mi>σ</mi><mn>2</mn></msup><mo stretchy="false">)</mo></mrow></math> . The fractional Itô SDE for NPKEs with temperature feedback effects [ 11 11 - 15 15 11-15 ] obtained from the central limit theorem can be written as follows [ 16 16 ]: <math display="block" id="M15"><mrow><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mi>t</mi><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>A</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mi>Q</mi><mo>+</mo><msup><mi>B</mi><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mfrac><mrow><mi>d</mi><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M16"><mi>α</mi></math> is the fractional order of the derivative and <math id="M17"><mrow><mn>0</mn><mo><</mo><mi>α</mi><mo>≤</mo><mn>1</mn></mrow></math> , <math id="M18"><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</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>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow></mrow></math> , <math id="M19"><mrow><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>W</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>W</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>W</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>W</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow></mrow></math> , and <math id="M20"><mrow><mi>Q</mi><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></mtable></mrow><mo>)</mo></mrow></mrow></math> . Here, <math id="M21"><mi>m</mi></math> is the total number of delayed neutron groups and <math id="M22"><mrow><msub><mi>W</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo><msub><mi>W</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo><msub><mi>W</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>W</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> are standard Wiener processes as defined in [ 8 8 ], <math id="M23"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is the neutron population, and <math id="M24"><mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is the precursor concentration of <math id="M25"><mi>i</mi></math> -group of delayed neutrons. The coefficient matrix <math id="M26"><mrow><mi>A</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is represented in the following form <math id="M27"><mrow><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>W</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>W</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>W</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>W</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow></mrow></math> <math display="block" id="M28"><mrow><mi>A</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><mi>ρ</mi><mo>−</mo><mi>β</mi></mrow><mi>l</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>m</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>1</mn></msub></mrow><mi>l</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></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mn>2</mn></msub></mrow><mi>l</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></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><mfrac><mrow><msub><mi>β</mi><mi>m</mi></msub></mrow><mi>l</mi></mfrac></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>…</mo></mtd><mtd><mrow><mo>−</mo><msub><mi>λ</mi><mi>m</mi></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo><mo>,</mo></mrow></mrow></math> where <math id="M29"><mi>ρ</mi></math> is the total reactivity, <math id="M30"><mrow><mi>β</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow></mstyle></mrow></math> is the total fraction of delayed neutrons, <math id="M31"><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow></math> is the fraction and <math id="M32"><mrow><msub><mi>λ</mi><mi>i</mi></msub></mrow></math> is the decay constant of <math id="M33"><mi>i</mi></math> -group of delayed neutrons, and <math id="M34"><mi>l</mi></math> is the prompt neutron generation time. The covariance matrix <math id="M35"><mrow><mi>B</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> , which is evaluated in reference [ 16 16 ], is represented as follows: <math display="block" id="M36"><mrow><mi>B</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>μ</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mo>…</mo></mtd><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mn>1</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mo>…</mo></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mrow><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mo>…</mo></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋮</mo></mtd><mtd><mo>⋱</mo></mtd><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msub><mi>μ</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>…</mo></mtd><mtd><mrow><msub><mi>μ</mi><mi>m</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M37"><mrow><msub><mi>μ</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>(</mo><mrow><mfrac><mrow><mi>ρ</mi><mo>+</mo><mi>β</mi></mrow><mi>l</mi></mfrac></mrow><mo>)</mo></mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><mrow><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> , <math id="M38"><mrow><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mi>l</mi></mfrac><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>λ</mi><mi>i</mi></msub><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> , and <math id="M39"><mrow><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>3</mn><mo>,</mo><mo>…</mo><mo>,</mo><mi>m</mi></mrow></math> . 3. Higher order approximation scheme 1 3.1 Fractional calculus 2 The Caputo derivative operator for <math id="M40"><mrow><mi>α</mi><mo>∈</mo><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> is defined as follows: <math display="block" id="M41"><mrow><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mi>t</mi><mi>α</mi></msubsup><mi>f</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mn>0</mn><mi>t</mi></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><mi>t</mi><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><msup><mi>f</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mi>d</mi><mi>s</mi></mrow></mrow></mstyle><mo>,</mo></mrow></math> in which <math id="M42"><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mo>⋅</mo><mo stretchy="false">)</mo></mrow></math> is the Euler gamma function. In this paper, a <math id="M43"><mrow><msup><mrow><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow><mrow><mi>t</mi><mi>h</mi></mrow></msup></mrow></math> order scheme for Caputo derivative <math id="M44"><mrow><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mi>t</mi><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with <math id="M45"><mrow><mi>α</mi><mo>∈</mo><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> has been discussed for obtaining the solutions of fractional SNPKEs. Let <math id="M46"><mrow><mn>0</mn><mo>=</mo><msub><mi>t</mi><mn>0</mn></msub><mo><</mo><msub><mi>t</mi><mn>1</mn></msub><mo><</mo><mn>...</mn><mo><</mo><msub><mi>t</mi><mi>M</mi></msub><mo>=</mo><mi>T</mi></mrow></math> and <math id="M47"><mrow><msub><mi>t</mi><mi>n</mi></msub><mo>=</mo><msub><mi>t</mi><mn>0</mn></msub><mo>+</mo><mi>n</mi><mi>τ</mi></mrow></math> be the equidistant discretized times with <math id="M48"><mrow><mi>τ</mi><mo>=</mo><msub><mi>Δ</mi><mi>n</mi></msub><mo>=</mo><mfrac><mrow><mo stretchy="false">(</mo><mi>T</mi><mo>−</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo></mrow><mi>M</mi></mfrac></mrow></math> for <math id="M49"><mrow><mi>M</mi><mo>∈</mo><mi>Z</mi></mrow></math> such that <math id="M50"><mrow><mi>τ</mi><mo>∈</mo><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> . Now using the Taylor expansion to <math id="M51"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo></mrow></math> , <math id="M52"><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></math> , and <math id="M53"><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></math> at point <math id="M54"><mrow><mi>t</mi><mo>=</mo><msub><mi>t</mi><mi>i</mi></msub></mrow></math> <math id="M55"><mrow><mo stretchy="false">(</mo><mn>0</mn><mo>≤</mo><mi>i</mi><mo><</mo><mi>n</mi><mo stretchy="false">)</mo></mrow></math> , we get <math display="block" id="M56"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>=</mo><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><msup><mi>Ψ</mi><mo>″</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mo>!</mo></mrow></mfrac><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></math> , <math display="block" id="M57"><mrow><mi>s</mi><mo>∈</mo><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo>,</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>,</mo></mrow></math> <math id="M58"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mi>τ</mi></mrow></mfrac><mo>−</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>!</mo></mrow></mfrac><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>4</mn></msup><mo stretchy="false">)</mo></mrow></math> , and <math display="block" id="M59"><mrow><msup><mi>Ψ</mi><mo>″</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><msup><mi>τ</mi><mn>2</mn></msup></mrow></mfrac><mo>−</mo><mfrac><mrow><msup><mi>Ψ</mi><mrow><mo stretchy="false">(</mo><mn>4</mn><mo stretchy="false">)</mo></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>12</mn></mrow></mfrac><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>4</mn></msup><mo stretchy="false">)</mo><mo>.</mo></mrow></math> Hence, we can obtain the following: <math id="M60"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mi>τ</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><msup><mi>τ</mi><mn>2</mn></msup></mrow></mfrac><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>−</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>!</mo></mrow></mfrac><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mo>!</mo></mrow></mfrac><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></math> , <math id="M61"><mrow><mn>0</mn><mo><</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo><</mo><mi>τ</mi></mrow></math> . Therefore, the Caputo derivative can be discretized as <math id="M62"><mtable columnalign="left"><mtr><mtd><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mi>t</mi><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><msub><mo>|</mo><mrow><mi>t</mi><mo>=</mo><msub><mi>t</mi><mi>n</mi></msub></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mn>0</mn><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mi>d</mi><mi>s</mi></mrow></mrow></mstyle></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo><mi>d</mi><mi>s</mi></mrow></mrow></mstyle></mrow></mstyle></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><mrow><mo>[</mo><mrow><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mi>τ</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>−</mo><mn>2</mn><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>+</mo><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow><mrow><msup><mi>τ</mi><mn>2</mn></msup></mrow></mfrac><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mrow></mrow></mrow></mstyle></mrow></mstyle></mtd></mtr><mtr><mtd><mrow><mrow><mo>−</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>!</mo></mrow></mfrac><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mo>!</mo></mrow></mfrac><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup></mrow><mo>]</mo></mrow><mi>d</mi><mi>s</mi><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">[</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mstyle><mo>+</mo><msup><mi>r</mi><mi>n</mi></msup><mo>.</mo></mtd></mtr></mtable></math> where <math id="M63"><mrow><mi>n</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mrow></math> and <math id="M64"><mrow><msup><mi>τ</mi><mi>n</mi></msup></mrow></math> is the truncation error. Therefore, the Caputo derivative has the following numerical approximation [ 17 17 ]: <math display="block" id="M65"><mrow><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">[</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mstyle><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">)</mo><mo>,</mo></mrow></math> where <math id="M66"><mrow><mn>0</mn><mo><</mo><mi>α</mi><mo><</mo><mn>1</mn></mrow></math> , <math id="M67"><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo>=</mo><mfrac><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow><mn>2</mn></mfrac><mo stretchy="false">[</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>i</mi><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>i</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">]</mo></mrow></math> , <math id="M68"><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo>=</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>i</mi><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>i</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mi>i</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup></mrow></math> , and <math id="M69"><mrow><msup><mi>r</mi><mi>n</mi></msup></mrow></math> is the truncation error in the following form: <math display="block" id="M70"><mrow><msup><mi>r</mi><mi>n</mi></msup><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">[</mo><mo>−</mo><msub><mi>C</mi><mi>Ψ</mi></msub><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mn>3</mn><msub><mi>C</mi><mi>Ψ</mi></msub><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo stretchy="false">]</mo><mi>d</mi><mi>s</mi><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mrow></mstyle></mrow></mstyle><mo>,</mo></mrow></math> where <math id="M71"><mrow><msub><mi>C</mi><mi>Ψ</mi></msub><mo>=</mo><mfrac><mrow><msup><mi>Ψ</mi><mo>‴</mo></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>!</mo></mrow></mfrac></mrow></math> is a constant. The right-hand side of eq. (3.2) can be expressed as follows: <math display="block" id="M72"><mrow><mfrac><mn>1</mn><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow><mrow><msub><mi>t</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></munderover><mrow><msup><mrow><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo>−</mo><mi>s</mi><mo stretchy="false">)</mo></mrow><mrow><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">[</mo><mo>−</mo><msub><mi>C</mi><mi>Ψ</mi></msub><msup><mi>τ</mi><mn>2</mn></msup><mo>+</mo><mn>3</mn><msub><mi>C</mi><mi>Ψ</mi></msub><msup><mrow><mo stretchy="false">(</mo><mi>s</mi><mo>−</mo><msub><mi>t</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mo stretchy="false">]</mo><mi>d</mi><mi>s</mi><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></mrow></mstyle></mrow></mstyle><mo>=</mo><mfrac><mrow><msub><mi>C</mi><mi>Ψ</mi></msub></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">(</mo><msub><mi>I</mi><mn>1</mn></msub><mo>+</mo><mn>.3</mn><msub><mi>I</mi><mn>2</mn></msub><mo stretchy="false">)</mo></mrow></mstyle><mo>,</mo></mrow></math> where http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10313427&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10313421&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10313427&type= Therefore, <math display="block" id="M73"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><msub><mi>C</mi><mi>Ψ</mi></msub></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">(</mo><msub><mi>I</mi><mn>1</mn></msub><mo>+</mo><mn>.3</mn><msub><mi>I</mi><mn>2</mn></msub><mo stretchy="false">)</mo></mrow></mstyle><mo>=</mo><mfrac><mrow><msub><mi>C</mi><mi>Ψ</mi></msub><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mrow><mo>{</mo><mrow><mo>−</mo><msup><mi>n</mi><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><mn>3</mn><mo stretchy="false">[</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>+</mo><mn>...</mn><mo>+</mo><msup><mn>2</mn><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>+</mo><mn>1</mn><mo stretchy="false">]</mo><mo>−</mo><mfrac><mn>6</mn><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></mfrac><mo stretchy="false">[</mo><msup><mrow><mo stretchy="false">(</mo><mi>n</mi><mo>−</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>+</mo><mn>...</mn><mo>+</mo><msup><mn>2</mn><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>+</mo><mn>1</mn><mo stretchy="false">]</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mrow><mo>+</mo><mfrac><mn>6</mn><mrow><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><msup><mi>n</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup></mrow><mo>}</mo></mrow><mo>.</mo></mtd></mtr></mtable></math> Let <math display="block" id="M74"><mtable columnalign="left"><mtr><mtd><mi>S</mi><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>=</mo><mo>−</mo><mn>3</mn><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mi>i</mi><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup></mrow></mstyle><mo>−</mo><mfrac><mn>6</mn><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msup><mi>i</mi><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></msup></mrow></mstyle><mo>+</mo><mfrac><mn>6</mn><mrow><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><msup><mi>n</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><msup><mi>n</mi><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>a</mi><mi>i</mi></msub></mrow></mstyle><mo>,</mo><mtext> </mtext><mi>n</mi><mo>≥</mo><mn>1.</mn></mtd></mtr></mtable></math> If <math id="M75"><mrow><mi>n</mi><mo>=</mo><mn>1</mn></mrow></math> , define <math id="M76"><mrow><msub><mi>a</mi><mn>0</mn></msub><mo>=</mo><mi>s</mi><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>6</mn><mrow><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mo>−</mo><mn>1</mn></mrow></math> . Then <math id="M77"><mrow><msub><mi>a</mi><mi>i</mi></msub></mrow></math> ( <math id="M78"><mrow><mi>i</mi><mo>≤</mo><mn>1</mn></mrow></math> ) can be defined as follows: <math display="block" id="M79"><mrow><msub><mi>a</mi><mi>i</mi></msub><mo>=</mo><mi>S</mi><mo stretchy="false">(</mo><mi>i</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo><mo>−</mo><mi>S</mi><mo stretchy="false">(</mo><mi>i</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>6</mn><mrow><mo stretchy="false">(</mo><mn>2</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mo stretchy="false">[</mo><msup><mrow><mo stretchy="false">(</mo><mi>i</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><msup><mi>i</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">]</mo><mo>−</mo><mfrac><mn>6</mn><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></mfrac><msup><mi>i</mi><mrow><mn>2</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><msup><mrow><mo stretchy="false">(</mo><mi>i</mi><mo>+</mo><mn>1</mn><mo stretchy="false">)</mo></mrow><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>−</mo><mn>2</mn><msup><mi>i</mi><mrow><mn>1</mn><mo>−</mo><mi>α</mi></mrow></msup><mo>.</mo></mrow></math> It can be proven that <math id="M80"><mrow><mo>|</mo><mi>S</mi><mo stretchy="false">(</mo><mi>n</mi><mo stretchy="false">)</mo><mo>|</mo></mrow></math> is bounded for <math id="M81"><mrow><mi>n</mi><mo>≥</mo><mn>1</mn></mrow></math> [ 17 17 - 19 19 17-19 ]. This proves that the series <math id="M82"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover><mrow><msub><mi>a</mi><mi>i</mi></msub></mrow></mstyle></mrow></math> converges. On further simplification we get, <math display="block" id="M83"><mtable columnalign="left"><mtr><mtd><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mo stretchy="false">[</mo><msub><mi>Ψ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle></mtd></mtr><mtr><mtd><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></mstyle><mo stretchy="false">]</mo><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">)</mo></mtd></mtr></mtable></math> In Eq. (3.2), if <math id="M84"><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow></math> , then <math id="M85"><mrow><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub></mrow></math> , which lies outside of <math id="M86"><mrow><mo stretchy="false">[</mo><mn>0</mn><mo>,</mo><mi>T</mi><mo stretchy="false">]</mo></mrow></math> . Various options have been provided to approach <math id="M87"><mrow><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub></mrow></math> . In numerical calculations, the neighboring function values have been used to approximate <math id="M88"><mrow><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub></mrow></math> , that is, <math id="M89"><mrow><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><mi>Ψ</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>−</mo><mi>τ</mi><msup><mi>Ψ</mi><mo>′</mo></msup><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>+</mo><mfrac><mrow><msup><mi>τ</mi><mn>2</mn></msup></mrow><mn>2</mn></mfrac><msup><mi>Ψ</mi><mo>″</mo></msup><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></math> . 1) When <math id="M90"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mtext>(</mtext><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><msup><mi>Ψ</mi><mo>″</mo></msup><mtext>(</mtext><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn></mrow></math> , then <math id="M91"><mrow><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>Ψ</mi><mn>0</mn></msub><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></math> ; the convergence order is <math id="M92"><mrow><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">)</mo></mrow></math> . 2) When <math id="M93"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mtext>(</mtext><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0</mn></mrow></math> , <math id="M94"><mrow><msup><mi>Ψ</mi><mo>″</mo></msup><mtext>(</mtext><mn>0</mn><mo stretchy="false">)</mo><mo>≠</mo><mn>0</mn></mrow></math> , then <math id="M95"><mrow><msub><mi>Ψ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msub><mo>=</mo><msub><mi>Ψ</mi><mn>0</mn></msub><mo>+</mo><mfrac><mrow><msup><mi>τ</mi><mn>2</mn></msup></mrow><mn>2</mn></mfrac><msup><mi>Ψ</mi><mo>″</mo></msup><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>+</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>3</mn></msup><mo stretchy="false">)</mo></mrow></math> ; the convergence order of Eq. (3.1) is <math id="M96"><mrow><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mn>2</mn></msup><mo stretchy="false">)</mo></mrow></math> . 3) When <math id="M97"><mrow><msup><mi>Ψ</mi><mo>′</mo></msup><mtext>(</mtext><mn>0</mn><mo stretchy="false">)</mo><mo>≠</mo><mn>0</mn></mrow></math> , then the convergence order is <math id="M98"><mrow><mi>O</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo></mrow></math> . 4. Solution of SPKE by higher order approximation method 1 In this section, higher order approximation to Caputo derivative has been applied to Eq. (2.1) as follow: <math display="block" id="M99"><mtable columnalign="left"><mtr><mtd><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mo stretchy="false">[</mo><msub><mi>Ψ</mi><mi>n</mi></msub><mo stretchy="false">(</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></mstyle><mo stretchy="false">]</mo></mtd></mtr><mtr><mtd><mo>=</mo><mi>A</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mi>Ψ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mi>Q</mi><mo>+</mo><msup><mi>B</mi><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mfrac><mrow><mi>d</mi><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mi>d</mi><mi>t</mi></mrow></mfrac></mtd></mtr></mtable></math> Therefore, the above expression can be simplified into an explicit numerical scheme as follows: <math display="block" id="M100"><mtable columnalign="left"><mtr><mtd><msub><mi>Ψ</mi><mi>n</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mo stretchy="false">(</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mrow></mfrac><mrow><mo>(</mo><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo><msup><mi>τ</mi><mi>α</mi></msup><mrow><mo>(</mo><mrow><mi>A</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><msub><mi>Ψ</mi><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>+</mo><mi>Q</mi><mo>+</mo><msup><mi>B</mi><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mi>Δ</mi><mi>W</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mrow><mo>−</mo><mrow><mo>(</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>−</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>2</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mstyle><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mtd></mtr></mtable></math> where <math id="M101"><mrow><mi>Δ</mi><mi>W</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>=</mo><msqrt><mi>h</mi></msqrt><msub><mi>S</mi><mi>n</mi></msub></mrow></math> and <math id="M102"><mrow><mi>n</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mn>...</mn><mo>,</mo><mi>M</mi></mrow></math> with initial condition <math display="block" id="M103"><mrow><mi>Ψ</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</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>l</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>l</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>l</mi><msub><mi>λ</mi><mi>m</mi></msub></mrow></mfrac></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M104"><mrow><msub><mi>S</mi><mi>n</mi></msub><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>S</mi><mrow><msub><mn>0</mn><mi>n</mi></msub></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>S</mi><mrow><msub><mn>1</mn><mi>n</mi></msub></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>S</mi><mrow><msub><mn>2</mn><mi>n</mi></msub></mrow></msub></mrow></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mrow><msub><mi>S</mi><mrow><msub><mi>m</mi><mi>n</mi></msub></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M105"><mrow><msub><mi>S</mi><mrow><msub><mn>0</mn><mi>n</mi></msub></mrow></msub><mo>,</mo><msub><mi>S</mi><mrow><msub><mn>1</mn><mi>n</mi></msub></mrow></msub><mo>,</mo><msub><mi>S</mi><mrow><msub><mn>2</mn><mi>n</mi></msub></mrow></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mi>S</mi><mrow><msub><mi>m</mi><mi>n</mi></msub></mrow></msub></mrow></math> are random variables from <math id="M106"><mrow><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo>,</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> with mean equal to 0 and variance equal to 1. Theorem 4.1 The local truncation error of the scheme is <math id="M107"><mrow><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">)</mo></mrow></math> . Proof: The local truncation error of the higher order scheme for Eq. (4.2) can be derived as follows: <math display="block" id="M108"><mtable columnalign="left"><mtr><mtd><msubsup><mi>R</mi><mi>j</mi><mi>k</mi></msubsup><mo>=</mo><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">[</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mstyle><mo>−</mo><mi>A</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><msub><mi>Ψ</mi><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>−</mo><mi>Q</mi><mo>−</mo><msup><mi>B</mi><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mi>Δ</mi><mi>W</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><msup><mi>τ</mi><mrow><mo>−</mo><mi>α</mi></mrow></msup></mrow><mrow><mi>Γ</mi><mo stretchy="false">(</mo><mn>3</mn><mo>−</mo><mi>α</mi><mo stretchy="false">)</mo></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>0</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><mo stretchy="false">[</mo><msub><mi>w</mi><mrow><mn>1</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>+</mo><msub><mi>w</mi><mrow><mn>2</mn><mo>,</mo><mi>n</mi><mo>−</mo><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>+</mo><mn>1</mn></mrow></msub><mo>−</mo><mn>2</mn><msub><mi>Ψ</mi><mi>i</mi></msub><mo>+</mo><msub><mi>Ψ</mi><mrow><mi>i</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">]</mo></mrow></mstyle><mo>−</mo><mmultiscripts><mi>D</mi><mprescripts/><mn>0</mn><mi>C</mi></mmultiscripts><msubsup><mrow/><mrow><msub><mi>t</mi><mi>n</mi></msub></mrow><mi>α</mi></msubsup><mi>Ψ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>n</mi></msub><mo stretchy="false">)</mo><mo>−</mo><mi>A</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo stretchy="false">(</mo><msub><mi>Ψ</mi><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></msub><mo>−</mo><msub><mi>Ψ</mi><mrow><mi>n</mi><mo>-</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>=</mo><mi>O</mi><mo stretchy="false">(</mo><msup><mi>τ</mi><mrow><mn>3</mn><mo>−</mo><mi>α</mi></mrow></msup><mo stretchy="false">)</mo><mo>.</mo></mtd></mtr></mtable></math> 5. Numerical results and discussions 1 Here the solutions of SPKE ( <math id="M109"><mrow><mi>i</mi><mo>=</mo><mn>6</mn></mrow></math> ) have been obtained using higher order approximation scheme. 5.1 Step reactivity 2 In this section, the numerical solutions of the fractional stochastic point kinetic model [ 10 10 ] using the following parameters have been obtained: <math id="M110"><mrow><msub><mi>λ</mi><mn>1</mn></msub><mo>=</mo><mn>0.0127</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M111"><mrow><msub><mi>λ</mi><mn>2</mn></msub><mo>=</mo><mn>0.0317</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M112"><mrow><msub><mi>λ</mi><mn>3</mn></msub><mo>=</mo><mn>0.115</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M113"><mrow><msub><mi>λ</mi><mn>4</mn></msub><mo>=</mo><mn>0.311</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M114"><mrow><msub><mi>λ</mi><mn>5</mn></msub><mo>=</mo><mn>1.4</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M115"><mrow><msub><mi>λ</mi><mn>6</mn></msub><mo>=</mo><mn>3.87</mn><mtext> </mtext><msup><mtext>s</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M116"><mrow><msub><mi>β</mi><mn>1</mn></msub><mo>=</mo><mn>0.000266</mn></mrow></math> , <math id="M117"><mrow><msub><mi>β</mi><mn>2</mn></msub><mo>=</mo><mn>0.001491</mn></mrow></math> , <math id="M118"><mrow><msub><mi>β</mi><mn>3</mn></msub><mo>=</mo><mn>0.001316</mn></mrow></math> , <math id="M119"><mrow><msub><mi>β</mi><mn>4</mn></msub><mo>=</mo><mn>0.002849</mn></mrow></math> , <math id="M120"><mrow><msub><mi>β</mi><mn>5</mn></msub><mo>=</mo><mn>0.000896</mn></mrow></math> , <math id="M121"><mrow><msub><mi>β</mi><mn>6</mn></msub><mo>=</mo><mn>0.000182</mn></mrow></math> , <math id="M122"><mrow><mi>β</mi><mo>=</mo><mn>0.007</mn></mrow></math> , <math id="M123"><mrow><mi>l</mi><mo>=</mo><mn>2.0</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><mtext> </mtext><mtext>s</mtext></mrow></math> , <math id="M124"><mrow><mi>v</mi><mo>=</mo><mn>2.5</mn></mrow></math> , and <math id="M125"><mrow><mi>q</mi><mo>=</mo><mn>0</mn></mrow></math> with <math id="M126"><mrow><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><msub><mi>N</mi><mn>0</mn></msub><mo>=</mo><mn>100</mn></mrow></math> (neutron) and <math id="M127"><mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow><mrow><mi>l</mi><msub><mi>λ</mi><mi>i</mi></msub></mrow></mfrac></mrow></math> . For different step reactivities, i.e., <math id="M128"><mrow><mi>ρ</mi><mo>=</mo><mn>0.003</mn></mrow></math> and <math id="M129"><mrow><mi>ρ</mi><mo>=</mo><mn>0.007</mn></mrow></math> , the mean peak values of <math id="M130"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with respect to time for fractional orders <math id="M131"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , 0.98, and 0.99 at step size <math id="M132"><mrow><mi>h</mi><mo>=</mo><mn>0.001</mn><mtext> </mtext><mtext>s</mtext></mrow></math> and for 500 trials are presented in Table 1 Table 1 . In addition, the results have been compared with those obtained by other methods, namely EM [ 8 8 ], Taylor 1.5 strong order [ 8 8 ], AEM [ 20 20 ], and efficient stochastic model (ESM) [ 10 10 ] with graphical representation in Fig. 1 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. 1 1 ρ 2 1 α = 0.96 2 1 α = 0.98 2 1 α = 0.99 1 1 EM α = 1[8] 1 1 Taylor 1.5 strong order α = 1 [8] 1 1 AEM α = 1[20] 1 1 ESM α = 1[10] 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Peak 1 1 Peak 1 1 Peak 1 1 0.003 1 1 180.796 1 1 0.1 1 1 180.033 1 1 0.095 1 1 180.819 1 1 0.083 1 1 208.6 1 1 199.408 1 1 186.30 1 1 179.93 1 1 0.007 1 1 128.655 1 1 0.001 1 1 128.248 1 1 0.001 1 1 124.113 1 1 0.001 1 1 139.568 1 1 139.569 1 1 134.54 1 1 134.96 5.2 Ramp reactivity 2 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 <math id="M141"><mrow><mi>ρ</mi><mo>=</mo><mn>0.1</mn><mi>β</mi><mi>t</mi></mrow></math> . For ramp external reactivity, i.e., <math id="M142"><mrow><mi>ρ</mi><mo>=</mo><mn>0.1</mn><mi>β</mi><mi>t</mi></mrow></math> , the mean peak values of <math id="M143"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with respect to time for fractional orders <math id="M144"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , 0.98, and 0.99 at step size <math id="M145"><mrow><mi>h</mi><mo>=</mo><mn>0.001</mn><mtext> </mtext><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo></mrow></math> and for 500 trials are listed in Table 2 Table 2 . In addition, the results have been compared with those obtained using other methods, namely AEM [ 20 20 ] and ESM [ 10 10 ] with graphical representation in Fig. 2 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. 1 1 α 1 1 σ 2 1 α = 0.96 2 1 α = 0.98 2 1 α = 0.99 1 1 AEM α = 1 [20] 1 1 ESM α = 1 [10] 1 1 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Peak 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Peak 1 1 0.1 βt 1 1 10 -11 1 1 113.563 1 1 0.998 1 1 113.275 1 1 186.30 1 1 113.045 1 1 0.1 1 1 113.267707 1 1 113.116433 5.3 Temperature feedback reactivity 2 In this section, solutions of fractional stochastic point kinetic model with <math id="M156"><mi>i</mi></math> -group of delayed neutrons ( <math id="M157"><mrow><mi>i</mi><mo>=</mo><mn>6</mn></mrow></math> ) 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 21 ] is represented in the following form: <math display="block" id="M158"><mrow><mi>ρ</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>ρ</mi><mrow><mi>e</mi><mi>x</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><msub><mi>ρ</mi><mi>f</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>,</mo><mo>　</mo><msub><mi>ρ</mi><mi>f</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>σ</mi><mstyle displaystyle="true"><mrow><munderover><mo>∫</mo><mn>0</mn><mi>t</mi></munderover><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">)</mo><mi>d</mi><mi>τ</mi></mrow></mrow></mstyle></mrow></math> where <math id="M159"><mrow><mi>σ</mi><mo>=</mo><mi>α</mi><msub><mi>K</mi><mi>c</mi></msub></mrow></math> , <math id="M160"><mrow><msub><mi>ρ</mi><mrow><mi>e</mi><mi>x</mi></mrow></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is the external reactivity, <math id="M161"><mrow><mi>T</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> is the temperature, <math id="M162"><mrow><msub><mi>T</mi><mn>0</mn></msub></mrow></math> is the initial temperature, <math id="M163"><mi>α</mi></math> is the temperature coefficient, and <math id="M164"><mrow><msub><mi>K</mi><mtext>c</mtext></msub></mrow></math> is the reciprocal of the thermal capacity. In the interval <math id="M165"><mrow><mo stretchy="false">[</mo><msub><mi>t</mi><mi>k</mi></msub><mo>,</mo><msub><mi>t</mi><mrow><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">]</mo></mrow></math> , the total reactivity can be expressed as follows [ 15 15 ]: <math display="block" id="M166"><mrow><mi>ρ</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>k</mi></msub><mo stretchy="false">)</mo><mo>≈</mo><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>t</mi><mi>k</mi></msub><mo stretchy="false">)</mo><mo>−</mo><mi>h</mi><mi>σ</mi><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>k</mi></munderover><mrow><mi>N</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>j</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><mo>.</mo></mrow></math> 5.3.1 Step external reactivity 3 In this section, the numerical solutions of the SPKE of <math id="M167"><mrow><msup><mi>U</mi><mrow><mn>235</mn></mrow></msup></mrow></math> nuclear reactor [ 3 3 ] using the following parameters have been obtained: <math id="M168"><mrow><msub><mi>λ</mi><mn>1</mn></msub><mo>=</mo><mn>0.0124</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M169"><mrow><msub><mi>λ</mi><mn>2</mn></msub><mo>=</mo><mn>0.0305</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M170"><mrow><msub><mi>λ</mi><mn>3</mn></msub><mo>=</mo><mn>0.111</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M171"><mrow><msub><mi>λ</mi><mn>4</mn></msub><mo>=</mo><mn>0.301</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M172"><mrow><msub><mi>λ</mi><mn>5</mn></msub><mo>=</mo><mn>1.13</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M173"><mrow><msub><mi>λ</mi><mn>6</mn></msub><mo>=</mo><mn>3.0</mn><mtext> </mtext><msup><mrow><mtext> s</mtext></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , <math id="M174"><mrow><msub><mi>β</mi><mn>1</mn></msub><mo>=</mo><mn>0.00021</mn></mrow></math> , <math id="M175"><mrow><msub><mi>β</mi><mn>2</mn></msub><mo>=</mo><mn>0.00141</mn></mrow></math> , <math id="M176"><mrow><msub><mi>β</mi><mn>3</mn></msub><mo>=</mo><mn>0.00127</mn></mrow></math> , <math id="M177"><mrow><msub><mi>β</mi><mn>4</mn></msub><mo>=</mo><mn>0.00255</mn></mrow></math> , <math id="M178"><mrow><msub><mi>β</mi><mn>5</mn></msub><mo>=</mo><mn>0.00074</mn></mrow></math> , <math id="M179"><mrow><msub><mi>β</mi><mn>6</mn></msub><mo>=</mo><mn>0.00027</mn></mrow></math> , <math id="M180"><mrow><mi>β</mi><mo>=</mo><mn>0.00645</mn></mrow></math> , <math id="M181"><mrow><mi>l</mi><mo>=</mo><mn>5.0</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><mtext> </mtext><mtext>s</mtext></mrow></math> , <math id="M182"><mrow><mi>α</mi><mo>=</mo><mn>5.0</mn><mo>×</mo><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>5</mn></mrow></msup><mtext> </mtext><msup><mtext>K</mtext><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> , and <math id="M183"><mrow><msub><mi>K</mi><mtext>c</mtext></msub><mo>=</mo><mn>0.05</mn><mtext> </mtext><msub><mrow><mtext> K/MW</mtext></mrow><mtext>s</mtext></msub></mrow></math> with <math id="M184"><mrow><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><msub><mi>N</mi><mn>0</mn></msub><mo>=</mo><mn>1</mn></mrow></math> (neutron) and <math id="M185"><mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo></mrow><mrow><mi>l</mi><msub><mi>λ</mi><mi>i</mi></msub></mrow></mfrac></mrow></math> . For different step external reactivities, i.e., <math id="M186"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.5</mn><mi>β</mi></mrow></math> , <math id="M187"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.75</mn><mi>β</mi></mrow></math> , and <math id="M188"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mi>β</mi></mrow></math> , the mean peak values of <math id="M189"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> for fractional orders <math id="M190"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , 0.98, and 0.99 and for different step external reactivities <math id="M191"><mrow><mn>0.5</mn><mi>β</mi></mrow></math> , <math id="M192"><mrow><mn>0.75</mn><mi>β</mi></mrow></math> , and <math id="M193"><mi>β</mi></math> using 500 trials are presented in Table 3(a) 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 15 ] in Table 3(b) Table 3(b) with graphical representation in Figs. 3 Figs. 3 , 4 4 , and 5 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. 1 1 ρ ex 2 1 α = 0.96 2 1 α = 0.98 2 1 α = 0.99 1 1 ρ ex 2 1 SSFEMM [15] (α = 1) 2 1 DFMM [15] (α = 1) 2 1 α = 0.98 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 (a) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 (b) 1 1 0.5β 1 1 42.6182 1 1 28.65 1 1 44.789 1 1 30.29 1 1 45.9708 1 1 29.25 1 1 0.5β 1 1 46.4939 1 1 28.34 1 1 46.2606 1 1 27.84 1 1 44.789 1 1 30.29 1 1 0.75β 1 1 159.21 1 1 8.875 1 1 160.124 1 1 8.895 1 1 162.99 1 1 9.305 1 1 0.75β 1 1 163.707 1 1 8.795 1 1 164.22 1 1 8.95 1 1 160.124 1 1 8.895 1 1 β 1 1 801.166 1 1 0.985 1 1 795.268 1 1 1.03 1 1 772.893 1 1 1.0625 1 1 β 1 1 760.589 1 1 1.065 1 1 769.238 1 1 1.0575 1 1 795.268 1 1 1.03 For <math id="M214"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , <math id="M215"><mrow><mn>0.98</mn></mrow></math> , and <math id="M216"><mrow><mn>0.99</mn></mrow></math> , the mean <math id="M217"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with two arbitrary sample paths have been shown for <math id="M218"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.5</mn><mi>β</mi></mrow></math> , <math id="M219"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.75</mn><mi>β</mi></mrow></math> , and <math id="M220"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mi>β</mi></mrow></math> . 5.3.2 Ramp external reactivity 3 The numerical solutions of the fractional SNPKE of <math id="M221"><mrow><msup><mi>U</mi><mrow><mn>235</mn></mrow></msup></mrow></math> nuclear reactor have been obtained using parameters similar to those in subsection 5.3.1. Here, the external reactivity is represented as <math id="M222"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.1</mn><mi>t</mi></mrow></math> and <math id="M223"><mrow><msub><mi>ρ</mi><mrow><mtext>ex</mtext></mrow></msub><mo>=</mo><mn>0.01</mn><mi>t</mi></mrow></math> , the nonlinear coefficient <math id="M224"><mi>σ</mi></math> takes <math id="M225"><mrow><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>11</mn></mrow></msup></mrow></math> or <math id="M226"><mrow><msup><mrow><mn>10</mn></mrow><mrow><mo>−</mo><mn>13</mn></mrow></msup></mrow></math> . For different ramp external reactivities, i.e., <math id="M227"><mrow><msub><mi>ρ</mi><mrow><mi>e</mi><mi>x</mi></mrow></msub><mo>=</mo><mn>0.1</mn><mi>t</mi></mrow></math> and <math id="M228"><mrow><mn>0.01</mn><mi>t</mi></mrow></math> , the mean peak values of <math id="M229"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with respect to time for fractional orders <math id="M230"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , 0.98, and 0.99 at step size <math id="M231"><mrow><mi>h</mi><mo>=</mo><mn>0.001</mn><mtext> </mtext><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo></mrow></math> and for 500 trials are listed in Table 4(a) Table 4(a) . These results have been compared with previously obtained results [ 10 10 , 15 15 ] in Table 4(b) Table 4(b) with graphical representation in Figs. 6(a)–(c) Figs. 6(a)–(c) , 7(a)–(c) 7(a)–(c) , and 8(a)–(c) 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. 1 1 α 1 1 σ 2 1 α = 0.96 2 1 α = 0.98 2 1 α = 0.99 1 1 α 1 1 σ 2 1 SSFEMM [15] (α = 1) 2 1 DFMM [15] (α = 1) 2 1 AEM [10] (α = 1) 2 1 α = 0.98 1 1 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 (a) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 (b) 1 1 0.01 1 1 10 -11 1 1 1.69869 E + 10 1 1 1.084 1 1 1.76031 E + 10 1 1 1.103 1 1 1.73211 E + 10 1 1 1.112 1 1 0.01 1 1 10 -11 1 1 1.68604 E + 10 1 1 1.118 1 1 1.69492 E + 10 1 1 1.118 1 1 1.673436 E + 10 1 1 0.854 1 1 1.76031 E + 10 1 1 1.103 1 1 1 1 10 -13 1 1 2.0663 E + 12 1 1 1.132 1 1 2.18262 E + 12 1 1 1.151 1 1 2.13422 E + 12 1 1 1.161 1 1 1 1 10 -13 1 1 2.12034 E + 12 1 1 1.169 1 1 2.12802 E + 12 1 1 1.168 1 1 2.082531 E + 12 1 1 0.877 1 1 2.18262 E + 12 1 1 1.151 1 1 0.1 1 1 10 -11 1 1 1.78236 E + 11 1 1 0.223 1 1 1.90873 E + 11 1 1 0.232 1 1 2.01965 E + 11 1 1 0.235 1 1 0.1 1 1 10 -11 1 1 1.88849 E + 11 1 1 0.235 1 1 1.89642 E + 11 1 1 0.235 1 1 1.790577 E + 11 1 1 0.142 1 1 1.90873 E + 11 1 1 0.232 1 1 1 1 10 -13 1 1 2.34211 E + 13 1 1 0.238 1 1 2.37627 E + 13 1 1 0.248 1 1 2.41795 E + 13 1 1 0.252 1 1 1 1 10 -13 1 1 2.24448 E + 13 1 1 0.251 1 1 2.26026 E + 13 1 1 0.251 1 1 2.143778 E + 13 1 1 0.150 1 1 2.37627 E + 13 1 1 0.248 5.4 Sinusoidal reactivity 2 In this section, the reactivity is in the form of sinusoidal change, i.e., <math id="M280"><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> . The numerical solution for this reactivity has been obtained using the following parameters: <math id="M281"><mrow><msub><mi>ρ</mi><mn>0</mn></msub><mo>=</mo><mn>0.005333</mn></mrow></math> , <math id="M282"><mrow><msub><mi>β</mi><mn>1</mn></msub><mo>=</mo><mi>β</mi><mo>=</mo><mn>0.0079</mn></mrow></math> , <math id="M283"><mrow><msub><mi>λ</mi><mn>1</mn></msub><mo>=</mo><mn>0.077</mn></mrow></math> , <math id="M284"><mrow><mi>Λ</mi><mo>=</mo><mn>0.001</mn></mrow></math> , <math id="M285"><mrow><mi>q</mi><mo>=</mo><mn>0</mn></mrow></math> , <math id="M286"><mrow><mi>N</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>=</mo><msub><mi>N</mi><mn>0</mn></msub><mo>=</mo><mn>1</mn></mrow></math> , and time period <math id="M287"><mrow><mi>T</mi><mo>=</mo><mn>100</mn><mtext> </mtext><mtext> </mtext><mo stretchy="false">(</mo><mi>s</mi><mo stretchy="false">)</mo></mrow></math> . The mean peak values of <math id="M288"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> with respect to time for fractional orders <math id="M289"><mrow><mi>α</mi><mo>=</mo><mn>0.96</mn></mrow></math> , 0.98, and 0.99 are presented in Table 5 Table 5 with graphical representation in Figs. 9 Figs. 9 . 1 1 ρ 2 1 α = 0.96 2 1 α = 0.98 2 1 α = 0.99 1 1 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 Peak 1 1 Time (s) 1 1 <math id="M293"><mrow><mn>0.005333</mn><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> 1 1 38.0005 1 1 38.28 1 1 45.8029 1 1 38.18 1 1 49.1345 1 1 38.99 6. Conclusion 1 In this study, fractional SPKEs have been solved using the higher order approximation scheme with different fractional orders <math id="M299"><mi>α</mi></math> . The obtained numerical solutions for mean <math id="M300"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> 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 10 , 15 15 , 17 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 <math id="M301"><mrow><mi>N</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> for different reactivities. References S. Saha Ray , Fractional calculus with applications for nuclear reactor dynamics , CRC Press , Taylor & Francis, New York, USA , 2015 . Fractional calculus with applications for nuclear reactor dynamics J.G. Hayes , E.J. Allen , Stochastic point-kinetics equations in nuclear reactor dynamics , Ann. Nucl. Energy . 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 G. Espinosa-Paredes , M.-A. Polo-Labarrios , Erick-G. Espinosa-Martinez et al ., Fractional neutron point kinetics equations for nuclear reactor dynamics , Ann. Nucl. Energy . 38 , 307 - 330 ( 2011 ). doi: 10.1016/j.anucene.2010.10.012 http://doi.org/10.1016/j.anucene.2010.10.012 Fractional neutron point kinetics equations for nuclear reactor dynamics R.-I. Cázares-Ramírez , V. A. Vyawahare , G. Espinosa-Paredes et al ., On the feedback stability of linear FNPK equations , Prog. Nucl. Energy . 98 , 45 - 58 ( 2017 ). doi: 10.1016/j.pnucene.2017.02.007 http://doi.org/10.1016/j.pnucene.2017.02.007 On the feedback stability of linear FNPK equations G. Espinosa-Paredes , M. A. Polo-Labarrios , J. Alvarez-Ramirez , Anomalous diffusion processes in nuclear reactors , Ann. Nucl. Energy . 54 , 227 - 232 ( 2013 ). doi: 10.1016/j.anucene.2012.11.024 http://doi.org/10.1016/j.anucene.2012.11.024 Anomalous diffusion processes in nuclear reactors G. Espinosa-Paredes , R. Vázquez-Rodríguez , E. V. Gallegos et al ., Fractional-space law for the neutron current density , Ann. Nucl. Energy . 55 , 120 - 125 ( 2013 ). doi: 10.1016/j.anucene.2012.12.009 http://doi.org/10.1016/j.anucene.2012.12.009 Fractional-space law for the neutron current density T. K. Nowak , K. Duzinkiewicz , R. Piotrowski , Fractional neutron point kinetics equations for nuclear reactor dynamics – Numerical solution investigations, , Ann. Nucl. Energy . 73 , 317 - 329 ( 2014 ). doi: 10.1016/j.anucene.2014.07.001 http://doi.org/10.1016/j.anucene.2014.07.001 Fractional neutron point kinetics equations for nuclear reactor dynamics – Numerical solution investigations, S. Saha Ray , Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor , Ann. Nucl. Energy . 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. Saha Ray , A. Patra , Numerical solution of fractional stochastic point kinetic equation for nuclear reactor dynamics , Ann. Nucl. Energy . 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 point kinetic equation for nuclear reactor dynamics 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 J. J. Duderstadt , L. J. Hamilton , Nuclear Reactor Analysis , John Wiley and Sons , New York , 1976 . Nuclear Reactor Analysis W. M. Stacey , Nuclear reactor physics , 2 nd Edition, Wiley-VCH Verlag GmbH & Co. KGaA , Weinheim , 2007 . Nuclear reactor physics A. Patra , S. Saha Ray , Solution of non-linear neutron point kinetics equation with feedback reactivity in nuclear reactor dynamics, Internat . J. Nucl. Energy Sci. Tech . 9 , 23 - 24 ( 2015 ). doi: 10.1504/IJNEST.2015.067803 http://doi.org/10.1504/IJNEST.2015.067803 Solution of non-linear neutron point kinetics equation with feedback reactivity in nuclear reactor dynamics, Internat A. Patra , S. Saha Ray , On the solution of the nonlinear fractional neutron point-kinetics equation with Newtonian temperature feedback reactivity , Nucl. Tech . 189 , 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 S. Singh , S. Saha 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. Energy . 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 A. E. Aboanber , A. A. Nahla , A. M. Edress , Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics , Nucl. Sci. Tech . 29 , 132 ( 2018 ). doi: 10.1007/s41365-018-0467-0 http://doi.org/10.1007/s41365-018-0467-0 Developed mathematical technique for fractional stochastic point kinetics model in nuclear reactor dynamics C. Li , R. Wu , H. Ding , High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations , Commun Appl. Ind. Math . 6 , e - 536 ( 2015 ). doi: 10.1685/journal.caim.536 http://doi.org/10.1685/journal.caim.536 High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations H. Li , J. Cao , C. Li , High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equations (III) , J. Comput. Appl. Math . 299 , 159 - 175 ( 2016 ). doi: 10.1016/j.cam.2015.11.037 http://doi.org/10.1016/j.cam.2015.11.037 High-order approximation to Caputo derivatives and Caputo-type advection–diffusion equations (III) J. X. Cao , C. P. Li , Y. Q. Chen , High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (II), Fract . Calc. Appl. Anal . 18 , 735 - 761 ( 2015 ). doi: 10.1515/fca-2015-0045 http://doi.org/10.1515/fca-2015-0045 High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (II), Fract A. A. Nahla , A. M. Edress , Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors , Nucl. Sci. Tech . 27 , 20 ( 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 , Stochastic model for the nonlinear point reactor kinetics equations in the presence Newtonian temperature feedback effects , J. Difference 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 10.1007/s41365-019-0566-6 2018-04-06 2018-08-07 2018-08-10 2019-02-18 2019-03-01 2021-04-10 © 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. 2019 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-019-0566-6 http://dx.doi.org/10.1007/s41365-019-0566-6 . 2019 Nuclear Science and Techniques 03 30