{'fraction_non_alphanumeric': 0.05114353511384864, 'fraction_numerical': 0.021759983844095523, 'mean_word_length': 4.317583892617449, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An efficient spectral-Galerkin method for fractional reaction-diffusion equations in unbounded domains

In this work, we apply a fast and accurate numerical method for solving fractional reaction-diffusion equations in unbounded domains. By using the Fourier-like spectral approach in space, this method can effectively handle the fractional Laplace operator, leading to a fully diagonal representation of the fractional Laplacian. To fully discretize the underlying nonlinear reaction-diffusion systems, we propose to use an accurate time marching scheme based on ETDRK4. Numerical examples are presented to illustrate the effectiveness of the proposed method.


Introduction
Progress in the last few decades has shown that many complex systems in science, economics and engineering are found to be more accurately described by fractional differential equations, see [6,8,28,45]. Usually, in these equations fractional derivatives are involved to deal with long range time dependence or long range spatial interactions. Among them, typical fractional operators are Reimann-Liouville fractional derivatives, Caputo fractional derivatives, Riesz fractional derivatives, fractional Laplacian and so on. Since analytic solutions to fractional differential equations are usually unknown, there has been a growing interest in developing effective and efficient numerical methods, see [2,5,10,13,25,44,46] and the references therein.
In this paper, we are going to deal with the fractional reaction-diffusion equations in unbounded domains where the second order space derivative is replaced with fractional Laplacian. To better illustrate our idea, we consider the following model equation for s ∈ (0, 1.5), where D > 0 is the diffusion coefficient and f (u) is the reaction term.
Reaction-diffusion equations model the change in space and time of the concentration of substances under the influence reaction and diffusion, and have found wide applications in biology, physics, chemistry and engineering, see, for example, [7,21,36]. Examples of particular interest include the Fisher-Kolmogorov equation when f (u) = ru(1 − u K ), see [34,47], and the Allen-Cahn equation when f (u) = u − u 3 , see [19,24,42]. For numerical methods, readers may see [4,31,35] and the references therein. For integer order reaction-diffusion equation, the use of integer order derivatives lies in the assumption that the random motion is stochastic Gaussian process. However, more and more studies has shown that certain systems exhibit anomalous diffusion, see, e.g., [12,37], and the references therein. In this paper, we study the case when the second-order space derivative is replaced with the fractional Laplace operator, which can account for anomalous diffusion phenomenon by allowing arbitrary jump length distribution, including the case with infinite variance. For time-fractional reaction-diffusion models, readers may refer to [3,16,18,37].
Fractional Laplacian can be defined as a pseudo-differential operator via the Fourier transform for s > 0: (

1.2)
Or equivalently, it can be defined for s ∈ (0, 1) by singular representation (cf. [30,Prop. 3.3]): where "p.v." stands for the principle value and the normalisation constant is given as Various numerical treatment of the fractional Laplace operator has been proposed, see, e.g., finite difference method [14,15,29,44], finite element method [1,5,9], and spectral method [27,41,43]. Recently, in [40], biorthogonal mapped Chebyshev functions are proposed as the basis for fractional differential equations in unbounded domain, which then lead to diagonal stiffness matrix by using the Dunford-Taylor formula. This provides an efficient tool for numerical approximations of fractional reaction-diffusion equations, especially in multidimensional problems. In this paper, we will follow the routines in [40] to deal with the fractional Laplace operator. The aim of this work is to propose a reliable and efficient approach to solve the fractional reaction-diffusion equation in unbounded domains involving the fractional Laplace operator. In order to obtain higher accuracy, we adapt the fourth-order exponential time differencing Runge-Kutta (ETDRK4) method for time stepping, which is based on the exact integration of the linear part and the approximation of an integral involving the nonlinear part.
The main contributions of this work are three folds: (1) Since fractional derivatives are usually nonlocal, here we study the unbounded domain directly, different from other previous works where domain truncation is usually used. (2) The biorthogonal mapped Chebyshev functions are used as the basis for our spectral-Galerkin method, which lead to diagonal stiffness matrix, and thus very efficient even in higher dimensions. (3) The ETDRK4 time stepping method is used for the resulted semidiscrete problem which is 4th-order accurate in time. This ensures that numerical treatment of the fractional reaction-diffusion equation is both more efficient and accurate.
The rest of this paper is organized as follows. In the next section, we briefly give some preliminaries on the properties of biorthogonal mapped Chebyshev functions and the Dunford-Taylor formula for s ∈ (0, 2). The numerical approach based on spectral-Galerkin method and ETDRK4 time stepping method for fractional reaction-diffusion equation is developed in section 3, where the linear stability of ETDRK4 is also analysed. In section 4, various numerical examples are provided to illustrate the effectiveness of proposed method, including the Fisher-Kolmogorov equation, Allen-Cahn equation, Gray-Scott equation and FitzHugh-Nagumo equation. The final section is for some concluding remarks.

Preliminary
The mapped Chebyshev functions (MCFs) based on the Fourier-like biorthogonal property was proposed and studied recently in [40]. The biorthogonal property enables us to make the fractional Laplacian fully diagonalised, and the complexity of solving an elliptic fractional PDE is quasi-optimal, i.e., O((N log 2N ) d ) with N being the number of modes in each spatial direction. The present section is to provide some basic preliminaries useful for the efficient spectral algorithms to be provided in the forthcoming section.
Definition 2.1. Introduce the one-to-one algebraic mapping and define the MCFs as T n (x) = 1 c n π/2 1 − y 2 T n (y) = 1 c n π/2

2)
for x ∈ R and integer n ≥ 0.
We have the following important property of the MCFs with respect to the mass matrix and stiffness matrix, see [39,Proposition 2.4].
Proposition 2.1. The MCFs are orthonormal in L 2 (R), and we have for the stiffness matrix 3) The Fourier-like biorthogonal MCFs { T n } N p=0 are defined as a linear combination of the mapped Chebyshev functions, where the expansion coefficients are the matrix formed by the orthonormal eigenvectors of S. Readers may refer to [40] for more details. Then they satisfy the following property: form an equivalent basis with {T n } N p=0 , denoted as V N , and they are biorthogonal in the sense that
which is the tensor product of d copies of V N . Define the d-dimensional tensorial Fourierlike basis and denote the vector of the corresponding eigenvalues in (2.4) by where the index set Υ N := n = (n 1 , · · · , n d ) : 0 ≤ n j ≤ N, 1 ≤ j ≤ d .
2.3. Dunford-Taylor formula. Besides the biorthogonal MCFs, the other key component of the spectral-Galerkin method in [40] is the following Dunford-Taylor formula for fractional Laplacian.
where I is the identity operator and Proof of the Dunford-Taylor formula can be found in [5]. In fact, similar results can be obtained for s ∈ (1, 2) by dividing s into s 1 + s 2 , where s 1 , s 2 ∈ (0, 1). (2.13) Proof. By Parseval's identity, we have (2.14) Note that Then we have

Spectral-Galerkin method
In this section, we conduct the numerical approximations for (1.1), where we first formulate the spectral-Galerkin method to deal with the space fractional derivatives, and later advance the resulting ODE in time with ETDRK4.
Then using the Dunford-Taylor formula in Lmm. 2.3 for s ∈ (0, 1), our spectral-Galerkin method is to find For s ∈ (1, 1.5), using the Dunford-Taylor formula in Thm 2.1, the spectral-Galerkin method is to find and v N = T n (x), n = (n 1 , n 2 , · · · , n d ) ∈ Υ N , then similar to [40, Thm 3.2], we have the following unified semi-discrete form for all s ∈ (0, 1.5) We remark that a diagonal representation of the fractional Laplace operator is obtained based on the biorthogonal MCFs in (2.6) and the Dunford-Taylor formula in Lmm. 2.3 and Thm. 2.1 yields a unified semi-discrete form for all s ∈ (0, 1.5).

Time discretisation.
We first write the resulted ODE system (3.6) in the following more general form: where L and N denote linear and nonlinear terms, respectively. In our case, L is diagonal matrix with the components −D|λ n | s 1 , and N(U ) represents the expansion coefficients of reaction term.
Typical numerical methods for solving problems of this kind include implicit-explicit method, integrating factor method and exponential time differencing (ETD) method, etc. To remove the stiffness caused by the linear part, we choose the ETD method which involves exact integration of the linear part followed by an approximation of the integral of the nonlinear term. Multiplying (3.8) by the term e −Lt and integrating over a single time step τ give For the approximation of the integral involving the nonlinear term, Runge-Kutta method is taken since it is easy to implement and requires only one previous evaluation when compared to multistep methods. ETDRK4 method was first proposed in [11], and then modified in [20] to increase numerical stability. In this work, we utilize the fourth-order Runge-Kutta scheme of Krogstad in [22], denoted as ETDRK4-B, which has smaller local truncation error and larger stability properties. More precisely, for the semi-discrete ODE system (3.8), the ETDRK4 scheme is where τ is the time stepsize and the stages µ 2 , µ 3 , µ 4 are defined as with the functions ϕ 1 , ϕ 2 , ϕ 3 given as 3.3. Linear stability analysis. The stability analysis of the ETDRK4 method is as follows (see also [22,32,33]). For the nonlinear ODE we suppose that there exists a fixed point U 0 such that λ U 0 + N(U 0 ) = 0. Linearizing about this fixed point leads to Here U is a perturbation to U 0 , and ρ = N ′ (U ) at U = U 0 . Then for this ODE, the fixed point U 0 is stable if Re(λ + ρ) < 0. The application of ETDRK4 method to (3.11) leads to a recurrence relation involving U n and U n+1 . Introducing the notation x = ρτ , y = λτ , we obtain the following amplification factor where c 0 = e y , We observe that as y → 0, then the fourth-order Runge-Kutta method is recovered where In Fig. 3.3, we give the stability region of ETDRK4 method for several different nonpositive y, we observe that stability region grows as the absolute value of y increases. We also present in the right figure of Fig. 3.3 the maximum and minimum eigenvalues of the stiffness matrix S for different N . It is observed that maximum eigenvalues increase with N approximately as N 2 while minimum eigenvalues decrease as N −2 . Then for the resulted ODE system (3.8) from our spectral-Galerkin method, the stiffness ratio is roughly N 4s . Thus proper treatment of the linear part is necessary to reduce the excessive restriction on step size and this is the advantage of the ETD method. At the same time, minimum eigenvalues go to 0 as N increases, then stability region will shrink to the standard fourth-order Runge-Kutta method as N goes to infinity.

Numerical examples
Several numerical examples will be carried out in this section to demonstrate the effectiveness of the spectral-Galerkin methods for reaction-diffusion equations in unbounded domains. Specifically, we study the fractional version of Allen-Cahn equation, Fisher-Kolmogorov equation, Gray-Scott equation and FitzHugh-Nagumo equation in one-, twoand three-dimensions.

Convergence test.
We first consider a modified one-dimensional Allen-Cahn equation: Here "modified" means we take a different nonlinear term f (u) = F ′ (u), where F (u) is a double-well potential free energy density. In this paper, we study the double-well potential given by as in [23,26]. Note that F (u) has energy minima at u = 0 or 1, different from the usual double-well potential at u = ±1, i.e., F (u) = (u 2 − 1) 2 /4. This is due to the homogenous boundary condition in (1.1) that u(x) → 0, as |x| → ∞. In order to compute with an exact solution, an extra source term g(x, t) is chosen where 1 F 1 a; b; x) is the confluent hypergeometric function defined as 3) The source term (4.2), together with an appropriate initial data, yields an exact solution u(x, t) = exp(−t) exp(−λ 2 x 2 ). This example aims to test the temporal convergence rate for the ETDRK4 scheme established for the fractional reaction-diffusion equation.
In this test, we take ǫ = 0.01, T = 6, λ = 1, N = 500 and several fractional parameters, i.e., s = 0.6, 0.9 and 1.2. In this example, the relative numerical error is computed via Numerical results with respect to time step size τ is given in Fig. 4.1, where the convergence order is shown to be approximately τ 4 , which is in good agreement with the theoretical prediction.
We also present in Fig. 4.1 numerical errors in two-dimensional space with same parameters, where the exact solution is chosen as u(x, y, t) = exp(−t) exp(−(x 2 + y 2 )).
Similar results are obtained.
where r is the intrinsic growth rate of a species and K is the environment carrying capacity, representing the maximum sustainable population density. The Fisher-Kolmogorov equation has attracted a lot of interest due to its wide application to model the growth and spreading of biological species. We first study this equation in one-dimension, and take D = 0.1, r = 0.25, K = 1 and N = 1000. Numerical results are given in Fig. 4.3. The other parameters are the same as in one-dimension. In the left figure of Fig. 4.4, we present the numerical result for the fractional Laplace operator (−∆) s , where s = 0.8. We observe that the contours of the solution are circles since fractional Laplacian is a symmetric operator. We also present in the right figure of Fig. 4.4 the numerical solution for Riesz fractional derivatives in two directions, i.e., we study the following equation

4.3.
Allen-Cahn equation. We consider in this example the spectral-Galerkin method for the fractional-in-space Allen-Cahn equation which is the classical phase field model. Here, the reaction term is given in (4.1). In our first numerical experiment, we set N x = N y = 500, ǫ = 0.02, τ = 0.1. With the random initial data we show the numerical results in Fig. 4.5 for different s. Although the solution domain is infinity, we display in Fig. 4.5 the numerical results in [−1, 1] 2 . In this experiment, we are primarily concerned with the effect of s on the speed of evolution. It is observed from Fig. 4.5 that the fractional Allen-Cahn models with larger values of s produce a rapid movement to larger bulk regions. We also note that reducing the fractional derivative order s leads to a thinner interfacial layer. In next test for fractional Allen-Cahn equation, we study the evolution of mean curvature for Here the initial condition is taken as We set ǫ = 0.01, τ = 10 −4 and N x = N y = 500, numerical results are given in Fig. 4.3.
As time evolves, the radius when u ≥ 0.5 is recorded for s = 0.8. Results in Fig.4.3 (a) shows the evolution of the numerical solution with the initial value given in (4.9) at t = 0, 0.2, 0.5, 0.8. The 0.5 level contour lines of u(x, y, t) are plotted in (b). Here contour of circles are observed since we use the fractional Laplace operator. At the same time, the circle shrinks with time. In (c), quantitative change of radius 2 , i.e., area, with respect to time for different s are given. We know for s = 1, the evolution of the mean curvature flow is radius = R 2 0 − 2t. In our case, here R 0 = 0.6. Thus, singularity for s = 1 happens at t = 0.18, which is the disappearing time, as is also shown in our numerical result in figure (c). In fact, linear decrease between area and time is observed in all cases, while the slope is different. When decreasing s, a smaller slope is obtained.

Gray-Scott equation. The space-fractional Gray-Scott equation is given by
where K u , K v are the two diffusion coefficients and F, κ are the dimensionless positive constants. This model corresponds to the following two reactions where U, V and P represent different chemical species. During the two reactions, U is continuously supplied at a given feed rate F and V is removed at a given kill rate F + κ. It can be easily obtained that (U * , V * ) = (1, 0) is a trivial stable point for this problem.
In order to satisfy the homogeneous boundary conditions as |x| → ∞, we need to modify the reaction terms in (4.10) a little bit. By letting u = 1 − U , v = V , we have for the first equation in (4.10) that i.e., ∂u ∂t Consequently, the Gray-Scott equation (4.10) is turned to (4.14) Now (u * , v * ) = (0, 0) is a stable point of (4.14). In our example, we choose the following initial condition as a permutation of the trivial stable state. Other parameters are K u = 2 × 10 −5 , K v = K u /2, F = 0.03, N = 800, τ = 0.1, while we vary s and κ. Here the ratio between diffusion coefficients K u /K v > 1 is chosen since the model is known to generate different mechanisms of pattern formation depending on the values of feed rate F and kill rate κ.

4.5.
FitzHugh-Nagumo equation. The FitzHugh-Nagumo equation represents one of the simplest models for studying excited media. We study the following space-fractional equation: where u is the excitation variable and v the recovery variable. The FitzHugh-Nagumo equations have been used to qualitatively model many biological phenomena. It is known that when δ = 0, (u * , v * ) = (0, 0) is a stable point for this problem. For our numerical experiments, we choose the initial condition  4.6. Allen-Cahn equation in three-dimensional space. We consider the spectral-Galerkin method for the fractional-in-space Allen-Cahn equation in three spatial dimensions. Here, the reaction term is given in (4.1). In our numerical experiment, we set N x = N y = N z = 200, other parameters are the same as the example for 2D problem in subsection 4.3. With the random initial data u 0 = 0.5 + 0.1(rand − 0.5), (4.19) we show the numerical results for isosurfaces of u(x, y, z) = 0.5 at final time T = 100 in

Concluding remarks
In this work, we applied the efficient and accurate Fourier-like method for fractional reaction-diffusion equations in unbounded domains as proposed in [40]. The method is particularly effective for problems involving the fractional Laplacian. To fully discretize the underlying systems, we propose to use an accurate time marching scheme based on ETDRK4. The main advantage of our method is that it yields a fully diagonal representation of the fractional Laplace operator, which provides great efficiency. Numerical examples are presented to illustrate the effectiveness of the proposed method. We also remark that this method can also be used to deal with problems involving Riesz derivatives in different directions and coupled systems of arbitrary n species.
There are still several possible extensions following the present method, e.g., space-and time-fractional problems involving the fractional Laplacian issues require future investigations. Other types of fractional derivatives such as Riemann-Liouville fractional derivatives and Caputo fractional derivatives will be the topic of our forthcoming research.

Acknowledgement
The work of the author is partially supported by the NSF of China (under the Grant No. 11731006).