Numerical analysis and pattern formation process for space-fractional superdiffusive systems

In this paper, we consider the numerical solution of fractional-in-space reaction-diffusion equation, which is obtained from the classical reaction-diffusion equation by replacing the second-order spatial derivative with a fractional derivative of order \begin{document}$ α∈(1, 2] $\end{document} . We adopt a class of second-order approximations, based on the weighted and shifted Grunwald difference operators in Riemann-Liouville sense to numerically simulate two multicomponent systems with fractional-order in higher dimensions. The efficiency and accuracy of the numerical schemes are justified by reporting the norm infinity and norm relative errors as well as their convergence. The complexity of the dynamics in the equation is theoretically discussed by conducting its local and global stability analysis and Numerical experiments are performed to back-up the theoretical claims.


1.
Introduction. The subject of fractional calculus is as old as classical calculus, which is based on a generalisation of integration and ordinary differentiation of arbitrary non-integer order. The idea of the fractional derivative has been an active subject of interest over the years, not only among the mathematicians but also among the engineers and physicists [56].
Due to the accuracy of integral and fractional differential equations in modelling various physical phenomena [46,61], the subject of fractional calculus has gain popularity and becomes the focus of many researchers. As a result, many researchers have shown their interest to study the properties of fractional calculus and provide robust and accurate analytical and numerical techniques for solving fractional differential equations, for example, the differential transform method [51], the finite element and finite difference methods [15,[33][34][35]43,66], the spectral collocation and tau methods [12][13][14]18,23,29,59], Chen [19] studied fractional diffusion equations by using the Kansa method, in which the MultiQuadrics and thin plate spline served as the radial basis function. Ray [58] used a two-step Adomian decomposition scheme for the analytic solution of space fractional diffusion equation. A multistep differential transform approach was introduced by Odibat and Momani [38] to address both chaotic and non-chaotic systems. Many other methods of solving the fractional derivative problems are classified in [27,[43][44][45]50].
A lot of definitions of fractional derivatives have been introduced among the suitable and commonly used ones are the definitions by Caputo, Coimbra, Davison and Essex, Erdely-Kober, Grünewald-Letnikov, Hadamard, Marchaud, Riesz, Weyl, and Riemann-Liouville [4,11,16,17,20,28], the most recent one is the Atangana and Baleanu derivative in Riemann-Liouville sense [5,8,9], which has sparked many scholars attention, see for instance [1,[6][7][8][9][24][25][26]. It is unfortunate that similar usages of different definitions have lead to different results [56]. For instance, Khan et al. [30] used the Caputo fractional derivative description and applied the homotopy analysis method to obtain both the analytical and numerical solutions of the nonlinear reaction-diffusion equations with time-fractional derivatives which they convert to a hierarchy of linear equations. Ashyralyev [3] obtained the formula for fractional difference derivative by using the connection between the fractional power with the fractional derivatives positive operators. Ding and Jiang [22] applied the spectral representation of the fractional Laplacian operator for the conversion of these equations into multi-term time fractional ordinary differential equations to analytically obtain the solution of the multi-term time-space fractional advection-diffusion equations on a finite domain with mixed boundary conditions. Atangana [4] adopts the concept of derivative of fractional variable order to generalise the time-fractional telegraph equation, which was numerically solved via the Crank-Nicholson scheme and confirm the stability and convergence of the method.
The system of equations described by the fractional derivative of order α is mostly tedious to solve analytically for all fractional orders bounded in the superdiffusive interval 1 < α ≤ 2, due to some fundamental challenges in numerically approximating the fractional derivatives. Hence, we seek an appropriate and efficient numerical technique to circumvent such difficulties. In this paper, we consider the time-dependent fractional-in-space reaction-diffusion system using a class of second-order approximations, known as the weighted and shifted Grünwald difference schemes in the sense of Riemann-Liouville fractional operators, subject to the Crank-Nicholson technique for the time discretization.
The aims of this paper are summarised as follows. We propose in Section 2, a class of discrete operators to approximate the Riemann-Liouville fractional derivatives with high order truncating errors. Further, by using the finite difference method whose formulation were based on the weighted and shifted Grünwald-Letnikov formulae to numerically solve the general fractional diffusion equation in one and two dimensions, their effectiveness and accuracy are theoretically proved and numerically verified by reporting the maximum (norm infinity) and L 2 (norm relative) errors. We analyse the linearized steady state and also study the global existence and bifurcation of non-constant positive steady states of the multicomponent model in Section 3. Some numerical results are presented in one and two dimensions in Section 4 to justify the accuracy and efficiency of the methods. Adequate Concluding remarks are drawn in the last Section.
2. Fractional order derivatives and numerical method of approximations in one and two dimensions. We start this section with some useful definitions and merits of using the Riemann-Liouville fractional derivatives [16,31,56,57,63].
Definition 2.1. The left and right R-L fractional derivatives of order α, (n − 1 < α < n) of the function u(x) on a closed interval [a, b] are given as [56] 1. left R-L fractional derivative: , Ω ∈ R then, the Fourier transforms of the left and right R-L fractional derivatives are found satisfying The shifted Grünwald operator [34] approximates the R-L fractional derivative uniformly with first-order accuracy, that is, where r is a positive integer and g (α) k = (−1) k ( α k ). Obviously, the coefficients g (α) k in equation (3) are the same as the coefficients of the power series function (1 − z) α , where for all |z| ≤ 1, which can be recursively evaluated using The coefficients in (3) are also found satisfying the properties for 1 ≤ α ≤ 2 [35,56].
Definition 2.2. Assume that f is a differentiable function, Let α be a real number bounded in the interval 0 ≤ α ≤ 2, then the derivative with order α with power law is defined as [9,56] RL

KOLADE M. OWOLABI
Motivated by the shifted operator in (3) and multistep scheme, we obtain the second-order approximation for the Riemann-Liouville fractional derivatives as follows.
x u, and let its Fourier transform be L 1 (R), the weighted and shifted Grünwald difference operator is defined by where r and s are symmetric integers, that is, L D α ,r,s u(x) = L D α ,s,r u(x), for r = s.
Proof. By using the definition of B α ,r in (3), conveniently, we can rewrite the weighted and shifted Grünwald difference operator as which by the Fourier transform results to where (7) and (8) that there exists By considering a function u(x) bounded interval [a, b], assuming u(a) = 0 or u(b) = 0, then the function u(x) can be zero extended for x < a or x > b. The L-R Riemann-Liouville fractional order derivatives of u(x) at point x can be approximated by the weighted and shifted Grünwald difference operators with second order accuracy where With the choice (r, s) = (1, 0), (1, −1), the discretized approximations (10) for the L-R Riemann-Liouville fractional derivatives are simplified into where (r, s) = (1, 0), ω The main advantage of using the Riemann-Liouville fractional derivative within the variation principles is due to its flexibility and the possibility of defining the integration by parts as the fractional reaction-diffusion problems become the classical case whenever α becomes an integer [50]. An additional merit is that it satisfies virtually all the mathematical principle under the scope of fractional calculus, especially when using Laplace transform to obtain initial condition with fractional power index (exponent) which is actually realistic in practical and mathematical point of views. Hence, most researchers [10,50,62] applauded the use of the Riemann-Liouville definition due to its suitability and adaptability. For details, we refer our readers to some classical books and research papers [16,31,32,50,56,[67][68][69][70].
Here, we consider a one-dimensional two-sided fractional-in-space diffusion equation of the form where α ∈ (1, 2], u(x, t) is the species concentration or population density in position x and time t. The terms a D α x and x D α b are the space fractional Riemann-Liouville operators of order α, for 1 < α ≤ 2. The parameters γ > 0 and β > 0 are the diffusion tensors or conductivities with γ 2 + β 2 = 0, and f represents a local source term. In what follows, we discretize the space fractional diffusion equation (12) by using the weighted and shifted Grünwald difference formula (11).
We divide the interval of integration [a, b] into subinterval of uniform mesh with step size = (b − a)/N and let the time step κ = T /M , for N > 0, M > 0 are both integers. The grid points are represented by x i = i and t n = nκ for 1 ≤ i ≤ N and 0 ≤ n ≤ M . We let t n+1/2 = (t n + t n+1 )/2 for 0 ≤ n ≤ M − 1, and in one dimension we adopt the notations i − u n i )/κ. Next, by using the technique of Crank-Nicolson scheme to discretize (12) in time, we have For the discretization in space, we use the weighted and shifted Grünwald difference operators L D α ,r,s u(x, t) and R D α ,r,s u(x, t) for the respective approximation of left and right Riemann-Liouvillve fractional derivatives a D α x u(x, t) and where | n i | ≤ê(κ 2 + h 2 ). Multiplying equation (13) by the time-step κ and simplifying the time layers, leads to Substitute for L D α ,r,s u and R D α ,r,s u in (11), we get By denoting the numerical approximation of u n i as U n i , we derive the Crank-Nicolson scheme using the weighted and shifted Grünwald difference scheme for problem (12) as For the sake of convenience and easy implementation, we write to describe the finite difference scheme (14) as with diagonals ω (α) k , k = 0, . . . , n − 1 are the coefficients given in (11) which correspond to (r, s) = (1, 0) or (r, r) = (1, −1), and (17) To justify the suitability and accuracy of our method in one dimension, as shown in Table 1, we compute the norm infinity and norm relative errors of the following fractional reaction-diffusion problem subject to initial condition u(x, 0) = x 1+α , x ∈ (0, 1], and the boundary conditions The exact solution is given as u(x, t) = e −t x 1+α . Obviously, the numerical experiment in Table 1 shows that as N is increasing, so also the accuracy of the scheme is improving for all values of α. Finally, Figure 1 shows the convergence result for one-dimensional problem (18) at some instances of α and final simulation time t = 2 and N = 200.
Following the 1D example discussed above, we establish the Crank-Nicolson method by using the weighted and shifted Grünwald difference formula (11) for problem (19). We likewise divide the computational domain Ω into subinterval of equal mesh with space steps x = (b − a)N x , y = (d − c)/N y and time-step κ = T /M , where N x , N y , M are positive integers. The grid points are given by x i = i x , y j = j y and t n = nκ for 1 ≤ i ≤ N x , 1 ≤ j ≤ N y and 0 ≤ n ≤ M . We also let t n+1/2 = (t n + t n+1 )/2 for 0 ≤ n ≤ M − 1, so as to adopt the notations For the space discretization, we use the weighted and shifted Grünwald difference operators L D α x ,r,s u, R D α x,r,s u and L D η y ,r,s u, R D η y ,r,s u to approximate the fractional terms a D α x u, x D α b u and c D η y u, y D η d u respectively. On multiplying (20) by the time-step κ, we have For the sake of simplicity, we choose equal step sizes as = x = y . By adopting the Taylor expansion, we get By adding (22) to (21) and simplifying further yields If we denote the numerical approximation to u n i,j by U n i,j , the finite difference approximation to 2D fractional-in-space diffusion problem (19) is given as Note that a simple calculation indicates From equations (23) and (25), we have

KOLADE M. OWOLABI
By getting rid of the truncation error and denoting the numerical approximation of u n i,j by U n i,j , we obtain To solve (24) efficiently in higher dimensions, the following second-order schemes can be adapted. The linear multistep implicit-explicit (IMEX) schemes [42] (2θ The choices of parameters (θ, ϕ) could results so various second-order schemes. For example, if we let (θ, ϕ) = ( 1 2 , 0) we obtain the most popular IMEX scheme called CNAB [2], which combines both second-order Adams-Bashforth for the nonlinear part and Crank-Nicolson method for the linear term. Details on derivation and stability properties can be found in [42]. The second-order exponential time differencing Runge-Kutta scheme of Cox and Matthews [21], which is denoted for brevity as the ETDRK2: a n = u n e Lh + N(e Lh − 1)/c, where N and L are the operators that represent the reaction and diffusion terms respectively.
In Table 2, we justify the accuracy and efficiency of the method in two dimensions by reporting the norm infinity and nor relative errors at some instances of fractional power index α at t = 1.0. We also calculate the timing results for the fractional diffusion problem u t = 0 D α x u(x, y, t) + x D α 1 u(x, y, t) + 0 D η y u(x, y, t) + y D η 1 u(x, y, t) + f (x, y, t), t > 0 (29) with the source term in the domain Ω = (0, 1) 2 , subject to initial condition u(x, y, 0) = x 3 (1 − x) 3 y 3 (1 − y) 3 and boundary conditions u(x, y, t)| ∂Ω = 0. The exact solution of (29) is u(x, y, t) = e −t x 3 (1 − x) 3 y 3 (1 − y) 3 . Convergence result at different instances of α for two-dimensional problem (29) is shown in Figure 2 with t = 2 and N = 200. Table 2. The norm infinity and norm relative of errors for one dimensional problem (29 ) obtained at some instances of fractional power α at final time t = 1.0 approximated with the Crank-Nicolson weighted and shifted Grünwald difference scheme for η = 1.8, κ = .   3. Multicomponent fractional reaction-diffusion system. The key role of diffusion in the modelling of several special phenomena that are largely encountered in engineering and science has been studied extensively. Diffusion and cross-diffusion in both classical [48,49] and fractional [44,46,55] scenarios, as the causes of the spontaneous emergence of ordered structures, known as patterns [40,41].
In this section, we are going study two fractional reaction-diffusion systems that modelled spatial interactions between a prey and the predators. We shall examine the local analysis at the constant positive steady state, and investigate the global existence and bifurcation of stationary patterns of the system. Assuming in this paper that predator-prey species spatial interactions is experimented in a confined and bounded domain Ω ∈ R N , thus we are led to consider the general fractional-in-space reaction-diffusion system where u i = (u 1 , . . . , u n ) is the vector of concentration or species density at time t and position x, the term F i (u i ), for i = 1, 2, . . . , n, accounts for the local kinetics of the system. The above equation is considered subject to initial condition u(x, 0) = u 0 (x), for i = 1, 2, . . . , n with either the Neumann or Dirichlet boundary conditions. For the analysis of n−species system, see [36,60,64]. For a two-component example, we consider the Gray-Scott system where ∆ α = ∂ α /∂x α is one-dimensional space fractional Laplacian operator, u = u(x, t) and v = v(x, t) are the concentrations of the two chemical species for the inhibitor u and activator v. The parameter A stands the rate at which the inhibitor u is fed from the reservoir into the reactor, while b represents the overall rate of decay of the activator due to draining. The terms δ 1 , δ 2 are the diffusion parameters. Derivation and details analysis of (31) can be found in [36]. For a three-component system, we consider where ∆ denotes the Riemann-Liouville fractional Laplacian operator of order α, in the super-diffusive scenarios 1 < α ≤ 2. The case α = 2 corresponds to the standard reaction-diffusion system. The parameter c is the increasing death rate of species w, a ∈ (0, 1) is the conversion coefficient that w is perfectly recovered and becomes v. The positive constants δ 1 , δ 2 and δ 3 are fractional diffusion coefficients, and the initial data u, v, w(x, 0) are continuous functions. Here, n represents the outward unit normal vector of the smooth boundary and the homogeneous Neumann boundary condition is considered. It is not difficult to see that the spatially homogeneous form of (32) obtained by setting δ 1 = δ 2 = δ 3 = 0, has positive steady state if and only if If (33) is satisfied, we can uniquely give the positive steady state bŷ then the steady state solutionẑ = (û,v,ŵ) T for the spatially homogeneous form of(32) is local asymptotically stable.
Proof. See Wang [65] for a similar proof to Theorem 3.1.     (32) obtained by setting δ 1 = δ 2 = δ 3 = 0 when a steady state structure (type II, Turing pattern ) is established in the domain. At initial time, the three species oscillate in phase, but as simulation time is increased, the dynamics enter into its steady state. Panel (d) is the species attractor, it further suggest a direction to gain useful insight to pattern formation.
3.1. Local analysis at steady state. In order to give a good working guidelines on the appropriate choice of parameters for the numerical simulation of full fractional reaction-diffusion system, it is mandatory to consider the local dynamics of the system. It is, naturally, the dynamics in the biologically meaningful cases u ≥ 0, v ≥ 0 and w ≥ 0 that are of great interest [43].
In the spirit of [52], we introduce a cross-diffusion term in the third equation in order to search for the diffusion-driven instability. We give the corresponding steady state of problem (32) by where k and σ are positive parameters, k is known as the cross-diffusion constant.
In this system, species w diffuses with the flux We realized in this study that, as 2kvw(σ + v −2 ) ≥ 0, the part {2kvw(σ + v −2 )}∇v of the flux results to increase in the population of species v. Readers are referred to [37,39] for details.

Global analysis and bifurcation of stationary patterns formation.
In this section, we discuss the global existence and bifurcation of stationary patterns to (34) with respect to diffusion constant δ 1 , as we fixed other parameters a, b, c, σ, δ 2 , δ 3 and k.
Proof. Since (34) has at least one non-constant positive solution, we shall establish the proof for any δ 1 ≥ P by contradiction. Consider P > 1/µ 2 , let δ * 3 > 0 satisfying and we consider the problem ∂ n z = 0, on ∂Ω.

KOLADE M. OWOLABI
We can say that z is a positive solution of (34) if and only if it is a solution of (37) for t = 1.Obviously, z is the unique solution of (37) for any 0 ≤ t ≤ 1. As previously observed, z is positive solution of (37) if and only if , and by direct calculation, we have It follows from (38), in view of conditions It is noticeable enough to see that zero is not an eigenvalue of the matrix dim E(µ i ) = d n , which is odd.
By assumption, we can see that both equations G(1; z) = 0 and G(0; z) = 0 possess only the positive solution z ∈ B(C). Thus, by (39) and (40) which is contradiction of (41). Hence, the proof is completed.
Proof. The proof is based on the topological degree argument as mentioned earlier, readers are referred to [53] for details to similar treatment.

Numerical experiments.
In this section, we shall focus primarily on the implementation of the numerical methods discussed in Section 3, when applied to solve full time-dependent fractional-in-space reaction-diffusion systems (31) and (32) for the superdiffusive process in the region 1 < α ≤ 2. In the simulation framework, numerical results in one dimensions are easily undertaken by using either method of lines in conjunction with spatial adaptive (implicit-explicit [42], exponential time differencing [47]) schemes, or even just simple Crank-Nicholson and finite difference scheme or finite element collocation methods [15]. It is in higher dimensions that the ideas presented in this paper really become of serious value. 4.1. Two-dimensional results. In two dimensions (2D), we experiment both (31) and (32) numerically on a fixed square domain [0, L] × [0, L], L > 0, subject to zero flux boundary condition and the random (smooth) initial conditions, and u, v, w ∈ R. It should be noted that L must be chosen large enough to give room for the waves to propagate. We let x, y ∈ R and ∆ α u = ∂ α u ∂x α + ∂ α u ∂y α denotes the fractional Laplacian operator of order α in two-dimensional space, with similar expressions apply for v(x, y, t) and w(x, y, t).
In Figure 4, we simulate with a = 0.035, b = 0.065, δ 1 = 2.05e−5 and δ 2 = 1.0e−5 using 128 Fourier nodes on domain size L = 2. Simulation runs for final time t = 5000. We also observed in the experiment that the spatial evolution of species are similar, as a result, we report only the distribution for v−species. It should be mentioned that apart from the spiral patterns reported here, other Turing patterns such as the mitotic-spot, pure spots as well as many other chaotic-spatiotemporal structures are obtainable, depending on the choice of the key parameters a, b, the initial function and variability of other parameters. See Pearson [36,54] for details.
Similarly for the three-component system (32) two-dimensional results showing spatiotemporal oscillations of the species at some instances of fractional order α is displayed in Figure 5. The first and second columns correspond to α = 1.35 and α = 1.95 respectively.

4.2.
Three-dimensional results. The dynamic richness of space fractional reaction-diffusion equations is further explored by considering the 3D experiments. Here, we let ∇ α = (∂ α /∂x α + ∂ α /∂y α + ∂ α /∂z α ), for u = (x, y, z, t), v = (x, y, z, t) with similar expression for w, subject to the smooth initial condition for the twocomponent system. The 3D chaotic spiral evolution of (31) at α = 1.15 (top-row), α = 1.55 (middle-row) and α = 1.75 (bottom-row) for both species is shown in Figure 6. The dots or spots patterns as reflected in the first column representing the distribution of u−species show the emergence of the irregular spot or mitotic-spot patterns which disintegrate from the boundary to the center of the domain to form the Turing spiral patterns, as displayed in the second column. Clearly, there is an agreement between the 2D and 3D results.
For the three-species case, we consider the initial functions [44,50] u(x, y, z, 0) = 1 − exp −10((x − p 2 ) 2 + (y − p 2 ) 2 + (z − p 2 ) 2 ) , v(x, y, z, 0) = exp −10((x − p 2 ) 2 + 2(y − p 2 w(x, y, z, 0) = 0.5 − exp −10((x − p 2 ) 2 + (y − p 2 ) 2 + 6(z − p 2 ) 2 ) , chosen to specifically induce a nontrivial spatiotemporal dynamics on a fixed but large domain size ±L, L = 5 to enable the waves to have enough room to propagate. In Figure 7, evolution of species u, u and w obtained at some different instances of α are shown in columns 1-3 respectively. In the 3D experiments for system (32), we observed that the distributions of the species are not very similar unlike the classical case (α = 2). Hence, we display for different instances of α in the regime 1 < α ≤ 2 as shown in rows 1-3, for α = 1.35, 1.50, 1.95 respectively, and the fourth-row corresponds to the classical case. We also noticed that the 3D pattern formation in the classical case is almost similar to the fractional scenarios in the supperdiffusive regime.

5.
Conclusion. The research work here proposed a better alternative second-order fractional derivative approximation to the commonly used finite difference scheme, via the weighted and shifted Gr unwald difference approximation in conjunction with the Crank-Nicholson scheme when applied to solve the general fractional-inspace partial differential equations. Effectiveness and accuracy of the method are proved theoretically and confirmed numerically by reporting the norm infinity and norm relative errors in one and two dimensions as well as their convergence results. We provide an extension to 3D by solving full fractional-in-space reaction-diffusion system. Numerical simulation results of the two fractional reaction-diffusion systems reveal an amazing variety of irregular spatiotemporal patterns. Further, we have been able to give enough proof by the computer simulations that pattern formation in the fractional scenario is practically the same as in the classical case, our assertion is evident in Figure 7. The availability of efficient numerical schemes for the integration of fractional-in-space reaction-diffusion systems is important, as these equations are becoming the key tools for describing a wide range of systems. Though we have discussed a numerical approach for the solution of space fractional super-diffusive order (1 < α < 2) equations in this work, readers should be noted that the procedures presented could be extended to time-space fractional integrodifferential equations. Furthermore, the discretization approach reported in this paper can be extended to solve any higher dimensional practical problems.