Stable finite difference method for fractional reaction–diffusion equations by compact implicit integration factor methods

In this paper we propose a stable finite difference method to solve the fractional reaction–diffusion systems in a two-dimensional domain. The space discretization is implemented by the weighted shifted Grünwald difference (WSGD) which results in a stiff system of nonlinear ordinary differential equations (ODEs). This system of ordinary differential equations is solved by an efficient compact implicit integration factor (cIIF) method. The stability of the second order cIIF scheme is proved in the discrete L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L^{2}$\end{document}-norm. We also prove the second-order convergence of the proposed scheme. Numerical examples are given to demonstrate the accuracy, efficiency, and robustness of the method.


Introduction
The space-fractional reaction-diffusion equations, in which the integer-order differential operator is replaced by a corresponding fractional one, have gained great popularity due to the wide application in physics [1], finance [2], and biology [3][4][5][6]. These new fractionalorder models provide a more adequate description of many processes than those of the integer-order for many processes with anomalous diffusion. This is because factionalorder operators possess the nonlocal nature property such that the models can describe the phenomena showing the effect of memory and hereditary properties [7].
In this paper, we consider the following two-dimensional space-fractional reactiondiffusion equation where 1 < α x , α y ≤ 2 are the space fractional orders, K x , K y > 0 are the diffusion constants. The Riesz fractional derivatives ∂ αx u ∂|x| αx and ∂ αy u ∂|y| αy are defined as follows: where a D α x x and x D α x b are left and right Riemann-Liouville fractional derivatives defined as Here (·) denotes the standard Gamma function.
Owing to the applications of fractional differential equations in many fields, it is important to seek effective method for these fractional models. The fractional differential equations are solved and analyzed by various analytical methods, such as monotone iterative method [8], Green function method [9], Laplace transform method [10], homotopy perturbation transform method [11], and other methods [12][13][14][15][16]. However, the analytic solutions of special fractional differential equations are in the form of trigonometric series, and it is very difficult to calculate them. For this reason, the study for the numerical methods has become an important research topic that offers great potential.
In recent years, various classical numerical methods have been extended to solve the fractional reaction-diffusion equations, such as the finite difference method (FDM) [17][18][19][20], finite element method (FEM) [21][22][23][24], and spectral method [25][26][27]. The finite difference methods are accepted as one of the most popular numerical methods for fractional reaction-diffusion equations because they allow easy implementation. Meerschaert and Tadjeran proposed the first-order shifted Grünwald formula in [18]. Based on the former work, Tian et al. developed weighted shifted Grünwald-Letnikov (WSGD) formulas [20]. Ortigueira firstly proposed the so-called fractional central difference scheme in [19]. Then Çelik and Duman analyzed this approximation and applied it to fractional diffusion equations [17].
The space discretization of fractional reaction-diffusion equation leads to a nonlinear ordinary differential equations (ODEs). Accurate and efficient simulations of such systems have lead many researcher to the design of various time discretization methods. The most notable method among them is the alternating direction implicit (ADI) technique. However, the discretization matrix obtained is full. It means that the Thomas algorithm for tridiagonal matrix is no longer applicable. It also should be noted that the ADI methods are limited to second-order accuracy in time. In this paper, an efficient compact implicit integration factor (cIIF) method [28,29] is developed. By introducing the compact representation for the matrix approximating the differential operator, the cIIF methods apply matrix exponential operations sequentially in every spatial direction. As a result, exponential matrices which are calculated and stored have small sizes, compared to those in the 1D problem. For two or three dimensions, the cIIF method is more efficient in both storage and CPU cost.
The main contributions of this work are as follows: • We provide an efficient approach to compute solutions of the two-dimensional fractional reaction-diffusion equations even with millions of difference points. • We prove that the second-order cIIF scheme is stable in the discrete L 2 -norm.
• We prove the second-order convergence of the proposed scheme, and the numerical tests verify the theoretical analysis. The rest of the paper is organized as follows: In Sect. 2 we present some notations and discretize the two-dimensional fractional reaction-diffusion equations into nonlinear ODEs in matrix formulation. In Sect. 3, we present the cIIF time discretization scheme. In Sect. 4 we compute various fractional reaction-diffusion equations to demonstrate the convergence rates and the efficiency of the proposed scheme. Finally, we summarize our conclusion in Sect. 5.

Numerical method 2.1 Compact finite difference method
The computation domain is discretized into grids described by the set (x j , We first introduce the second-order WSGD method derived in [20] to approximate the Riemann-Liouville space fractional derivatives. The essential idea of this approximation is using the weighted average to eliminate the low-order leading terms in asymptotic expansions for the truncation errors. The WSGD formulas for the left and right Riemann-Liouville fractional derivatives in the x-direction are defined as for j = 1, . . . , N x -1, k = 1, . . . , N y -1. The coefficients ω (α) k of the second order in (5) are given as with We refer the reader to the paper [20] for the detailed proof of the truncation error expression.
Based on the discretization of Riemann-Liouville fractional derivatives, the Riesz fractional derivatives in the x-direction can be discretized as where Similarly, Riesz fractional derivatives in the y-direction can be discretized as follows: where δ α y u(x j , y k ) = -1 . Defining the grid function U j,k (t) = u(x j , y k , t), the semidiscrete compact difference scheme for problem (1) can be written as denotes the truncation error in space. Omitting the truncation term R j,k (t) and replacing the grid function U j,k (t) with its numerical approximation u j,k (t), we obtain the following difference scheme: To apply the time discretization method for the ODEs (10), we now rewrite them in a matrix form. Define the numerical solution in the following matrix forms: For a solution matrix U ∈ R N x -1×N y -1 , the discrete L 2 -norm of U is defined as Then we have where U F is Frobenius norm. The 2-norm for matrix U is defined as U 2 = σ max (U) where σ max (U) represents the largest singular value of matrix U. In the case of a normal matrix U, U 2 = λ max (U), with λ max (U) being the absolute value of the largest eigenvalue of U.
Then introducing the M × M difference matrix D α M as follows: the difference scheme (8) and (9) can be written in the following matrix form: Let K x D α N x -1 = A and K y D α N y -1 = B. Then the ODEs (10) can be written in the following matrix form: where F(U) is an (N x -1) × (N y -1) matrix with each element defined as F(u j,k ).

Compact implicit integration factor method
Now we will apply the cIIF method to solve (15). Assume the final time is t = T and let time step t = T/N , t n = n t, 0 ≤ n < N . The numerical approximation matrix for u j,k (t) at time t n is defined as U n . To construct the cIIF methods for (15), we multiply it by the integration factor e -At from the left, and by e -Bt from the right, and then integrate over one time step from t n to t n+1 to obtain Then we approximate the integrand in (16) by using the (r -1)th order Lagrange interpolation polynomial with interpolation points at t n+1 , t n , . . . , t n-r+2 , namely where The local truncation error is O( t r ). In this paper we use the following second-order compact implicit integration factor (cIIF2) scheme: The local truncation error for (18) can be written in matrix form as E = (E jk ) N x -1×N y -1 with each element being E jk = O( t 3 + th 2 ).

Stability and convergence
We will study the linear stability of the numerical scheme for equation (18) with reaction term F(u) satisfying the Lipschitz continuity condition. To obtain the linear stability analysis, we need the following lemmas. (13) is a symmetric and negative definite matrix.

Lemma 3.2 Let
A and B be arbitrary N × N square matrices. We have the following inequalities: Proof Taking square roots completes the proof. (1) is globally Lipschitz continuous, i.e., there exists a real constant L such that |f (u)f (v)| ≤ L|u -v|. The fourth-order cFDM coupled with cIIF2 scheme (18) is unconditionally stable in the discrete L 2 -norm.

Theorem 3.3 Assume that function f (u) in
Proof Assume that U 1 n+1 and U 2 n+1 are two different solutions of (18) with different initial data sets. Let n+1 = U 1 n+1 -U 2 n+1 , then it satisfies By Lemma 3.2 and Lipschitz continuity, taking Frobenius norm for (19) Since the matrices A and B are negative definite, the eigenvalues of matrix e A t are smaller than 1 such that e A t 2 ≤ 1 and e B t 2 ≤ 1. We get the following result from (20): With time iterations, we have From the above result (22), it can be concluded that Then (23) yields the following result: Now we complete the proof by (12). (1) is globally Lipschitz continuous. The fourthorder cFDM coupled with cIIF2 scheme (18) converges to the solution of fractional problem

Theorem 3.4 Assume that function f (u) in
Proof We denote by U and U the exact value and the approximation of solution matrix of (1), respectively. Substituting U into (18), we have where E = (E jk ) N x -1×N y -1 and E jk = O( t 3 + th 4 ). Denote = U -U. Subtracting (25) from (18) yields the following error equation: It is obvious that 0 = 0 for the initial condition. Taking Frobenius norm for (26) and applying Theorem 3.3, we obtain by Lipschitz continuity. Here C denotes a positive constant independent of discretization parameters, which may take different values at different occurrences. As in the proof of Theorem 3.3, we have As time evolves, we have Multiplying by h x h y on both sides of the above inequality, we get . This completes the proof.

Numerical experiments
In this section, we will demonstrate the performance of the proposed scheme on some test problems. Firstly, in order to verify our theoretical analysis, we test our scheme for a nonlinear reaction-diffusion equation with exact solution. Then we apply the scheme to the two-dimensional space Riesz fractional Fitzhugh-Nagumo model which represents one of the simplest models for studying excitable media. In all the numerical experiments, we used a uniform spatial step size along each direction, i.e., h x = h y = h. All the computations were performed in Matlab based on a ThinkPad desktop with i3-3110 CPU and 4 GB memory.
Example 1 We consider the following nonlinear reaction-diffusion equation on a square = [0, 1] 2 : with F(u) = -u 2 . The exact solution of this equation is u = e -t x 2 (1x) 2 y 2 (1y) 2 . The final computation time is T = 1. The function f (x, y, t) is determined by the exact solution. For details of f (x, y, t), please refer to [31].
We use the second-order WSGD method coupled with cIIF scheme to demonstrate the accuracy of space and time discretization. The orders of convergence in space and time are computed as q = log 2 (e(h)/e(h/2)) and p = log 2 (e(τ )/e(τ /2)), respectively. Table 1 displays the L 2 errors of the WSGD scheme with different values of α x and α y . The order of convergence is computed using a very small time step t = 1E -3. In Table 2, the temporal errors and convergence order of the cIIF2 scheme are given for different time steps and α x and α y , with h = 1/256. From Tables 1 and 2, we conclude that the convergence rate in space and time is of second order. The numerical results are well in line with the theoretical analysis.
with homogeneous Dirichlet boundary condition.
The Fitzhugh-Nagumo model can generate the stable patterns in the form of spiral waves. Initially, we set the state in the lower-left quarter of the domain as u = 1 and in the upper-half part as w = 0.1. Then the trivial state (u, w) = (0, 0) was perturbed and further curved and rotated clockwise, generating the spiral pattern.
In our simulation, the computation domain is discretised using N x = N y = 256 points in each spatial coordinate. The time step is chosen as t = 0.2 and the final time is set to be T = 1000. We first consider the integer-order Fitzhugh-Nagumo model (i.e., α x = α y = 2). The spiral wave solutions for two diffusion coefficients K x = K y = 1E -4 and K x = K y = 1E -5 are shown in Fig. 1. We find that the number of wavefronts increases with decreasing the diffusion coefficient. Keeping the coefficients K x = K y = 1E -4 and reducing the fractional powers α x and α y , we then get the stable rotating solution for α x = α y = 1.8 ( Fig. 2(a)) and α x = α y = 1.5 ( Fig. 2(b)). As expected, the width of the excitation wavefront (red areas) reduces with decreasing α x and α y , so does the number of the wavefronts. However, it is noted that the role of reducing the fractional power is not equivalent to the influence of a decreased diffusion coefficient in the integral order case. This can be clearly observed by comparing Figs. 1 and 2.
We take the different diffusion coefficients in the x-and y-directions as K x = 1E -4, K y = 2.5E -5. Figure 3(a) displays the wave propagation at T = 1000. The smaller diffusion coefficient in the y-direction leads to an elliptical pattern. Next we consider the anisotropic fractional power α x = 2, α y = 1.7 with K x = K y = 1E -4. A similar pattern is displayed in

Concluding remarks
In this paper, the weighted shifted Grünwald difference method is developed to solve the two-dimensional Riesz space fractional nonlinear reaction-diffusion equation, in which the temporal component is discretized by the compact implicit integration factor method. The advantage of the method is that computing exponentials of large matrices is reduced to the computation of exponentials of matrices of significantly smaller sizes. The stability and convergence are strictly proven, which shows that the compact implicit integration factor method is stable in L 2 -norm and convergent with second-order accuracy. The WSGD method coupled with cIIF method is applied to solve the fractional Fitzhugh-Nagumo model. Numerical experiments are provided to verify the theoretical analysis. The numerical results confirm that the proposed method is a powerful and reliable method for the fractional reaction-diffusion equations. Given its efficiency and good stability conditions, the cIIF method can be extended to the fourth-order fractional parabolic equations (e.g., Cahn-Hilliard equations), which will also be further explored in the future work.