# Numerical Approximation for Fractional Neutron Transport Equation

Authors: Zhengang Zhao; Yunying Zheng
Journal: Journal of Mathematics
 (2021)

---

## Abstract

Fractional neutron transport equation reflects the anomalous transport processes in nuclear reactor. In this paper, we will construct the fully discrete methods for this type of fractional equation with Riesz derivative, where the generalized WENO5 scheme is used in spatial direction and Runge–Kutta schemes are adopted in temporal direction. The linear stabilities of the generalized WENO5 schemes with different stages and different order ERK are discussed detailed. Numerical examples show the combinations of forward Euler/two-stage, second-order ERK and WENO5 are unstable and the three-stage, third-order ERK method with generalized WENO5 is stable and can maintain sharp transitions for discontinuous problem, and its convergence reaches fifth order for smooth boundary condition.


---

## Body

## 1. Introduction

The fractional differential operations, such as the natural generalizations of the classical operators of integer order, have become more and more popular. Many phenomena in physics, engineering, biology, chemistry, and other sciences have been described successfully by fractional calculus and fractional differential equations [1]. Recently, several numerical methods to solve fractional differential equations have been developing rapidly. For example, the finite difference method [2–5], the finite element method [6–13], the discontinuous finite element method [14–20], the spectral method [21–24], the finite volume method [25, 26], etc. However, most of them are focusing on the fractional parabolic equation. It seems that there are few works regarding the numerical methods for fractional hyperbolic equation, where the physical phenomena often refer to vibration and oscillation of the complex material or complex progress. Compared to the classical hyperbolic equation, the solutions of fractional hyperbolic equation may also develop discontinuously even when the initial condition is smooth. It is impossible to deduce the analytical solutions in most cases, and one thus has to rely on numerical methods. So, it is very meaningful to study the numerical methods for fractional hyperbolic equation.As a typical representative of hyperbolic equation, the classical neutron transport equation studies the motions and interactions of neutrons with materials in a nuclear reactor. Besides, the transport equation also plays an important role within the transport theory. Compared to the classical neutron transport equation, the fractional neutron transport equation denotes the anomalous transport processes [27], which was studied by many researchers. Das and Biswas [28] developed a constitutive neutron diffusion equation for describing the neutron flux profile by using the concept of fractional divergence. Vyawahare et al. [29] studied a fractional-order model of a nuclear reactor, which predicts subdiffusive behaviour for long time. Vyawahare and Nataraj [30] dealt with a fractional-order model of neutron transport process and pointed out that it is a more faithful and realistic representation of neutron movements than the classical model.Weighted essentially nonoscillatory (WENO) schemes [31–34] are the effective numerical methods to solve the hyperbolic equation, specially, the fifth-order WENO (WENO5), which is a robust scheme for solving the classical neutron transport equation. Naturally, applying the WENO schemes to solve the fractional differential equation numerically is an interesting work. However, it seems that there are few works about the WENO schemes used to approximate the fractional differential equation except [3]. In [3], Deng et al. designed a WENO6 scheme for the fractional diffusion equation in the sense of Caputo derivative.In this work, we will design the generalized WENO5 schemes to solve the following spatial fractional neutron transport equation:(1)∂u∂t+λ∂αu∂xα=fx,t,where λ is a given constant coefficient and ∂αu/∂xα represents the Riesz derivative given by(2)∂αu∂xα=CαΓ1−αddx∫axx−τ−αuτ,tdτ−∫xbτ−x−αuτ,tdτ,in which Cα=−1/2cosπ/2α,0<α<1.The basic idea of numerical methods in this paper is to split theα0<α<1 order fractional derivative into a fractional integral and a classical first derivative. Next, we apply the Gauss–Jacobi quadrature to approximate the fractional integral and WENO5 to approximate the classical first derivative separately. Then, the three-stage, third-order SSP (3,3) explicit Runge–Kutta (ERK) method or the third-order TVD ERK method [35] is selected as the discrete time method to solve the ODE systems. The numerical experiments confirm fifth-order convergence order and show the maintenance of sharp discontinuous transitions. Furthermore, the linear stability of various ERK methods together with generalized WENO5 to approximate the fractional neutron transport equation is also discussed. From the linear stability analyses, we get that the combinations of the forward Euler scheme with WENO5 and any two-stage, second-order ERK schemes with WENO5 are linearly unstable in the fractional neutron transport equation. And, the combination of any three-stage, third-order ERK method with WENO5 for the fractional neutron transport equation is linearly stable. All of these theoretical results are demonstrated by numerical experiments.The structure of this paper is as follows. InSection 2, the generalized WENO5 schemes of the spatial fractional neutron transport equation are designed, including the Gauss quadrature formula relative to Jacobi weights, to compute the fractional integral terms. Then, the linear stability analyses of the designed schemes are discussed in Section 3. Finally, the numerical examples are presented to verify the convergence and stability of the designed schemes in Section 4.


## 2. The Generalized WENO5 Schemes

To achieve our method, we need to introduce the fractional integrals and fractional derivatives.Firstly, we define theαth order left- and right-fractional integrals of function ux as follows:(3)Ja,xαux=1Γα∫axusx−s1−αds,Jx,bαux=1Γα∫xbuss−x1−αds,where α>0. Next, we give the αth order left- and right-Caputo derivatives:(4)Da,xαCux=Ja,xn−αunx,Dx,bαCux=−1nJx,bn−αunx,in which n−1<α<n∈Z+.According to [36], when 0<α<1, we can easily derive the following relations between two fractional derivatives:(5)∂αu∂xα=CαDa,xαCux+Dx,bαCux+uax−a−α−ubb−x−αΓ1−α.Therefore, equation (1) can be rewritten as(6)∂ux,t∂t+λCαDa,xαC+Dx,bαCux,t+λCαua,tx−a−α−ub,tb−x−αΓ1−α=fx,t.By using the definition of the Caputo derivative, we have(7)∂ux,t∂t+λJa,x1−αvx,t−Jx,b1−αvx,t+λCαua,tx−a−α−ub,tb−x−αΓ1−α=fx,t,where vx,t=∂/∂xux,t.In order to get a high-order accuracy of the numerical solution, we apply Gauss-type quadrature formula relative to Jacobi weights to compute the fractional integral terms in (7). Therefore, we transfer the integral intervals a,x and x,b to the standard interval −1,1, respectively:(8)∫axx−ξ−αvξ,tdξ=x−a21−α∫−111−η−αva+x2+x−a2η,tdη,(9)∫xbξ−x−αvξ,tdξ=b−x21−α∫−111+η−αvb+x2+b−x2η,tdη.Now, we denoteh as the spatial step size, h=b−a/N, where N is the numbers of subintervals divided in the spatial directions. Then, by using N+1 points Gauss–Jacobi quadrature and the Lagrange interpolation, the left- and right-fractional integrals in the right-hand side of (8) and (9) at the mesh points x=xj can be, respectively, approximated by(10)∫axjxj−ξ−αvξ,tdξ∼xj−a21−α∑j=0Nvxj,t∑k=0NFja+xj2+xj−a2ηkωk,∫xjbξ−xj−αvξ,tdξ∼b−xj21−α∑j=0Nvxj,t∑k=0NFjb+xj2+b−xj2ηkω˜k,where ηkk=0N is the Gauss nodes, ωkk=0N and ω˜kk=0N are Jacobi weights, and Fj is the Lagrange interpolation basis function associated with the points xjj=0N, which can be replaced by any other kind of interpolation basis functions.In this paper, we denoteuxj,tn and ujn as the exact solution and the numerical solution at mesh point xj,tn, respectively. For the intermediate function vxj,t, we use WENO5 to do the discretization as follows:(11)vxj,t=∂u∂xxj,t=1hu^j+1/2−u^j−1/2+Oh5,j=0,1,…,N,where the term u^j+1/2 is the numerical flux and can be split into positive and negative parts:(12)u^j+1/2=u^j+1/2++u^j+1/2−.On stencil,Srj=Ij−r,…,Ij+s, r+s+1=3,r=0,1,2.Foru^j+1/2+ and u^j+1/2−, we define them on each stencil Srj as follows:(13)u^j+1/2+=ω0+13uj−2−76uj−1+116uj+ω1+−16uj−1+56uj+13uj+1+ω2+13uj+56uj+1−16uj+2,(14)u^j+1/2−=ω0−13uj+3−76uj+2+116uj+1+ω1−−16uj+2+56uj+1+13uj+ω2−13uj+1+56u1−16uj−1,where, for parameters ωk±, stencil Srj, and the other details of the WENO5 method, one can refer to [31, 34, 37].Next, we will use the forward Euler method, SSP (3,3) schemes, and some other Runge–Kutta schemes to solve the ODE system in the temporal direction after doing the spatial discretization.


## 3. Linear Stability Analysis

In this section, we will discuss the linear stability properties of the generalized WENO5 method when coupled with different kinds of ERK methods for solving the fractional neutron transport equation (1). Here, we define CFL number as σ=λΔt/h and only consider the corresponding homogeneous problem:(15)∂ux,t∂t+∂αu∂xαx,t=0,under homogeneous boundary conditions ua,t=ub,t=0.In the following, we will show that the combination of the generalized WENO5 and forward Euler scheme for the fractional neutron transport equation (15) is linearly unstable.Lemma 1 (see [37]).
When the forward Euler method is applied with the classical WENO5 scheme, the functionu^j+1/2+,n−u^j−1/2+,n/ujn, which is defined as zϕ, can be approximated by(16)zϕ=760ϕ6+Oϕ8+ıϕ+Oϕ7,where the Taylor expansions of sinϕ and cosϕ are used, ϕ is a small positive number, only related to the spatial step h, and ı is the imaginary unit, i.e., ı2=−1.Theorem 1.
The combination of the generalized WENO5 method and the forward Euler method for the fractional neutron transport equation (15) is linearly unstable.Proof.
Here, we only consider the Lax–Frichrichs flux splitting, which means that the negative part of the numerical flux is zero [37]. And, the same analysis and conclusions can be shifted to the other flux splitting.
According to the Fourier transform, we assume(17)uj=u^eıjϕ,where ϕ=ξh, and ξ is the spatial frequency. By using (13) and (17), we further conclude that(18)u^j+1/2+=ujω0+13e−2ıϕ−76e−ıϕ+116+ω1+−16e−ıϕ+56+13eıϕ+ω2+13+56eıϕ−16e2ıϕ.
When the forward Euler method is adopted, the fully discrete generalized WENO5 method becomes(19)ujn+1=ujn−σCαΓ1−α∑j=0NLuj−3n,uj−2n,…,uj+2nxj−a21−α∑k=0NFja+xj2+xj−a2ηkωk−b−xj21−α∑k=0NFjb+xj2+b−xj2ηkω˜k,where Luj−3n,uj−2n,…,uj+2n=u^j+1/2+,n−u^j−1/2+,n, which can be defined as zϕujn.
Then, equation (15) becomes(20)ujn+1=ujn−σCαΓ1−α∑j=0Nzϕujnxj−a21−α∑k=0NFja+xj2+xj−a2ηkωk−b−xj21−α∑k=0NFjb+xj2+b−xj2ηkω˜k=ujn1−σCαΓ1−α∑j=0Nzϕxj−a21−α∑k=0NFja+xj2+xj−a2ηkωk−b−xj21−α∑k=0NFjb+xj2+b−xj2ηkω˜k.
Define the amplification factor:(21)gσz=1−σCαΓ1−α∑j=0Nzϕxj−a21−α∑k=0NFja+xj2+xj−a2ηkωk−b−xj21−α∑k=0NFjb+xj2+b−xj2ηkω˜k.
By Lemma1, we get the amplification factor as follows:(22)gσz=1−σCαΓ1−αOϕ6+ιϕ+Oϕ7∑j=0Nxj−a21−α∑k=0NFja+xj2+xj−a2ηkωk−b−xj21−α∑k=0NFjb+xj2+b−xj2ηkω˜k=1−σOϕ6+ιϕ+Oϕ7∑j=0N∑k=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k,where Cα/Γ1−αxj−a/21−α and −Cα/Γ1−αb−xj/21−α are denoted by cj and dj correspondingly.
According to a simple calculation, we can get that(23)g2=1+σ2∑j=0N∑k=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k2ϕ2+Oϕ6>1,∀σ>0.
The proof is complete.
Next, we consider the followings−stage ERK method to replace the forward Euler method:(24)yn1=yn,ynk=yn+Δt∑i=1k−1ak,iftn+ciΔt,yni,k=2,…,s,yn+1=yn+Δt∑i=1sbiftn+ciΔt,yni,where the coefficients satisfy ck=∑j=1k−1ak,j, for k=1,…,s, and ∑i=1sbi=1.Lemma 2 (see [38, 39]).
The amplification factor of thes−stage, p−order ERK method is(25)gz^=1+∑l=1p1l!z^l+∑l=p+1s∑j=lsbj∑j>m1>⋯>ml−1≤1aj,m1am1,m2⋯aml−2,ml−1z^l=1+∑l=1p1l!z^l+∑l=p+1sz^lbTAl−1e,where b=b1,b2,…,bs,A=aij1≤i,j≤s are the coefficients of the s−stage, p−order ERK method, e=1,1,…,1, and z^=−σz, where z is caused by the character of the spatial operator.Remark that, from Theorem1, we can see that the operator z^ should be adjusted to −σ∑j=0Nzϕ∑k=0NcjFja+xj/2+xj−a/2ηkωk+djFjb+xj/2+b−xj/2ηkω˜k in the fractional neutron transport equation, which is caused by the left and right spatial fractional integral terms.Theorem 2.
The combination of the generalized WENO5 method with any two-stage, second-order ERK method for the fractional neutron transport equation (15) is linearly unstable.Proof.
According to Lemmas1 and 2, the amplification factor can be written as(26)gz^=1+z^+12z^2=1−σ∑j=0Nzϕ∑k=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k+12σ2∑j=0Nzϕ∑k=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k2=1−σ760ϕ6+Oϕ8+1ϕ+Oϕ7∑k=0N∑j=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k+12σ2760ϕ6+Oϕ8+1ϕ+Oϕ72∑k=0N∑j=0NcjFja+xj2+xj−a2ηkωk+djFjb+xj2+b−xj2ηkω˜k2=1−κσ760ϕ6+Oϕ8+1ϕ+Oϕ7+12κ2σ2760ϕ6+Oϕ8+1ϕ+Oϕ72,where the constant ∑k=0N∑j=0NcjFja+xj/2+xj−a/2ηkωk+djFjb+xj/2+b−xj/2ηkω˜k is denoted by κ.
Then, we have that(27)g2=1−12κ2σ2ϕ2+Oϕ62+−κσϕ+Oϕ72=1+14κ4σ4ϕ4+Oϕ6>1,∀σ>0.
The proof is finished.Theorem 3.
The combination of the generalized WENO5 method with any three-stage, three-order ERK method for the fractional neutron transport equation (15) is linearly stable.Proof.
According to Lemmas1 and 2 and the proof of Theorem 2, we can get the amplification factor as follows:(28)gz^=1+z^+12z^2+16z^3=1−κσ760ϕ6+Oϕ8+ıϕ+Oϕ7+12κ2σ2760ϕ6+Oϕ8+ıϕ+Oϕ72−16κ3σ3760ϕ6+Oϕ8+ıϕ+Oϕ73.
Then, we have that(29)g2=1−12κ2σ2ϕ2−760κσϕ6+Oϕ82+−κσϕ+16κ3σ3ϕ3+Oϕ72=1−112κ4σ4ϕ4+−760κσ+136κ6σ6ϕ6+Oϕ82<1,∀σ>0.
The proof is finished.
From the above proof, we can see that the conclusion of the linear stabilities for the combination of the generalized WENO5 method with the ERK method in the fractional neutron transport equation is the same as in the classical neutron transport equation. The reason is that the added fractional integral terms are approximated by a linear combination, in which the weights are bounded in the Gauss–Jacobi quadrature formulas. Therefore, all of the conclusions about the linear stabilities for the classical hyperbolic conservation law [37] can be shifted to problem (1), where the proofs are omitted.


## 4. Numerical Examples and Results

In this section, we will design three numerical examples to examine the effectiveness of the numerical algorithms. The first example with periodic boundary conditions is designed to verify the linear stability of the proposed schemes. The second example with two points’ boundary conditions is computed to verify the convergence order. And, the third example is designed for numerical solving piecewise smooth solution and discontinuous initial condition.Example 1.
Consider the following fractional neutron transport equation:(30)∂ux,t∂t+D0,xαRLux,t=fx,t,x∈0,1,where ux,t=e−tsin4πx−2sin2πx is chosen as the exact solution.
Firstly, we use the forward Euler method with the WENO5 scheme. Next, we use the two-stage, second-order SSP ERK method with the WENO5 scheme. We compute the two cases forα=0.9andσ=0.3,0.5,0.8,1.1, which are shown in Figures 1 and 2 separately. According to the theoretical analyses, spurious oscillations due to linear instabilities are present.
Then, we use the SSP (3,3) ERK method with the WENO5 scheme.L∞ error and convergence order are shown in Table 1. From Table 1, we can see that the numerical scheme performs well and shows the fifth-order convergence accuracy.Figure 1
Absolute error of Example1 with the forward Euler method and WENO5 scheme for different CFL numbers at T=0.1,N=2560,andα=0.9. (a) σ=0.3, (b) σ=0.5, (c) σ=0.8, and (d) σ=1.1.
(a)(b)(c)(d)Figure 2
Absolute error of Example1 with the two-stage, second-order ERK method and WENO5 scheme for different CFL numbers at T=0.1,N=10240,andα=0.9. (a) σ=0.3, (b) σ=0.5, (c) σ=0.8, and (d) σ=1.1.
(a)(b)(c)(d)Table 1
Errors and convergence rates of Example1 with the third-order TVD ERK method and WENO5 scheme for α=0.9 and σ=0.4 at T=0.1.
NL∞cv.rate404.6260e − 003—801.0090e − 0032.19681601.8011e − 0055.80793203.4944e − 0075.6877Example 2.
Consider the following fractional neutron transport equation:(31)∂ux,t∂t+D0,xαRLux,t=fx,t,x∈0,1,with the initial condition, two points’ boundary conditions,(32)ux,0=x41−x3,u0,t=u1,t=0,and the source term is(33)fx,t=e−t−x4+3x5−3x5+x7+ax4−α+bx5−α+cx6−α+dx7−α,where a=−24/Γ5−α,b=360/Γ6−α,c=−2160/Γ7−α,andd=5040/Γ8−α. Then, the exact solution is ux,t=e−tx41−x3.
Applying the SSP (3,3) method with the WENO5 scheme, we getL∞ error and convergence order in Table 2. And the fifth-order convergence of the combination of the SSP (3,3) method with the WENO5 scheme is still observed.Table 2
TheL∞ errors and convergence rates of Example 2 with the third-order TVD ERK method and WENO5 scheme for α=0.9 and σ=0.4 at T=0.1.
NL∞cv.rate402.3533e − 006—801.0556e − 0074.47861607.5624e − 0093.80313202.6494e − 0104.8351Example 3.
Suppose thata,b=−10,10. Consider the fractional neutron transport equation as follows:(34)∂ux,t∂t+∂αu∂xαx,t=0,x∈−10,10,with the boundary conditions u−10,t=1,u10,t=0, and the initial condition(35)ux,0=1,x≤0,0,x>0.
In the simulation of Figure3, we apply the SSP (3, 3) ERK method with the WENO5 scheme and take N=200andσ=0.1. Figure 3 shows numerical solutions with different orders of fractional derivative, α=0.2,0.4,0.6,0.8 at T=0.1. It can be seen that sharp transitions are retained in the numerical simulations for different α.
Figure4 shows the numerical solutions and their absolute errors for α=0.9,N=100 and 200, and σ=0.1 at T=1. It can be seen that the bigger the N value, the sharper the transitions in the simulations.
Figure5 shows the numerical solutions when the parameters are chosen as T=0.1,α=0.9,N=200, and σ=0.1,0.3,0.5,0.7 separately.
Figure6 shows the numerical solutions when the parameters are chosen as α=0.9andN=200, at T=1,T=15,andT=30 separately.Figure 3
Numerical solutions of Example3 with N=200andσ=0.1 for different orders of fractional derivatives at T=0.1. (a) α=0.2, (b) α=0.4, (c) α=0.6, and (d) α=0.8.
(a)(b)(c)(d)Figure 4
The numerical solutions and their absolute errors of Example3 with α=0.9andσ=0.1 for different N at T=1. (a) N = 100, (b) N = 100, (c) N = 200, and (d) N = 200.
(a)(b)(c)(d)Figure 5
The numerical solutions of Example3 with α=0.9andN=200 and for different CFL numbers at T=0.1. (a) σ=0.1, (b) σ=0.3, (c) σ=0.5, and (d) σ=0.7.
(a)(b)(c)(d)Figure 6
The numerical solutions of Example3 with α=0.9andN=200 at T=1,T=15,andT=30.


## 5. Conclusion

This paper develops the generalized WENO5 methods for numerically solving the spatial fractional neutron transport equation with smooth solutions or piecewise smooth solutions, which is a new attempt in the work of numerical solution for fractional differential equation. By using the classical WENO5 schemes and the Gauss–Jacobi quadrature, we perform the spatial discretization. Then, we get the nonlinear ODE systems, which are computed by the SSP (3,3) Runge–Kutta method and some other Runge–Kutta method. The linear stabilities of the numerical schemes are proved. And, the numerical results are present to verify the effectiveness and reliability of the numerical schemes.Generally speaking, the WENO scheme is widely used to solve the problem that the solution of the equation is not very smooth and can maintain a high-order convergence rate as well. Therefore, for the upcoming work, we will mainly discuss the discontinuous problems and the nonlinear problems with the WENO method in the future. Maybe, by using the nonuniform mesh method or meshless method together with the WENO scheme for the fractional problems is a good choice.



---

*Source: 6676640-2021-03-15.xml*