Fourier Spectral Method for Solving Fractional-order System

Abstract In this paper, we have studied a new fractional reaction-diffusion two-species system as an extension to the Rosenzweig-MacArthur reaction-diffusion di-trophic food chain system which models the spatial interactions between a prey and predator. To guarantee good working guidelines when numerically simulating the model, we first show that the system is locally asymptotically stable, as it provides good conditions and correct choice of ecological parameters to enhance a biologically meaningful result. We propose a fast and accurate method for numerical solutions of space fractional reactiondiffusion equations. The technique is based on Fourier spectral method in space and exponential integrator scheme in time. The complexity of fractional derivative index in fractional reaction diffusion model is numerically formulated and graphically displayed in one-, twoand three-dimensions.


Introduction
Historically, the earlier mathematical approach for two-variable system was reported to model spatial interactions of predator-prey dynamics of the Lotka-Volterra type. It had since become a good testing ground for describing population dynamics. Years later, more realistic and interesting two-species models were proposed. For instance, May [16] proposed in a model that the prey population grows logistically. The linear stability analysis result of (1) [5] indicates that the coexistence equilibrium, dependent of the carrying capacity κ of the prey population density is always a stable point. Owing to the assumption that the predator-prey model underlies a saturation, we have a model attributed to Rosenzweig-MacArthur di-trophic food chain system [6,29] du dt = τ u 1 − where u and v denote the population of prey and predator respectively for t > 0, measure as a number of individuals, biomass or density, τ represents the prey biotic potential or intrinsic growth rate, κ is the environmental carrying capacity of the prey, 76 Fourier Spectral Method for Solving Fractional-order System σ 2 defines the predator maximum predation rate, γ 2 is the predator efficiency and ρ, is known to be half saturation amount (the rate of prey required to get one-half of σ 2 ). The death rate of predator is denoted by δ 2 .
In this paper, we begin our study with space fractional diffusion equation where U(℘, t) is the species population density distribution at positions ℘ = ℘(x, y, z) and time t, it could also be regarded as the concentration of chemical species. The space fractional operator δ ∂ β ∂℘ β (℘, t), 0 < β < 2 is given on the basis of the Fourier transformation, with variable δ > 0 as the diffusion coefficient. Here, we define the Fourier transformation and its corresponding inverse asF Next, we apply the space fractional operator on the species population density U( Prior to definition in (5), Fourier transform of (5) results to an ordinary differential equation (ODE) of the form subject to initial function The solution of (6) is sought asÛ which results toÛ with corresponding inverse Fourier transform Readers are referred to classical books [13,27,28,42] for details of the preliminary theories and definitions of fractional derivatives. The aims of the present paper are in folds; In Section 2, we give an extension to the commonly two-variable Rosenzweig-Macarthur reaction-diffusion system and examine the linear stability analysis of the new model. In Section 3, we formulate a good versatile numerical techniques based on Fourier spectral method in space and exponential time differencing method for simulating the fractional reaction-diffusion problem of Rosenzweig-MacArthur type with delay in one and two dimensions. Numerical experiments and concluding remarks are given in Sections 4 and 5 respectively. In this work, we give an extension to model (2) in the form of fractional reaction-diffusion system subject to the initial condition u(℘, 0) = u 0 (℘), v(℘, 0) = v 0 (℘), ℘ ∈ Ω and the zero-flux boundary condition where ν is the outward normal to to the boundary ∂Ω, u(℘, t) ∈ R and v(℘, t) ∈ R are the respective population densities of prey u and predators v at time t and position x. Parameters σ, τ, γ and δ are strictly positive, D u > 0 and D v > 0 are the prey-predator diffusion coefficients. The parameter β > 0 is the fractional power that distinguishes between the classical and fractional reaction-diffusion systems. If β = 1, equation (7) reduces to the classical reaction-diffusion system. If 0 < β < 1, we have a sub-diffusive case and if β > 1, a situation that corresponds to super-diffusive system is attained.
In order to give a good working versatile guidelines when numerically simulating the full fractional reaction-diffusion system, it is important to provide details of the local dynamics of such system. local analysis provides appropriate conditions and good choice of parameters for the solutions to have a biologically meaningful equilibria. We shall focus on the dynamics of interest in the region u ≥ 0, v ≥ 0 that is biologically meaningful. The nullclines are the solution curves that correspond to the equations f (u, v) = 0, g(u, v) = 0.
Linear stability analysis reveals that system (7) is trivially satisfied if τ = 0, a point that corresponds to total extinction of both species, so the trivial equilibrium point (0, 0) is a saddle. We are not interested in the trivial result because it is of no interest biologically. Hence, we are primarily concerned with an equilibrium state which is either stable or unstable corresponding to the coexistence of both species. The point is given by The equilibrium point E s is locally asymptotically stable if κ < (σ + δ)/(σ − δ), the point becomes unstable when a Hopf bifurcation occurs at H b = (σ + δ)/(σ − δ). Consider the case τ > 0, for σ > δ exp(γτ ) and κ < δ/(σ exp(−γτ ) − δ), a point corresponding to the existence of only prey, that is E 1 = (κ, 0), which is stable for all values of τ ≥ 0. So also, for σ > δ exp(γτ ) and κ > δ/(σ exp(−γτ ) − δ), there is another semi-trivial state E 2 = (0, κ) which corresponds to the predators existence. There is also one nontrivial steady state where the two species coexist, and whose stability depends on parameter τ . We find the parameter set in which stability can occur as: We know that the stability of point E s depends on the value of the delay τ ≥ 0. Stability properties of E s are found by the real part of the characteristic roots, such that and By simplifying (8), we obtain .

78
Fourier Spectral Method for Solving Fractional-order System It should be noted that both τ and κ can be chosen as bifurcation parameters for system (7). for instance, for τ = 0 and κ ∈ (0, H a ) the point E 2 is stable, and the coexistence point E s is locally asymptotically stable for κ ∈ (H a , H b ). Whenever κ H a , a limit cycle occurs. It is noticeable that the destabilization of the equilibria is shift forward due to inclusion of a delay. That is, if τ > 0 and κ ∈ (0, Eû(τ )), with Eû(τ ) ≤ Eû(0) = H a , then the point E 2 is stable, but if κ ∈ (Eû(τ ), Eŝ(τ )), with Eŝ(τ ) ≤ Eŝ(0) = H b , then the coexistence equilibrium state E s is locally asymptotically stable. it should also be noted that an increase in the delay may give rise to damped oscillations whenever κ > Eŝ(τ ), as depicted in Figure 1.

Numerical technique for the fractional reaction-diffusion system
The fractional reaction-diffusion Rosenzweig-MacArthur system (7) supports the use of two classical methods, since the problem can be divided into the linear and nonlinear parts [33,37]. As a result, we formulate a Fourier spectral method for the spatial discretization and employ higher order exponential time differencing scheme to advance the resulting system of ordinary differential equations in time.

Formulation of adaptive methods in space and time
In recent years, credit has been given to spectral methods over the conventional finite differences [7,8,20], due to its capability of removing the stiffness issue associated with the of a reaction-diffusion problems. From the known theory of integrating factor technique, we shall formulate the spectral method in two spatial dimensions. In the spirits of [20,25], we apply the integrating factor technique to the Fourier transform systems (7), to obtain where U and V are the double Fourier transforms for the species densities u(x, y, t), v(x, y, t). In other words, To explicitly circumvent the issue of stiffness in the second partial derivatives, we let Ω β = (χ β x + χ β y ), and set Next, we need to discretize the square domain by considering the equispaced points N x and N y in the spatial directions of x and y. We use the discrete fast Fourier transform (DFFT) [33] to transform (9) to a system of ordinary differential equations (ODEs) where Boundary conditions are now set at extremes of the domain size [−L, L]. Now, the system has been converted to ODEs, the stiffness issue has equally been removed. It should be noted that any explicit higher-order time stepping integrators can be used.

Numerical experiments
In this section, we want to consider the nontrivial example of nonlinear Rosenzweig-MacArthur fractional partial differential equations for pattern formation in one and two dimensions. Numerical simulations of the di-trophic food-chain models are provided to demonstrate and compare the asymptotic nature of the time-dependent fractional reaction-diffusion system.

One-dimensional experiment
We first illustrate the asymptotic stability of equilibrium E s result which correspond to the existence of the two species through numerical simulations of model (7) in absence of diffusion, that is D u = D v = 0. In Figure 1, effect of the delay on the solution of the system (7) is shown at different instants of delay τ = 0.5, 1, 1.5 for columns 1,2,3 respectively. The first and second rows indicate the time series solutions and their corresponding limit cycles at time t = 100. Clearly, increase in the delay τ has given rise to damped oscillations.
In Figure 2, we illustrate also the stability of the equilibrium E s via numerical simulations of the fractional reactiondiffusion system (7) at instants value of delay τ , with ecological parameters: {κ = 0.8, σ = 4, γ = 0.4, δ = 0.6, D u = 0.2, D v = 0.05}. These parameters and initial populations u(0) = 0.5, v(0) = 0.5 are found satisfying the conditions for asymptotic stability in Section 2. It is noticeable that as τ is increased from 0.2 to 1.5, both species oscillate in phase. So, the numerical simulation results in both Figures 1 and 2 agree with the theoretical result in Section 3.   To demonstrate the reliability and efficiency of our proposed method, we report the convergence results in space as well as the computational time for the fractional Rosenzweig-MacArthur system (7) for varying values of the fractional power β. Table 1 shows the convergence results of the fractional Fourier transform with different instants values of α at final time t = 1. Reference solutions were calculated by evaluating (7) with 2 10 Fourier modes. Obviously, the fractional Fourier approach is able to gain spectral convergence up to machine precision regardless of chosen β values.
In the simulation experiments, various type of dynamics are observed and we realiazed that both species have similar distributions. As a result, we streamlined our analysis of pattern formation to only one distribution (we report that of species v in this paper). In Figure 3, we perform numerical simulations in 2D for different instants of fractional power β = 1, 1.  Figure 4 at β = 2 which we clasify as sub-diffusive for 0 < β < 1, diffusive for β = 1 and super-diffusive in the region 1 < β < 2.   Figure 4. Contour plots of (7) in 2D showing effect of fractional power at different instants of β = 0.5, 1, 1.5 which correspond to sub-diffusive, diffusive and super-diffusive scenarios at t = 1000.
In Figure 6, we increase the simulation time from t = 5 to t = 50 to obtain various object patterns such as diamond, star, oval or egg-like and dice shapes. It should be mentioned that, with increasing time and perturbation of the initial data, other Turing dynamics can be observed.

Conclusion
In this work, we have demonstrate the use of Fourier spectral method in conjunction with the exponential time differencing scheme to numerically simulate the Rosenzweig-MacArthur fractional reaction-diffusion system with delay in the presence of fractional derivative index β. Our mathematical analysis results tallies with the numerical results as shown in Figures 1 and  2, which justify the case for when the delay τ is increasing, both species oscillate in phase which can result to homogeneous and non-homogeneous oscillations. The 2D computer simulations of the fractional reaction-diffusion system (in the regions β > 1) has given a strong proof that pattern formations in fractional case is practically the same as in the standard reactiondiffusion scenario at β = 1. With the variation in the choice of parameters, one can observe the rich variety of spatiotemporal behaviour. The mathematical ideas presented in this work can be extended to solve both the integer and non-integer order fractional reaction-advection problems in one and high dimensions.