High-order solvers for space-fractional differential equations with Riesz derivative

This paper proposes the computational approach for fractional-in-space reaction-diffusion equation, which is obtained by replacing the space second-order derivative in classical reaction-diffusion equation with the Riesz fractional derivative of order \begin{document}$ α $\end{document} in \begin{document}$ (0, 2] $\end{document} . The proposed numerical scheme for space fractional reaction-diffusion equations is based on the finite difference and Fourier spectral approximation methods. The paper utilizes a range of higher-order time stepping solvers which exhibit third-order accuracy in the time domain and spectral accuracy in the spatial domain to solve some fractional-in-space reaction-diffusion equations. The numerical experiment shows that the third-order ETD3RK scheme outshines its third-order counterparts, taking into account the computational time and accuracy. Applicability of the proposed methods is further tested with a higher dimensional system. Numerical simulation results show that pattern formation process in the classical sense is the same as in fractional scenarios.


1.
Introduction. The concept of fractional differential equations has been considered as a useful tool in the modeling of many physical phenomena. In the study of groundwater resources, non-integer order models give the most meaningful and adequate description of contaminant and chemical transport in heterogeneous aquifers [5,12,45,53].
Fractional differential equations have also been applied in transport dynamics to describe the complex systems that are governed by non-exponential and anomalous diffusion patterns [16]. Nowadays, Fractional differential equations are frequently applied in all application areas of engineering and science, such as: aerodynamics, electrical systems, signal and image processing, economics analysis of the stock pricing, finance, biomedical engineering, polymer rheology, electrodynamics, control theory, biophysics, viscoelasticity theory and blood flow phenomena [15,35,41,42,48,[63][64][65][66].
Moreover, fractional calculus models provide a relatively simple way of describing the physical and electrical properties of complex, heterogeneous, and composite biomaterials [43]. Non-integer order calculus can yield a concise model for the description of the dynamic events that arise in biological phenomena. In this paper, we shall consider the general space fractional (non-integer) reaction-diffusion equation u i (x, 0) = ψ 0 (x), x ∈ (a, b), u i (a, t) = u i (b, t) = 0, t ∈ (0, T ], i = 1, 2, . . . , n, where α ∈ (0, 1) for subdiffusive or α ∈ (1, 2] for the superdiffusive processes, and R D α x is the Riesz fractional derivative operator defined by with c α = 1 2 cos(απ/2) is satisfied for all u, and u i (x, t) is regarded as either species concentration or population size and n the total number of species (components) taken into account in the model. The term ∆ α is called the laplacian operator with order α in the domain 0 < α ≤ 2 given by the Riesz representation, which is classified into subdiffusive for 0 < α < 1 and superdiffusive for 1 < α ≤ 2, and the case α = 2 corresponds to the integer order problems [54], such systems arise in numerous application in chemical kinetics and populations dynamics. The term d i > 0 is the diffusivity of the ith spaecies, and f = f i (u i (x, t)) as nonlinear function that accounts for the source term, with the homogeneous Neumann (u i (a) = 0 u i (b) = 0) and homogeneous Dirichlet boundary conditions (u i (a) = Λ 1 , u i (b) = Λ 2 ), where Λ 1 , Λ 2 are constants. Qualitative behaviour of system (1) for i = 1, 2 has been reported in [60] for α = 2. Generally speaking, a coupled multi-species fractional reaction-diffusion models n (for n ≥ 1) is an equation of the form ∂u i (ψ) ∂t = d iR D α x u i (ψ) + f i (u 1 , . . . , u n ), i = 1, . . . , n, subject to zero-flux boundary conditions and the initial conditions u i (ψ, 0) = u i0 (ψ), i = 1, 2, . . . , n where d i and u i are often regarded as the diffusivity and n−dimensional concentration or population density of the ith species respectively, ψ is the position in space and f i (u i ) which account for the nonlinear interactions and local kinetics that characterized the various growth and interaction features of the coupled system, with f i (u 1 , u 2 , . . . , u n ) denoting the overall growth rate for the ith species. Analysis of the n− species system (3) still remain tactical in literature, the stability of the steady states is obtained in exactly the same manner as given in [60,62] by linearizing about the steady states u i =û i , i = 1, . . . , n. In absence of diffusion and at α = 2 for the sake of simplicity, the solution u i =û i of the system (3) is said to be Turing unstable if it is locally stable. If a steady state is unstable then the solution u may evolve into another steady state or into a stable oscillatory structure like a limit cycle [52,61] or spatiotemporal chaos patterns [55,62]. Hence We then examine the eigenvalues λ of the stability or community matrix defined by where the elements b i,j determine the effects of the j−species on the i−species near the steady state. By adopting the Routh-Hurwitz conditions [68], the necessary and sufficient condition for the eigenvalues of B is that, solutions of |B − λI| = 0, for Reλ > 0. For instance, if u i is the perturbation to the steady stateÛ i , the equation for u i yields If b i,j > 0 then U j improves U i s growth but if b i,j < 0 it endangers (extinction) its growth. If both b i,j > 0 and b j,i > 0, there exists a mutual interaction as both U i and U j enhance each other's growth. We say that the species are into competition if both b i,j < 0 and b j,i < 0. For detail analysis on multi-species models, we refer readers to some classical books [6,21,36,52]. Even though the discussion on derivatives of non-integer order could be dated back as far as the development of the classical theory of integer order differential calculus, many researchers have investigated fractional calculus concepts in various ways as well as the numerical solutions of fractional-order differential equations using different numerical techniques. A lot of numerical methods have been used to discretize the fractional Laplacian operator in (1) [13,17,24,62].
Several numerical methods have been used to discretize the Laplacian operator of equation (1) [3,4,7,10,11,22,62]. However, knowing well that the fractional differential operator is nonlocal, which has posed some numerical and computational challenges [14], other than the issue of stiffness encountered in the classical reactiondiffusion equations.
The Riesz fractional derivative has been adopted in analysis of both space and time fractional diffusion problems [44,67], where the solution in terms of they expressed the solutions in terms of Fox H−function is expressed. It is well known that the space fractional diffusion equation with fractional Riesz space derivative yields same results as those obtained from the continuous time random walk theory for Lévy flights [49,50]. A numerical method for solving fractional diffusion problem with Hilfer composite fractional time derivative and Riesz space fractional derivative has been reported in [67]. Luchko et al. [40] have applied the quantum fractional Riesz derivative to study the Schrödinger equation for a free particle.
In this paper, we examine the solution of a fractional reaction-diffusion system of the form (1) in which the space fractional derivative is replaced by the Riesz derivative x D α θ , and the resulting systems will be advanced in time with the family of third order schemes based on implicit-explicit and the exponential time differencing methods. Many time-dependent space fractional partial differential equations combine low-order nonlinear terms with higher-order fractional or nonlocal terms, such as the space fractional Fisher, Allen-Cahn, Schrödinger, Gray-Scott, Schnakenberg and Swift Hohenberg equations [56,62]. It is desirable to use highorder approximations for such problems in both space and time to obtain accurate numerical results. Most numerical methods that have been used are restricted to second order in time [22,46,47]. In this paper, we discretize the space fractional derivative with the Fourier spectral methods [57,58,62]. This approach is known to give full diagonal representation of the fractional operator and achieves spectral convergence regardless of the value of fraction power used in the computation [22,62]. For the temporal discretization, we formulate a range of fractional third order implicit-explicit (IMEX) schemes, based on the linear multistep, predictorcorrector and Runge-Kutta schemes, as well as the third order exponential time differencing Runge-Kutta method [25]. To the best of our knowledge, these third order schemes have never been used to study the space fractional reaction-diffusion problems in high dimension. The reason for using higher-order schemes is that we expect our spectral method to provide exponential accuracy, which is usually lack by lower order temporal integration methods. These methods are extended to solve high-dimensional models.
The paper is organized as follows. In Section 2, we give some of the useful definitions based on one-variable and two-variable system, and properties of various fractional derivative operators. Different approximation methods in time and space are discussed in Section 3. Efficiency and effectiveness of the proposed schemed are justified with some numerical experiments in Section 4. We conclude the paper with Section 5.

Preliminaries and definitions.
We begin this section with preliminary definitions of fractional derivatives for one and several variable, see for instance, [1,35,51,63].
Definition 2.1. Let f be an analytic function on stripe C ∈ Ω which covers the real axis. If f is α−differentiable at the point x 0 ∈ R, there exists L ∈ R and a function δ : R → R such as lim with N α f (x 0 ) is called the Nishimoto's fractional derivative of order α in (0, 1] [1].
Definition 2.2. Let f be an analytic function on stripe C n ∈ Ω which covers R n . For any α−differentiable function u : R n → R p at a point x = (x 1 , x 2 , . . . , x n ) ∈ R n , n, p ∈ N, there exists a homogeneous function of order α−degree L : R n → R p such as for ∈ (R + ) n , if α = 1 2n+1 ∈ (0, 1], then for a given x, the map ∈ x : R n → R p verifies lim →0 ∈ x ( α ) = 0 and α = ( α 1 , α 2 , . . . , α n ) for = ( 1 , 2 , . . . , n ). The function L is known as the α−differential of u at x, denoted as d α u x . The α−th transform of u at point x = is given by for j = 1, 2, . . . , p. Assume E is an open set in R. If u is α−differentiable at x ∈ E, we conclude that the function u is α−differentiable in E.

Definition 2.3. (High-order PDE) Assume Ω is an open set in
x − j at the point σ. hence, the real PD of the 2α−th of the function u is given by G, while that of the pα−th of the function u is denoted by In what follows, we report some useful definitions of fractional derivative. The Liouville-Weyl (LW) fractional integrals of order 0 < α < 1 are defined by The LW fractional derivatives of order 0 < α < 1 are given as the left inverse operators that correspond to the LW fractional integrals Alternatively, we can write the definitions in (8) as which can be derived simply by A specific choice of linear combination of the LW fractional integrals leads to the Riesz fractional integral I α R : Therefore, the Riesz fractional derivative is given as (12) On the other hand, Feller [28,29] proposed the generalization of the Riesz fractional derivative in the form The Feller fractional derivative is given by which tallies with the Riesz fractional definition in (11).
Again, setting θ = 1 yields another special case Hence, as a linear combination of D α 1 and D α 0 one can write the Feller fractional derivative as where which reads The parameter θ here is interpreted as a rotational parameter for the Feller fractional derivative, instead of a skewness parameter.
In the spirit of Feller [28], we can write the Riesz-Feller fractional-in-space derivative of order α and skewness θ conveniently in the form where stands for the Fourier transform and its inverse respectively. The function ω θ α (κ) is given by When θ = 0, we get a symmetric operator w.r.t. x that can be defined as For 0 < α ≤ 2 and |θ| ≤ min{α, 2 − α}, the Riesz-Feller derivative can be shown in the x domain with integral representation When θ = 0, the Riesz-Feller fractional derivative reduces to the Riesz fractional derivative of order α for 1 < λ ≤ 2 bounded by the analytic continuation in the domain 0 < α ≤ 2 for α = 1. In what follows, and for the remainder part of this paper, we use the Riesz derivative operator.
3. Algorithm for the Riesz space-fractional derivative. The Riesz spacefractional derivative with superdiffusive order α ∈ (1, 2) is given as [35,38,70] where R D α x , RL D α a,x and RL D α x,b are the respective Riesz derivative, left and right Riemann-Liouville derivatives. We use the mesh points x n = a+nh, n = 0, 1, 2, . . . , N , and h = (b − a)/N is the step size in spatial direction.

3.2.
Fourier spectral spatial discretization method. This discretization technique presented with the FDM above is a low-order accuracy type. Here, we give high-order discretization of the initial-boundary value problem in (1) via the Fourier spectral method (FSM).
Then for any u ∈ U α , the Laplacian operator (−∆) α is given by where Φ i and Λ i in the computation depend on the specified or given boundary conditions: i The homogeneous Neumann boundary condition ii The homogeneous Dirichlet boundary condition It is obvious from the Definition 3.1 that one can approximate the exact solution where N > 0 is an integer. Moreover, by combining with equation (30), the i−th Fourier mode of the fractional reaction-diffusion equation wheref i is the i−th Fourier coefficient of the local kinetic or reaction term. Spectral methods have been applied recently in various ways to different fractional order problems, see for example, [18-20, 22, 56-58, 62] and references therein.

Time-stepping integrators for the fractional reaction-diffusion equation.
In this section, we report briefly various time integrators of third-order schemes formulated to advance the resulting nonlinear system of ordinary differential equations (ODE) in time. To start with and for the sake of convenience, we write (1) in a general compact form where F represents a generic nonlinear source term, and A = ∆ α is the Laplace operator.
3.3.1. Exponential time differencing Runge-Kutta method. The idea of solving ODE system that are stiff was introduced by Hairer and wanner [32], which has since gained a lot of grounds and attract some scholars in the study of computational dynamics, which include for example, linearly implicit schemes [2], integrating factors (IF), projection method, time-splitting schemes [69], multiscale methods [26,30] and the exponential time differencing (ETD) [25]. The ETD scheme, first studied by Cox and Matthews [25], is one of the various schemes which have developed for solving stiff differential equations, see [32] for details discussion. Our interest in the present paper is the modified ETD schemes designed for the diffusive type equations [33], which we are formulating to handle general space fractional reaction-diffusion systems.
Discretizing equation (32) in the spatial variables, for example, using either the spectral approximation or finite difference techniques, results to a system of ODEs The derivation of the s-step exponential time differencing (ETD) schemes is followed directly from [25]. We begin by multiplying the ODE form of (33) by the integrating factor e −At , and then integrate over a single time step from t = t n and t = t n+1 = t n + h to have Equation (34) is known to be exact, various ETD schemes arise from the approximation of its integral. A third-order exponential time differencing Runge-Kutta scheme denoted as the ETD3RK is given by a n = u n e Ah/2 + e Ah/2 − 1 F (u n , t n )/A, b n = u n e Ah + (e Ah − 1)[2F (a n , t n + h/2) − F (u n , t n )]/A, +4F (a n , t n + h/2)[2 + hA + e Ah (−2 + hA)] The values of u at t n + h/2 and t n + h are approximated by the terms a n and b n respectively. The analytical studies of the ETD and the modified ETD methods can be found in [25,27]. Its stability region is presented in Figure 1 (a).

3.3.2.
Implicit-explicit predictor-corrector scheme. Implicit-explicit (IMEX) schemes are another family of well studied schemes with an established history for solving stiff partial differential equation in literature.
The IMEX schemes consist of using an explicit multistep method, for instance, the second order Adams-Bashforth method, to integrate the nonlinear part of the problem and an implicit scheme, for instance, the second order Adams-Moulton formula, to integrate the linear part. There are other kinds of formulations, for example, for derivations based on Runge-Kutta rather than Adams-Bashforth schemes, see [8,23,31]. In this paper, we use an IMEX scheme that belongs to a class of implicit-explicit predictor-corrector methods IMEX-PC, which has been successfully applied to tackle stiff PDEs. The main objective here is to split the known multistep IMEX method into predictor-corrector schemes. A family of higher order IMEX-PC was constructed by [39] based on the strategy obtained in [9].
By adopting the techniques of Ascher et al. [8,9], the general s−step IMEX scheme when applied to the ODE (33) can be written as The split form of (36) results to the following IMEX predictor-corrector schemes (−a j u n+j + kb j Au n+j + kc j F (u n+j , t n+j )), predictor, (−a j u n+j + kb j Au n+j + kb j F (u n+j , t n+j )) +kb s g(ū n+s , t n+s ), corrector.
The third-order IMEX predictor-corrector scheme used in this paper is denoted by IMEX3PC. The IMEX-PC(3, m) is drawn from a family of 3-step, three parameters (µ, ψ, η) schemes of order m given as: The choice (µ, ψ, η) = (1, 0, 0) leads to IMEX3PC scheme of step 3 and order 3. It should be noted that the order of the predictor-corrector method is determined by the predictor and the corrector. The error constant of the proposed scheme is dependent on u, b s , F (u, t) and its derivatives. This assertion is validated with the following theorem. This theorem shows that the order of the predictor is not necesaarily the same as that of the corrector. The proof of similar theorem is readily given in [39]. Following the stability analysis of Li et al. [39], we plot the stability region of IMEX3PC when applied to (33) in Figure 1 (b).

3.3.3.
The third-order implicit-explicit Runge-Kutta scheme. Implicit-explicit Runge-Kutta (IMEXRK) schemes are another family of linear multistep methods that have been used by various authors in different forms. It should be mentioned that among the IMEXRK schemes developed for (33), for instance, Ascher et al. [8] constructed a family of L-stable two-,three-stages diagonally implicit Runge-Kutta (DIRK) and four-stage, third-order combination schemes whose formulations were based on implicit-explicit Runge-Kutta methods for solving convection-diffusion problem. In another development, Kennedy and Carpenter [34] constructed a family of higher-order, L-stable explicit and singly diagonally implicit Runge-Kutta using IMEX schemes for solving one-dimensional convection-diffusion-reaction equations. Years later, Koto [37] constructs IMEX Runge-Kutta schemes for reaction-diffusion equations with an established convergence and stability. In the spirit of [60], an explicit Runge-Kutta (ERK) scheme is used to solve the non-stiff part F(u) and a diagonally implicit Runge-Kutta (DIRK) scheme is employed to solve the stiff part L(u) of equation (33) to overcome the issue of stability restriction inherent in the explicit method.
The s−stage implicit-explicit Runge-Kutta (IMEXRK) scheme from t n−1 to t n can be given in general form as: where 1 ≤ l ≤ s, σ τ = t n − t n−1 , and g 1 , g 2 treated the implicit and explicit parts respectively.

Numerical experiments.
Some numerical examples are considered in this section to demonstrate the applicability, and justify the effectiveness as well as the accuracy of the various proposed methods. The accuracy of the methods when applied to one-and two-dimensional fractional problems is measured by reporting the maximum norm errors where N is the number of points utilized for the discretization in one and two dimensions, u . and u(·, T ) are the numerical and exact approximations to the given problem respectively.

Example 1.
For the purpose of comparison, we first consider the fractional reaction-diffusion equation (1) in one-dimension, subject to homogeneous Dirichlet boundary conditions. We employ the exact analytical solution and its corresponding local kinetic in [22] in x ∈ (0, 1) given by The results in Table 1 shows the maximum norm error for the numerical solution of (1) using both the finite difference and spectral methods in conjunction with the IMEX3RK at different instances of α, with d = 0.5, N = 100, T = 1 and initial condition u(x, 0) = exp(−100(x − 0.5) 2 ). It is obvious that the error of the FSM is smaller to that of the FDM regardless of the value of α, though both methods have almost the same timing results of ≈ 0.17s. Figure 2 illustrates the convergence of the third-order schemes at t = 0.1 and t = 2 for α = 1.45. When simulation time is increased from 0.5 to 2.0, the ETD3RK has the better convergence compared to other two methods.
Similarly, we extend numerical experiment on problem (1) to two-dimensional case to actually ascertain which of the schemes performs better than the others. The exact and the local kinetic considered here are given by The two-dimensional results are presented in Table 2. However, we realized in both one-and two-dimensional experiments that the schemes give a better accuracy in sub-diffusion region 0 < α < 1 than the super-diffusion case 1 < α < 2 for increasing value of N . Due to the remarkable performance displayed by the ETD3RK method when applied to solve one-and two-dimensional problems over other methods considered in this paper, the remainder experiments will be carried out with only the ETD3RK.  In this subsection, we consider a fractional reaction-diffusion system that is characterized by the disributions of two chemical species u and v given in the form The origin (0,0) is the only spatially uniform steady state. For diffusion-driven (Turing space) instability to occur, the conditions; +β < 0, (β +1) > 0, Dβ + > 0 and (Dβ + ) 2 − 4D (β + 1) > 0 must be satisfied. For the numerical simulation in 2D, we choose a periodic boundary conditions ±L subject to the axisymmetric initial conditions used in [59] as: where r 2 = x 2 + y 2 . The choice of using axisymmetric initial conditions above allow us to examine the 2D computations in a nontrivial manner, since the system in the presence of diffusion is unstable. To give room for the waves to propagate, spatial dimensions ±L, L = 20 are chosen large enough. Numerical simulations reveal that the distributions of the two chemical species here are not almost the same, hence we display results for species u and v as seen in Figure 3 at some instances of α. The upper-row represent both plot and surface mesh results at α = 0.5 for both species, while the lowerrow typifies the distribution of the species when α = 1.55, see the figure caption for other parameters. Various complex and chaotic spatiotemporal structures are observed especially in the superdiffusive case, depending on the choice of the initial data and parameter values, see for instance, Figure 4 for mitotic patterns in the supperdiffusive scenarios. Applicability of the proposed numerical method is tested when applied to model (42) in three-dimensions. We set the initial conditions as: u(x, y, z, 0) = 1 − exp(−10(x − 0.5) 2 + (y − 0.5) 2 + (z − 0.5) 2 ), v(x, y, z, 0) = exp(−10(x − 0.5) 2 + 2(y − 0.5) 2 + (z − 0.5) 2 ), (44) subject to the homogeneous Dirichlet boundary conditions clamped at the ends of domain size ±L 3 , for L = 5. A keen look at the 3D results displayed in Figure 5 have shown that the evolution of both species in space is almost similar. Hence, we restrict our pattern analysis to only one distribution. A 3D surface plot illustrating the distribution patterns of v−species due caused by diffusion-driven instability is shown in Figure 6. The first and second columns correspond to subdiffussive (0 < α < 1) and superdiffusive (1 < α < 2) scenarios. Clearly, we can conclude based on the results in Figure 6 that as the value of fractional power α is increasing or tending to an integer in both sub-and superdiffussive regions, so also the patterns are changing from being irregular to regular shapes.
The spatiotemporal structures obtained here are based purely on the species spatial interactions in the presence of diffusion, and not by any underlying heterogeneity in the domain. Depending on the choice of parameters and initial data, it should be mentioned that formation of other regular and irregular structures are possible. Figure 6. Three dimensional results for system (42) at different instances of fractional power α, with τ 3 = 0.26 and final time t = 5. The first and second columns correspond to subdiffusive and superdiffusive cases. Other parameters are given in Figure 3 caption. 5. Conclusion. In this paper, we have considered three third-order time stepping schemes for solving a range of fractional-in-space reaction-diffusion equations. The Riesz derivatives are approximated with both fractional finite difference method and Fourier spectral methods. Numerical simulation results have shown that the Fourier spectral method gained a significant accuracy over its finite difference counterpart. In this work, we also make a numerical comparison of the third-order schemes (IMEX3PC, IMEX3RK, ETD3RK) which have been applied to solve phase field models and reaction-diffusion problems. The comparison here is based on the simulation of space fractional reaction-diffusion equations in one and two spatial dimension. Numerical results obtained in sub-diffusive region 0 < α < 1 are more accurate than those in the super-diffusive case 1 < α < 2. In higher dimensions, we realized that pattern formation process in the fractional scenario is as the same as those obtained in the classical case. When applying to integer and noninteger order problems, the schemes adapted in this work have been shown to be accurate and efficient. In the future, these schemes will be applied to more real-world fractional reaction-diffusion equations that are of recurring interests in science and technology.