Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors

Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors 1 corresp first-author Corresponding author, a.nahla@science.tanta.edu.eg Corresponding author, a.nahla@science.tanta.edu.eg a.nahla@science.tanta.edu.eg 1 Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt The stochastic point kinetics equations with a multi-group of delayed neutrons, which are the system of a couple stiff stochastic differential equations, are presented. The analytical exponential model is used to solve the stochastic point kinetics equations in the dynamical system of the nuclear reactor. This method is based on the eigenvalues and corresponding eigenvectors of the coefficient matrix. The analytical exponential model calculates the mean and standard deviations of neutrons and precursor populations for the stochastic point kinetics equations with step, ramp, and sinusoidal reactivities. The results of the analytical exponential model are compared with published methods and the results of the deterministic point kinetics model. This comparison confirms that the analytical exponential model is an efficient method for solving stochastic stiff point kinetics equations. Stochastic differential equation Nuclear reactor dynamics Multi-group precursor concentration 1 Introduction 1 The point reactor kinetics equations are the most essential model in the field of nuclear engineering. This system is a coupled linear differential equation and describes the neutron population density and the precursor concentration of delayed neutrons at the center point of homogenous reactors. A point reactor is a homogenous reactor in which the spatial effects have been eliminated. This is obviously possible if the homogenous reactor’s length is infinite in all spatial dimensions. The dynamical process explained by the point kinetics equation is stochastic in nature. The neutron population density and delayed neutron precursor concentrations differ randomly with respect to time. At the levels of high power, the random behavior is imperceptible. But at low-power levels, such as at the beginning, random fluctuation in the neutron population density and neutron precursor concentrations can be crucial. The aim of this work presents an accurate method for stochastic point kinetics equations with step, ramp, and sinusoidal reactivities. There are some techniques which are used for stochastic point reactor kinetics equations. The first of these techniques is the stochastic Piecewise Constant Approximation (stochastic PCA) method and Monte Carlo computations, which were used to calculate the neutron population density and sum of the precursors concentration population density for different values of step reactivity [ 1 1 , 2 2 ]. A simplified stochastic model based on the forward stochastic model in the stochastic kinetics theory and the Itô stochastic differential equations was developed for treating monoenergetic space-time nuclear reactor kinetics in one dimension [ 3 3 , 4 4 ]. Simulation and experimental study of a random neutron analyzing system with a 252 Cf neutron source was presented in Ref. [ 5 5 ]. The Euler-Maruyama and Taylor 1.5 strong order methods were presented for solving stochastic point kinetics equations with step and sinusoidal reactivities [ 6 6 , 7 7 ]. Finally, simplified stochastic point kinetics equations (SSPK) were modeled with a system of Itô stochastic differential equations. This approach does not require computing the square root of a matrix, which is a great computational advantage [ 8 8 ]. In this work, the system of stochastic point kinetics model in the dynamic nuclear reactor is derived in Sec. 2. The analytical exponential model is presented and applied to solve the stochastic point kinetics equations in Sec. 3. The numerical results of the proposed method are discussed and compared with the traditional methods in Sec. 4. Finally, the conclusions and future work are discussed in Sec. 5. 2 Stochastic Model 1 The one speed neutron diffusion equations with a multi-group of delayed precursor concentrations are written as Ref. [ 9 9 - 11 11 9-11 ]: <math display="block" id="M1"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mn>1</mn><mi>v</mi></mfrac><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>D</mi></mrow></mtd><mtd columnalign="left"><mrow><msup><mo>∇</mo><mn>2</mn></msup><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>−</mo><mo stretchy="false">(</mo><msub><mi>Σ</mi><mi>a</mi></msub><mo>−</mo><msub><mi>Σ</mi><mi>f</mi></msub><mo stretchy="false">)</mo><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mo stretchy="false">[</mo><mi>ν</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mo>−</mo><mn>1</mn><mo stretchy="false">]</mo><msub><mi>Σ</mi><mi>f</mi></msub><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><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><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>S</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M2"><mrow><mfrac><mo>∂</mo><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><msub><mi>β</mi><mi>i</mi></msub><mi>ν</mi><msub><mi>Σ</mi><mi>f</mi></msub><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><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><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>i</mi><mo>=</mo><mn>1,2,3,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi></mrow></math> where Φ ( r , t ) is the neutron population, r is the position (cm), t is the time (s), C i ( r , t ) is the i -group of delayed precursor concentration, S 0 ( r , t ) is the external neutron source, D is the diffusion coefficient, ∑ a is the absorption cross sections, ∑ f is the fission cross sections, ν is the neutron fission, v is the neutron speed, λ i is the decay constant of i -group of delayed neutrons, β i is the fraction of i -group delayed neutrons, and <math id="M3"><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. The neutron population and the i -group of the delayed precursor concentration, using separation of variables, can be written as: <math display="block" id="M4"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>=</mo><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>C</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>=</mo><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>S</mi><mn>0</mn></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>=</mo><mi>q</mi><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> where n ( t ) is the neutron population density as a function of time only and c i ( t ) is the precursor concentration density of delayed neutrons. The function <math id="M5"><mrow><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo></mrow></math> is the fundamental function, which can be determined from the following <math display="block" id="M6"><mrow><msup><mo>∇</mo><mn>2</mn></msup><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mo>+</mo><msubsup><mi>B</mi><mi>g</mi><mn>2</mn></msubsup><mover accent="true"><mi>Ψ</mi><mo stretchy="true">^</mo></mover><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo stretchy="false">)</mo><mo>=</mo><mn>0,</mn></mrow></math> where <math id="M7"><mrow><msubsup><mi>B</mi><mi>g</mi><mn>2</mn></msubsup></mrow></math> is the geometric buckling. Substituting Eqs. (3) and (4) into Eqs. (1) and (2) we get <math display="block" id="M8"><mtable columnalign="left"><mtr><mtd><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><mo>=</mo><mo>−</mo><mo stretchy="false">[</mo><mi>D</mi><msubsup><mi>B</mi><mi>g</mi><mn>2</mn></msubsup><mo>+</mo><msub><mi>Σ</mi><mi>a</mi></msub><mo>−</mo><msub><mi>Σ</mi><mi>f</mi></msub><mo stretchy="false">]</mo><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><mo>+</mo><mo stretchy="false">[</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mi>ν</mi><mo>−</mo><mn>1</mn><mo stretchy="false">]</mo><msub><mi>Σ</mi><mi>f</mi></msub><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><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><msub><mi>c</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><mi>q</mi><mn>,</mn></mtd></mtr></mtable></math> <math display="block" id="M9"><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><mo>=</mo><msub><mi>β</mi><mi>i</mi></msub><mi>ν</mi><msub><mi>Σ</mi><mi>f</mi></msub><mi>v</mi><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><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>i</mi><mo>=</mo><mn>1,2,3,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mn>.</mn></mrow></math> According to Ref. [ 1 1 , 2 2 , 11 11 ], Eqs. (5) and (6) can be separated into four terms: deaths, births, decay, and external sources. Therefore, <math display="block" id="M10"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><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><mo>=</mo><mo>−</mo></mrow></mtd><mtd columnalign="left"><mrow><munder><munder><mrow><mo stretchy="false">[</mo><mi>D</mi><msubsup><mi>B</mi><mi>g</mi><mn>2</mn></msubsup><mo>+</mo><msub><mi>Σ</mi><mi>a</mi></msub><mo>−</mo><msub><mi>Σ</mi><mi>f</mi></msub><mo stretchy="false">]</mo><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Deaths</mtext></mrow></munder></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><munder><munder><mrow><mo stretchy="false">[</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mi>ν</mi><mo>−</mo><mn>1</mn><mo stretchy="false">]</mo><msub><mi>Σ</mi><mi>f</mi></msub><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Births</mtext></mrow></munder></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><munder><munder><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><mo stretchy="true">︸</mo></munder><mrow><mtext>Decay</mtext></mrow></munder><mo>+</mo><munder><munder><mi>q</mi><mo stretchy="true">︸</mo></munder><mrow><mtext>Source</mtext></mrow></munder><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M11"><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><mo>=</mo><munder><munder><mrow><msub><mi>β</mi><mi>i</mi></msub><mi>ν</mi><msub><mi>Σ</mi><mi>f</mi></msub><mi>v</mi><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Births</mtext></mrow></munder><mo>−</mo><munder><munder><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><mo stretchy="true">︸</mo></munder><mrow><mtext>Decay</mtext></mrow></munder><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>i</mi><mo>=</mo><mn>1,2,3,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mn>.</mn></mrow></math> Equations (7) and (8) can be rewritten as deterministic point kinetics equations using the reactivity <math id="M12"><mrow><mi>ρ</mi><mo>=</mo><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mo>∑</mo><mi>a</mi></msub><mo>+</mo><mi>D</mi><msubsup><mi>B</mi><mi>g</mi><mn>2</mn></msubsup></mrow><mrow><mi>ν</mi><msub><mo>∑</mo><mi>f</mi></msub></mrow></mfrac></mrow></math> , which is a function of time, the generation time of neutrons <math id="M13"><mrow><mi>Λ</mi><mo>=</mo><mfrac><mn>1</mn><mrow><mi>ν</mi><msub><mo>∑</mo><mi>f</mi></msub><mi>v</mi></mrow></mfrac></mrow></math> , and the constant <math id="M14"><mrow><mi>α</mi><mo>=</mo><mfrac><mn>1</mn><mi>ν</mi></mfrac></mrow></math> . <math display="block" id="M15"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><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><mo>=</mo><mo>−</mo></mrow></mtd><mtd columnalign="left"><mrow><munder><munder><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mn>1</mn><mo>−</mo><mi>ρ</mi><mo>−</mo><mi>α</mi></mrow><mi>Λ</mi></mfrac></mrow><mo>)</mo></mrow><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Deaths</mtext></mrow></munder><mo>+</mo><munder><munder><mrow><mrow><mo>(</mo><mrow><mfrac><mrow><mn>1</mn><mo>−</mo><mi>α</mi><mo>−</mo><mi>β</mi></mrow><mi>Λ</mi></mfrac></mrow><mo>)</mo></mrow><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Births</mtext></mrow></munder></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover></mstyle><munder><munder><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><mo stretchy="true">︸</mo></munder><mrow><mtext>Decay</mtext></mrow></munder><mo>+</mo><munder><munder><mi>q</mi><mo stretchy="true">︸</mo></munder><mrow><mtext>Source</mtext></mrow></munder><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M16"><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><mo>=</mo><munder><munder><mrow><mfrac><mrow><msub><mi>β</mi><mi>i</mi></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mo stretchy="true">︸</mo></munder><mrow><mtext>Births</mtext></mrow></munder><mo>−</mo><munder><munder><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><mo stretchy="true">︸</mo></munder><mrow><mtext>Decay</mtext></mrow></munder><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>i</mi><mo>=</mo><mn>1,2,3,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mn>.</mn></mrow></math> Note that, n ( t ) is the population size of neutrons. c i ( t ) is the population size of the i th group of delayed precursor concentration. Equations (9) and (10) are the deterministic point kinetics equations where terms are separated into births, deaths, decay, and external source q . Indeed, the neutron birth rate due to fission is <math id="M17"><mrow><mi>b</mi><mo>=</mo><mfrac><mrow><mn>1</mn><mo>−</mo><mi>α</mi><mo>−</mo><mi>β</mi></mrow><mrow><mi>Λ</mi><mo stretchy="false">[</mo><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mi>ν</mi><mo>−</mo><mn>1</mn><mo stretchy="false">]</mo></mrow></mfrac><mo>=</mo><mfrac><mn>1</mn><mrow><mi>ν</mi><mi>Λ</mi></mrow></mfrac></mrow></math> , where [(1- β ) ν -1] is the number of new born neutrons in each fission. The neutron death rate due to captures and leakage is <math id="M18"><mrow><mi>d</mi><mo>=</mo><mfrac><mrow><mn>1</mn><mo>−</mo><mi>ρ</mi><mo>−</mo><mi>α</mi></mrow><mi>Λ</mi></mfrac></mrow></math> . Also, λ i c i is the rate of the i th group of delayed precursor decay. To derive the stochastic dynamical system, let us take a small time interval size h where the probability of more than one event occurring during this time interval is small. During time interval h , there are ( I +3) different possibilities for an event. The change in populations n and c i during this time interval is <math display="block" id="M19"><mrow><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow><mo>⇔</mo><mo>〈</mo><mi>Δ</mi><mi>P</mi><mo>|</mo><mo>=</mo><mrow><mo>(</mo><mrow><mi>Δ</mi><mi>n</mi><mtext> </mtext><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub><mtext> </mtext><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub><mtext> </mtext><mo>⋯</mo><mtext> </mtext><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> In the present investigation let us assume that the changes are approximately normally distributed. The ( I +3) possibilities and its probabilities are Ref. [ 1 1 , 2 2 ]: <math display="block" id="M20"><mrow><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mn>1</mn></msub><mo>〉</mo><mo>=</mo><msub><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mn>1</mn></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mo>−</mo><mn>1</mn></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><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mn>1</mn></msub><mo>=</mo><mi>h</mi><mi>d</mi><mi>n</mi><mn>,</mn></mrow></math> <math display="block" id="M21"><mrow><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mn>2</mn></msub><mo>〉</mo><mo>=</mo><msub><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mn>2</mn></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo><mi>ν</mi><mo>−</mo><mn>1</mn></mtd></mtr><mtr><mtd><msub><mi>β</mi><mn>1</mn></msub><mi>ν</mi></mtd></mtr><mtr><mtd><msub><mi>β</mi><mn>2</mn></msub><mi>ν</mi></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mi>β</mi><mn>1</mn></msub><mi>ν</mi></mtd></mtr></mtable><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mn>2</mn></msub><mo>=</mo><mi>h</mi><mi>b</mi><mi>n</mi></mrow></math> <math display="block" id="M22"><mrow><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mn>3</mn></msub><mo>〉</mo><mo>=</mo><msub><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mn>3</mn></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mo>−</mo><mn>1</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mn>3</mn></msub><mo>=</mo><mi>h</mi><msub><mi>λ</mi><mn>1</mn></msub><msub><mi>c</mi><mn>1</mn></msub><mn>,</mn></mrow></math> <math display="block" id="M23"><mrow><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mn>4</mn></msub><mo>〉</mo><mo>=</mo><msub><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mn>4</mn></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>−</mo><mn>1</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mn>4</mn></msub><mo>=</mo><mi>h</mi><msub><mi>λ</mi><mn>2</mn></msub><msub><mi>c</mi><mn>2</mn></msub><mn>,</mn></mrow></math> <math display="block" id="M24"><mtable><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mrow><mi>I</mi><mo>+</mo><mn>2</mn></mrow></msub><mo>〉</mo><mo>=</mo><msub><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow><mrow><mi>I</mi><mo>+</mo><mn>2</mn></mrow></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mo>−</mo><mn>1</mn></mtd></mtr></mtable><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mrow><mi>I</mi><mo>+</mo><mn>2</mn></mrow></msub><mo>=</mo><mi>h</mi><msub><mi>λ</mi><mi>I</mi></msub><msub><mi>c</mi><mi>I</mi></msub><mn>,</mn></mtd></mtr></mtable></math> <math display="block" id="M25"><mrow><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mrow><mi>I</mi><mo>+</mo><mn>3</mn></mrow></msub><mo>〉</mo><mo>=</mo><msub><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mi>Δ</mi><mi>n</mi></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mi>Δ</mi><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mrow><mi>I</mi><mo>+</mo><mn>3</mn></mrow></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mn>1</mn></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><mo>)</mo></mrow><mn>,</mn><mtext> </mtext><mtext> </mtext><msub><mi>p</mi><mrow><mi>I</mi><mo>+</mo><mn>3</mn></mrow></msub><mo>=</mo><mi>h</mi><mi>q</mi><mn>.</mn></mrow></math> In the present analysis, it is assumed that the extraneous source produces neutrons randomly following a Poisson process with intensity q . According to these assumptions, the changes in neutron population and precursor concentration are approximately normally distributed with a mean of <math display="block" id="M26"><mtable columnalign="left"><mtr><mtd><mi>E</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>I</mi><mo>+</mo><mn>3</mn></mrow></munderover></mstyle><msub><mi>p</mi><mi>k</mi></msub><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mi>k</mi></msub><mo>〉</mo><mo>=</mo><mi>h</mi><mrow><mo>(</mo><mtable><mtr><mtd><mfrac><mrow><mi>ρ</mi><mo>−</mo><mi>β</mi></mrow><mi>Λ</mi></mfrac><mi>n</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><msub><mi>c</mi><mi>i</mi></msub><mo>+</mo><mi>q</mi></mtd></mtr><mtr><mtd><mfrac><mrow><msub><mi>β</mi><mn>1</mn></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>−</mo><msub><mi>λ</mi><mn>1</mn></msub><msub><mi>c</mi><mn>1</mn></msub></mtd></mtr><mtr><mtd><mfrac><mrow><msub><mi>β</mi><mn>2</mn></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>−</mo><msub><mi>λ</mi><mn>2</mn></msub><msub><mi>c</mi><mn>2</mn></msub></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><mfrac><mrow><msub><mi>β</mi><mi>I</mi></msub></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>−</mo><msub><mi>λ</mi><mi>I</mi></msub><msub><mi>c</mi><mi>I</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow><mn>.</mn></mtd></mtr><mtr><mtd><mrow/></mtd></mtr></mtable></math> According to Ref. [ 8 8 ], the variance takes the form <math display="block" id="M27"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo><mo>≃</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>I</mi><mo>+</mo><mn>3</mn></mrow></munderover></mstyle><msub><mi>p</mi><mi>k</mi></msub><mo>|</mo><mi>Δ</mi><msub><mi>P</mi><mi>k</mi></msub><mo>〉</mo><mo>〈</mo><mi>Δ</mi><msub><mi>P</mi><mi>k</mi></msub><mo>|</mo><mo>=</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mn>,</mn></mrow></math> where B is the diagonal matrix as <math display="block" id="M28"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mfrac><mrow><mo>−</mo><mn>1</mn><mo>−</mo><mi>ρ</mi><mo>+</mo><mn>2</mn><mi>β</mi><mo>+</mo><msup><mrow><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>β</mi><mo stretchy="false">)</mo></mrow><mn>2</mn></msup><mi>ν</mi></mrow><mi>Λ</mi></mfrac><mi>n</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><msub><mi>c</mi><mi>i</mi></msub><mo>+</mo><mi>q</mi></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mrow><mfrac><mrow><msubsup><mi>β</mi><mn>1</mn><mn>2</mn></msubsup><mi>ν</mi></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>+</mo><msub><mi>λ</mi><mn>1</mn></msub><msub><mi>c</mi><mn>1</mn></msub></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mrow><mfrac><mrow><msubsup><mi>β</mi><mn>2</mn><mn>2</mn></msubsup><mi>ν</mi></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>+</mo><msub><mi>λ</mi><mn>2</mn></msub><msub><mi>c</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><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mo>⋯</mo></mtd><mtd><mrow><mfrac><mrow><msubsup><mi>β</mi><mi>I</mi><mn>2</mn></msubsup><mi>ν</mi></mrow><mi>Λ</mi></mfrac><mi>n</mi><mo>+</mo><msub><mi>λ</mi><mi>I</mi></msub><msub><mi>c</mi><mi>I</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> Using central limit theorem, the random variate <math id="M29"><mrow><mfrac><mrow><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo>−</mo><mi>E</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo></mrow><mrow><msqrt><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo></mrow></msqrt></mrow></mfrac><mo>=</mo><mo>|</mo><mi>η</mi><mo>〉</mo></mrow></math> follows standard normal distribution [ 1 1 , 2 2 , 12 12 ]. This implies <math display="block" id="M30"><mrow><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo>=</mo><mi>E</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo><mo>+</mo><msqrt><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mo>|</mo><mi>Δ</mi><mi>P</mi><mo>〉</mo><mo stretchy="false">)</mo></mrow></msqrt><mo>|</mo><mi>η</mi><mo>〉</mo><mn>,</mn></mrow></math> where <math id="M31"><mrow><mo>|</mo><mi>η</mi><mo>〉</mo><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mi>η</mi><mn>1</mn></mtd></mtr><mtr><mtd><mi>η</mi><mn>2</mn></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><msub><mi>η</mi><mrow><mi>I</mi><mo>+</mo><mn>1</mn></mrow></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow></math> , and η 1 , η 2 , η I +1 N (0,1). Substituting Eqs. (18) and (19) into Eq. (20), we have http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339366&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339364&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339366&type= This Eq. (21) gives, as h →0, the following Itô stochastic differential equation system <math display="block" id="M32"><mrow><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>|</mo><mi>P</mi><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><mi>P</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mo>|</mo><mi>Q</mi><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><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>|</mo><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mn>,</mn></mrow></math> where, http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339377&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339373&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=10339377&type= and <math id="M33"><mrow><mo>|</mo><mi>Δ</mi><mi>W</mi><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> where, W 1 ( t ), W 2 ( t ),, W I +1 ( t ) are Wiener processes [ 6 6 , 7 7 ]. Equation (22) represents the stochastic point reactor kinetics model. Since B=0, it reduces to the deterministic point kinetics model. This model was solved analytically and numerically in many references, for example, Ref. [ 13 13 - 18 18 13-18 ]. 3 Analytical Exponential Model 1 Using the integration factor, differential Eq. (22) in matrix form is rewritten as <math display="block" id="M34"><mtable columnalign="left"><mtr><mtd><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mi>t</mi><mo stretchy="false">)</mo><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mi>t</mi><mo stretchy="false">)</mo><mo>|</mo><mi>Q</mi><mo>〉</mo></mtd></mtr><mtr><mtd><mo>+</mo><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mi>t</mi><mo stretchy="false">)</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mfrac><mtext>d</mtext><mrow><mtext>d</mtext><mi>t</mi></mrow></mfrac><mo>|</mo><mi>W</mi><mo stretchy="false">(</mo><mi>t</mi><mo stretchy="false">)</mo><mo>〉</mo><mn>.</mn></mtd></mtr></mtable></math> Let us divide the time into M , very small time intervals with step size h . The matrices A and B are constants during the time interval [ t m , t m +1 ] where t m +1 = t m + h , and m =0,1,2,⋅s, M -1. Equation (24) can then take the following form <math display="block" id="M35"><mtable columnalign="left"><mtr><mtd><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msub><mi>t</mi><mrow><mi>m</mi><mo>+</mo><mn>1</mn></mrow></msub><mo stretchy="false">)</mo><mo>|</mo><mi>P</mi><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><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>〉</mo></mtd></mtr><mtr><mtd><mo>=</mo><mi>h</mi><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>|</mo><mi>Q</mi><mo>〉</mo><mo>+</mo><msqrt><mi>h</mi></msqrt><mi>exp</mi><mo stretchy="false">(</mo><mo>−</mo><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo>|</mo><mi>η</mi><mo>〉</mo><mn>,</mn></mtd></mtr></mtable></math> and consequently, we get <math display="block" id="M36"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mo>|</mo><mi>P</mi><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><mi>exp</mi><mo stretchy="false">(</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">)</mo><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>〉</mo></mrow></mtd><mtd columnalign="left"><mrow><mo>+</mo><mi>h</mi><mi>exp</mi><mo stretchy="false">(</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">)</mo><mo>|</mo><mi>Q</mi><mo>〉</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msqrt><mi>h</mi></msqrt><mi>exp</mi><mo stretchy="false">(</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">)</mo><msup><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo>|</mo><mi>η</mi><mo>〉</mo><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> The mathematical treatment of this system can be found by calculating all the eigenvalues and corresponding eigenvectors of matrix A and performing straightforward computations. However, this is an expensive scheme, especially when the reactivity varies with time, since the eigenvalues of matrix A are calculated by solving the Inhour equation, a ( I +1)th-order algebraic equation, at each time step. The eigenvectors of A , denoted by ket vectors | U j ⟩ and the corresponding eigenvalues denoted by ω j , obey the relation <math display="block" id="M37"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><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><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mtext> </mtext><mi>j</mi><mo>=</mo><mn>0,1,2,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mn>.</mn></mrow></math> The eigenvalues ω j of matrix A are the roots of the inhour formula <math display="block" id="M38"><mrow><mi>ρ</mi><mo>=</mo><mi>Λ</mi><mi>ω</mi><mo>+</mo><mi>ω</mi><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></mrow><mrow><msub><mi>λ</mi><mi>i</mi></msub><mo>+</mo><mi>ω</mi></mrow></mfrac><mn>.</mn></mrow></math> For j , an arbitrary function f ( h A ) satisfies the following expression [ 13 13 - 16 16 13-16 ] <math display="block" id="M39"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><mo stretchy="false">)</mo><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>=</mo><mi>f</mi><mo stretchy="false">(</mo><mi>h</mi><msub><mi>ω</mi><mi>j</mi></msub><mo stretchy="false">)</mo><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mn>,</mn><mtext> </mtext><mtext> </mtext><mtext> </mtext><mo>∀</mo><mi>j</mi><mo>=</mo><mn>0,1,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi><mo>;</mo></mrow></math> and, consequently, we get <math display="block" id="M40"><mrow><mi>exp</mi><mo stretchy="false">(</mo><mi>h</mi><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><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><msup><mi>e</mi><mrow><mi>h</mi><msub><mi>ω</mi><mi>j</mi></msub></mrow></msup><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</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></mrow></math> with the properties <math display="block" id="M41"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>A</mi></mstyle><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><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mn>,</mn><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></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mtext>and</mtext><mtext> </mtext><mo>〈</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>k</mi></msub><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</mi></mstyle><mi>j</mi></msub><mo>〉</mo><mo>=</mo><msub><mi>δ</mi><mrow><mi>k</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>k</mi><mtext>=j;</mtext></mtd></mtr><mtr><mtd><mtext>0, </mtext><mi>k</mi><mo>≠</mo><mtext>j</mtext><mtext>.</mtext></mtd></mtr></mtable></mrow></mrow></mtd></mtr></mtable></mrow></math> The ket eigenvectors | U j ⟩ and the bra eigenvectors U j | are calculated analytically [ 13 13 - 16 16 13-16 ] <math display="block" id="M42"><mtable columnalign="left"><mtr><mtd><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><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><mtable><mtr><mtd><mn>1</mn></mtd></mtr><mtr><mtd><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></mtd></mtr><mtr><mtd><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></mtd></mtr><mtr><mtd><mo>⋮</mo></mtd></mtr><mtr><mtd><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></mtd></mtr></mtable><mo>)</mo></mrow><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><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><mtext> </mtext><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><mtext> </mtext><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><mn>.</mn></mrow></mtd></mtr></mtable></mtd></mtr><mtr><mtd><mrow/></mtd></mtr></mtable></math> For the normalized condition ⟨ U j | U j ⟩=1, we get <math id="M43"><mrow><msub><mi>Ω</mi><mi>j</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><msqrt><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></msqrt></mrow></mfrac><mn>,</mn><mtext> </mtext><mtext> </mtext><mo>=</mo><mn>0,1,2,</mn><mo>⋯</mo><mn>,</mn><mi>I</mi></mrow></math> . Hence, the general solution of the stochastic point kinetics equation takes the following form <math display="block" id="M44"><mrow><mo>|</mo><mi>P</mi><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><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>0</mn></mrow><mi>I</mi></munderover></mstyle><msup><mi>e</mi><mrow><mi>h</mi><msub><mi>ω</mi><mi>j</mi></msub></mrow></msup><mo>|</mo><msub><mstyle mathvariant="bold" mathsize="normal"><mi>U</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 stretchy="false">[</mo><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mi>m</mi></msub><mo stretchy="false">)</mo><mo>〉</mo><mo>+</mo><mi>h</mi><mo>|</mo><mi>Q</mi><mo>〉</mo><mo>+</mo><msqrt><mi>h</mi></msqrt><msup><mstyle mathvariant="bold" mathsize="normal"><mi>B</mi></mstyle><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac></mrow></msup><mo>|</mo><mi>η</mi><mo>〉</mo><mo stretchy="false">]</mo><mn>,</mn></mrow></math> where the initial condition is <math id="M45"><mrow><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><mn>0</mn><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><mo>|</mo><mi>P</mi><mo stretchy="false">(</mo><msub><mi>t</mi><mn>0</mn></msub><mo stretchy="false">)</mo><mo>〉</mo><mo>=</mo><msup><mrow><mrow><mo>(</mo><mrow><msub><mi>n</mi><mn>0</mn></msub><mtext> </mtext><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><mtext> </mtext><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><mtext> </mtext><mo>⋯</mo><mtext> </mtext><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><mo>)</mo></mrow></mrow><mi>T</mi></msup></mrow></math> . 4 Numerical Results and Discussion 1 In this section, many of the examples are presented to measure the accuracy of the analytical exponential model (AEM) for stochastic point kinetics equations. The mean and standard deviation of the neutron and precursor population are calculated by solving the stochastic point kinetics equations with three different cases: step, ramp, and sinusoidal reactivities, and are compared with the published stochastic methods. The first benchmark problem does not model an actual physical nuclear reactor problem, but this problem provides a simple computational solution for comparing the stochastic model. This model simulates a step reactivity insertion and assumes one neutron precursor [ 1 1 ]. The parameters of this benchmark are: the neutron generation time <math id="M46"><mrow><mi>Λ</mi><mo>=</mo><mfrac><mn>2</mn><mn>3</mn></mfrac></mrow></math> s, reactivity <math id="M47"><mrow><mi>ρ</mi><mo>=</mo><mo>−</mo><mfrac><mn>1</mn><mn>3</mn></mfrac></mrow></math> , decay constant λ 1 =0.1 s -1 , fraction delayed neutron β 1 = β =0.05, number of neutrons per fission ν =2.5, external source q =200 s -1 , and the initial condition assumes n (0)=400 and c (0)=300. Table 1 Table 1 presents the mean and standard deviations of neutron and precursor populations using 40 time intervals for a time interval of length of 2. The results of the AEM are compared with the results of the Monte Carlo, stochastic PCA method [ 1 1 ], Euler-Maruyama [ 6 6 ], Taylor 1.5 strong order [ 6 6 ] and the deterministic point kinetics model (DPKM) i.e B=0 at time t =2 s. All the stochastic methods use 5000 trails. Table 1 Table 1 shows the accuracy of the analytical exponential model such that the results of the AEM are in good agreement with the deterministic point kinetics model, more than other methods under the same conditions. This example confirms that the AEM is an efficient method for solving stochastic point kinetics equations in nuclear reactor dynamics. 1 1   1 1 Monte Carlo 1 1 Stochastic PCA 1 1 Euler-Maruyama 1 1 Taylor 1.5 Strong order 1 1 AEM 1 1 DPKM(B=0) 1 1 E [ n (2)] 1 1 400.03 1 1 395.32 1 1 412.23 1 1 412.10 1 1 396.28 1 1 n (2)=396.63 1 1 [ n (2)] 1 1 27.311 1 1 29.411 1 1 34.391 1 1 34.519 1 1 31.212 1 1   1 1 E [ c (2)] 1 1 300.00 1 1 300.67 1 1 315.96 1 1 315.93 1 1 300.42 1 1 c (2)=300.40 1 1 [ c (2)] 1 1 7.8073 1 1 8.3564 1 1 8.2656 1 1 8.3158 1 1 7.9576 1 1   The second benchmark problem simulates a step reactivity insertion for an actual nuclear reactor with six groups of delayed neutrons [ 1 1 ]. The parameters of this reactor are as follows: β i =[0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182], β =0.007, λ i =[0.0127, 0.0317, 0.115, 0.311, 1.4, 3.87] s -1 , Λ =0.00002 s, ν =2.5, n 0 =100, and no external source q =0 s -1 . Two cases of step reactivity insertion, ρ =0.003 and ρ =0.007, are presented. Calculation results of the mean and standard deviations can be seen in Table 2 Table 2 at time 0.1 s for the first case ρ =0.003 and at time 0.001 s for the second case ρ =0.007. The results of the AEM are compared with the results of Monte Carlo, Stochastic PCA method [ 1 1 ], Euler-Maruyama [ 6 6 ], Taylor 1.5 Strong order [ 6 6 ], simplified stochastic point kinetics (SSPK) method [ 8 8 ], and DPKM (B=0). All the stochastic methods use 5000 trails. The agreement is seen between the results of the AEM, deterministic point kinetics model, Monte Carlo, Stochastic PCA, and SSPK methods in Table 2 Table 2 for two cases of reactivity. 1 1   1 1 Monte Carlo 1 1 Stochastic PCA 1 1 Euler-Maruyama 1 1 Taylor 1.5 Strong order 1 1 SSPK 1 1 AEM 1 1 DPKM(B=0) 1 1 ρ =0.003 1 1   1 1   1 1   1 1   1 1   1 1   1 1   1 1 E [ n (0.1)] 1 1 183.04 1 1 186.31 1 1 208.6 1 1 199.408 1 1 184.8 1 1 186.30 1 1 179.95 1 1 σ [ n (0.1)] 1 1 168.79 1 1 164.16 1 1 255.95 1 1 168.547 1 1 186.96 1 1 164.14 1 1   1 1 <math id="M48"><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>0.1</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 4.478E+5 1 1 4.491E+5 1 1 4.498E+5 1 1 4.497E+5 1 1 4.489E+5 1 1 4.490E+5 1 1 4.489E+5 1 1 <math id="M49"><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>0.1</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 1495.7 1 1 1917.2 1 1 1233.38 1 1 1218.82 1 1 982.64 1 1 1911.91 1 1   1 1 ρ =0.007 1 1   1 1   1 1   1 1   1 1   1 1   1 1   1 1 E [ n (0.001)] 1 1 135.67 1 1 134.55 1 1 139.568 1 1 139.569 1 1   1 1 134.54 1 1 135.0 1 1 σ [ n (0.001)] 1 1 93.376 1 1 91.242 1 1 92.042 1 1 92.047 1 1   1 1 91.234 1 1   1 1 <math id="M50"><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>0.001</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math> 1 1 4.464E+5 1 1 4.464E+5 1 1 4.463E+5 1 1 4.463E+5 1 1   1 1 4.464E+5 1 1 4.464E+5 1 1 <math id="M51"><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>0.001</mn><mo stretchy="false">)</mo></mrow></mstyle></mrow><mo>]</mo></mrow></mrow></math>   1 1 7.8073 1 1 8.3564 1 1 8.2656 1 1 8.3158 1 1   1 1 19.235 1 1   The mean neutron population and two individual neutron sample paths are given in Fig. 1 Fig. 1 for a step reactivity insertion ( ρ =0.003) and in Fig. 2 Fig. 2 for a step reactivity insertion ( ρ =0.007). The third benchmark problem simulates a ramp reactivity insertion for the above actual nuclear reactor with six groups of delayed neutrons. The parameters of this reactor are the same as in the second benchmark problem, except the reactivity is ρ =0.1 βt . The mean neutron population and two individual neutron sample paths are shown in Fig. 3 Fig. 3 using a time interval size of h =0.01 s after 5000 trails. The final benchmark problem simulates a sinusoidal reactivity insertion, <math id="M52"><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> [ 7 7 ]. The parameters of the reactor with one group of delayed neutrons are as follows: ρ 0 =0.005333 (0.68 $), λ 1 =0.077 s -1 , β 1 = β =0.0079, Λ =10 -3  s, q =0 s -1 , n 0 =1, and a half period time of T =50 s. The mean neutron population and two individual neutron sample paths are seen in Fig. 4 Fig. 4 using the time interval size h =0.1 s after 5000 trails. 5 Conclusion 1 In this work, stochastic point kinetics equations were introduced, which represent the generalization of the deterministic point kinetics equations. The analytical exponential model was described, which is based on the eigenvalues and corresponding eigenvectors. This method is an efficient approximate for the stiff system of stochastic point kinetics differential equations. The mean and standard deviations of the neutron and precursor populations using the analytical exponential model are in good agreement with these by the deterministic point kinetics equations, Monte Carlo, and stochastic PCA methods more than other references methods. This agreement confirms that the analytical exponential model is an efficient method for solving stochastic point kinetics equations in nuclear reactor dynamics using a step, ramp, and sinusoidal reactivities. Possible future work may include the derivation and study of stochastic multi-energy group reactor kinetics equations where spatial effects are included. References J.G. Hayes , E.J. Allen , Stochastic point kinetics equations in nuclear reactor dynamics . Ann Nucl Energ , 2005 , 32 : 572 - 587 . 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 , Stochastic point kinetics equations in nuclear reactor dynamics. Thesis , Texas Tech University , 2005 . https://repositories.tdl.org/ttu-ir/bitstream/handle/2346/22276/HayesThesis.pdf https://repositories.tdl.org/ttu-ir/bitstream/handle/2346/22276/HayesThesis.pdf 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 , 2010 , 47 : 705 - 711 . 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 , 2011 , 43 : 523 - 530 . 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 P. Feng , S. Liu , B. Wei , et al . Simulation and experimental study of a random neutron analyzing system with 252Cf neutron source . Nucl Sci Tech , 2011 , 22 : 39 - 46 . http://www.j.sinap.ac.cn/nst/CN/Y2011/V22/I1/39www.j.sinap.ac.cn/nst/CN/Y2011/V22/I1/39 http://www.j.sinap.ac.cn/nst/CN/Y2011/V22/I1/39www.j.sinap.ac.cn/nst/CN/Y2011/V22/I1/39 Simulation and experimental study of a random neutron analyzing system with 252Cf neutron source S.S. Ray , Numerical simulation of stochastic point kinetic equation in the dynamical system of nuclear reactor . Ann Nucl Energ , 2012 , 49 : 154 - 159 . 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 stochastic point kinetics equations with sinusoidal reactivity in dynamical system of nuclear reactor . Int J Nucl Energ Sci Technol , 2013 , 7 : 231 - 242 . DOI: 10.1504/IJNEST.2013.052165 http://doi.org/10.1504/IJNEST.2013.052165 Numerical solution for stochastic point kinetics equations with sinusoidal reactivity in dynamical system of nuclear reactor S.M. Ayyoubzadeh , N. Vosoughi , An alternative stochastic formulation for the point kinetics . Ann Nucl Energ , 2014 , 63 : 691 - 695 . 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 W.M. Stacey , Nuclear reactor physics . 2nd Edition , Weinheim (Germany) : WILEY-VCH Verlag GmbH & Co. KGaA , 2007 . Nuclear reactor physics J.J. Duderstadt , L.J. Hamilton , Nuclear Reactor Analysis . New York (USA) : John Wiley and Sons , 1976 . Nuclear Reactor Analysis M.M.R. Williams , Random processes in nuclear reactors . Oxford (UK) : Pergamon Press , 1974 . Random processes in nuclear reactors E.J. Allen , Modeling with Itô stochastic differential equations . Dordrecht (Netherlands) : Springer , 2007 . Modeling with Itô stochastic differential 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 , 2002 , 35 : 3245 - 3263 . 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 , 2002 , 35 : 9609 - 9627 . 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 , 2008 , 238 : 2648 - 2653 . 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 , 2011 , 38 : 2810 - 2817 . 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 H.F. Li , X.L. Shang , W.Z. Chen , An accurate solution of point kinetics equations of one-group delayed neutrons and an extraneous neutron source for step reactivity insertion . Chinese Sci Bull , 2010 , 55 : 4116 - 4119 . DOI: 10.1007/s11434-010-4220-2 http://doi.org/10.1007/s11434-010-4220-2 An accurate solution of point kinetics equations of one-group delayed neutrons and an extraneous neutron source for step reactivity insertion W.Z. Chen , H. Xiao , H.F. Li , et al . Explicit appropriate basis function method for numerical solution of stiff systems . Ann Nucl Energ , 2015 , 75 : 353 - 357 . DOI: 10.1016/j.anucene.2014.08.040 http://doi.org/10.1016/j.anucene.2014.08.040 Explicit appropriate basis function method for numerical solution of stiff systems 10.1007/s41365-016-0025-6 2015-05-31 2015-07-22 2015-07-23 2016-03-02 2016-02-20 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. 2016 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-016-0025-6 http://dx.doi.org/10.1007/s41365-016-0025-6 . 2016 Nuclear Science and Techniques 01 27