A stabilized semi-implicit Fourier spectral method for nonlinear space-fractional reaction-diffusion equations

TL;DR

This paper introduces a second-order stabilized semi-implicit Fourier spectral method for solving nonlinear space-fractional reaction-diffusion equations, ensuring stability without CFL constraints and providing optimal error estimates.

Contribution

The work develops a novel stabilized semi-implicit spectral method for fractional reaction-diffusion equations with proven stability and error bounds, applicable to complex biological and physical models.

Findings

Method is unconditionally stable and second-order accurate.

Successfully applied to fractional Allen-Cahn, Gray-Scott, and FitzHugh-Nagumo models.

Revealed distinct pattern formation behaviors in fractional models.

Abstract

The reaction-diffusion model can generate a wide variety of spatial patterns, which has been widely applied in chemistry, biology, and physics, even used to explain self-regulated pattern formation in the developing animal embryo. In this work, a second-order stabilized semi-implicit time-stepping Fourier spectral method is presented for the reaction-diffusion systems of equations with space described by the fractional Laplacian. We adopt the temporal-spatial error splitting argument to illustrate that the proposed method is stable without imposing the CFL condition, and we prove an optimal L2-error estimate. We also analyze the linear stability of the stabilized semi-implicit method and obtain a practical criterion to choose the time step size to guarantee the stability of the semi-implicit method. Our approach is illustrated by solving several problems of practical interest, including the fractional Allen-Cahn, Gray-Scott and FitzHugh-Nagumo models, together with an analysis of the properties of these systems in terms of the fractional power of the underlying Laplacian operator, which are quite different from the patterns of the corresponding integer-order model.

Tables

Table 1. Table 1: Case I of Example 5.1 : the convergence rate in time for different values of κ 𝜅 \kappa .

$\alpha$	$\tau$	$\kappa =1$		$\kappa =3$		$\kappa =4$
		Error	Order	Error	Order	Error	Order
	1/10	5.2105e-3	-	1.3160e-2	-	1.7031e-2	-
	1/20	1.3441e-3	1.9548	3.4301e-3	1.9398	4.4529e-3	1.9354
1.3	1/40	3.4169e-4	1.9759	8.7904e-4	1.9643	1.1447e-3	1.9598
	1/80	8.5820e-5	1.9933	2.2188e-4	1.9862	2.8956e-4	1.9830
	1/160	2.1142e-5	2.0212	5.4810e-5	2.0173	7.1619e-5	2.0155
	1/10	3.1801e-3	-	8.1250e-3	-	1.0553e-2	-
	1/20	8.0715e-4	1.9782	2.0745e-3	1.9696	2.6943e-3	1.9697
1.7	1/40	2.0396e-4	1.9845	5.2881e-4	1.9720	6.8893e-4	1.9675
	1/80	5.1102e-5	1.9968	1.3329e-4	1.9881	1.7410e-4	1.9844
	1/160	1.2576e-5	2.0228	3.2914e-5	2.0178	4.3059e-5	2.0155

Table 2. Table 2: Case II of Example 5.1 : the convergence rate in time for different values of κ 𝜅 \kappa .

$\alpha$	$\tau$	$\kappa =1$		$\kappa =3$		$\kappa =4$
		Error	Order	Error	Order	Error	Order
	1/10	1.6244e-2	-	3.6935e-2	-	4.6592e-2	-
	1/20	4.3963e-3	1.8856	1.0319e-2	1.8397	1.3196e-2	1.8200
1.3	1/40	1.1374e-3	1.9506	2.7041e-3	1.9321	3.4785e-3	1.9235
	1/80	2.8767e-4	1.9832	6.8752e-4	1.9757	8.8665e-4	1.9720
	1/160	7.1082e-5	2.0168	1.7028e-4	2.0135	2.1984e-4	2.0119
	1/10	1.6522e-2	-	3.7636e-2	-	4.7481e-2	-
	1/20	4.4665e-3	1.8871	1.0508e-2	1.8406	1.3442e-2	1.8206
1.7	1/40	1.1550e-3	1.9513	2.7525e-3	1.9327	3.5420e-3	1.9241
	1/80	2.9206e-4	1.9835	6.9973e-4	1.9759	9.0270e-4	1.9723
	1/160	7.2160e-5	2.0170	1.7329e-4	2.0136	2.2381e-4	2.0120

Equations

$$\displaystyle\left\{\begin{array}[]{ll}{\partial_{t}}u(x,t)=-K_{\alpha}(-\Delta)^{\alpha/2}u(x,t)+G(u(x,t)),\ (x,t)\in\Omega\times I,\\ u(x,0)=u^{0}(x),\ x\in\Omega,\\ u(\cdotp,t)\mathrm{\ is\ 2\pi-periodic},\ x\in\partial\Omega,\ \forall t\in I,\\ \end{array}\right.$$

$$(-\Delta {)}^{\alpha /2}u(x)=\sum _{k,l\in \mathbb{Z}}(k^2+l^2{)}^{\alpha /2}{\hat{u}}_{kl}{e}^{ik{x}_1+il{x}_2},$$

$${\hat{u}}_{kl}=(u,e^{ik{x}_1+il{x}_2})=\frac{1}{(2\pi {)}^2}{\int }_{\Omega }u{e}^{-(ik{x}_1+il{x}_2)}dx.$$

$$u(x)=\sum _{k,l\in \mathbb{Z}}{\hat{u}}_{kl}{e}^{ik{x}_1+il{x}_2}.$$

$$|u{|}_r^2=\sum _{k,l\in \mathbb{Z}}{\hat{u}}_{kl}^2(k^2+l^2{)}^r,\Vert u{\Vert }_r^2=\sum _{k,l\in \mathbb{Z}}{\hat{u}}_{kl}^2(1+k^2+l^2{)}^r.$$

$$X_N=\mathrm{span}\left\{e^{ik{x}_1+il{x}_2}:-N/2\le k,l\le N/2-1\right\}.$$

$$(u-P_{N}u,v)=0,\forall v\in X_N,u\in L_{per}^2(\Omega ).$$

$$D_1{u}^n=\frac{u^n-u^{n-1}}{\tau },D_2{u}^n=\frac{3}{2}{D}_1{u}^n-\frac{1}{2}{D}_1{u}^{n-1}.$$

$$D_2{u}^n+K_{\alpha }(-\Delta {)}^{\alpha /2}{u}^n=G(u^n)+O({\tau }^2),n\ge 2.$$

$$D_2{u}^n+K_{\alpha }(-\Delta {)}^{\alpha /2}{u}^n=2G(u^{n-1})-G(u^{n-2})-\tau \kappa (D_1{u}^n-D_1{u}^{n-1})+R^n,2\le n\le K,$$

$$\{\begin{array}{rl}& {\mathcal{A}}^n(U_N,v)=(P_N(2G(U_N^{n-1})-G(U_N^{n-2})),v)-\tau \kappa (D_1{U}_N^n-D_1{U}_N^{n-1},v),\\ & U_N^0=P_N{u}^0,\\ & U_N^1=P_N{u}^1,\end{array}$$

$${\mathcal{A}}^n(u,v)=(D_2{u}^n,v)+K_{\alpha }((-\Delta {)}^{\alpha /2}{u}^n,v).$$

$$U_N^1=P_N(u_0+\tau {\partial }_{t}u(0)).$$

$$\Vert u^n-U_N^n\Vert \le C({\tau }^2+N^{-r}).$$

$${\mathcal{A}}^n(U_N,v)=\rho (2U_N^{n-1}-U_N^{n-2},v)-\kappa (U_N^n-2U_N^{n-1}+U_N^{n-2},v).$$

$$\left(\frac{3}{2}+{\mu }_{k,l}\tau +\kappa \tau \right){\hat{U}}_{k,l}^n-\left(2+2\rho \tau +2\kappa \tau \right){\hat{U}}_{k,l}^{n-1}+\left(\frac{1}{2}+\rho \tau +\kappa \tau \right){\hat{U}}_{k,l}^{n-2}=0,$$

$$\left(\frac{3}{2}+{\mu }_{k,l}\tau +\kappa \tau \right){\eta }^2-\left(2+2\rho \tau +2\kappa \tau \right)\eta +\left(\frac{1}{2}+\rho \tau +\kappa \tau \right)=0.$$

$${\eta }^{\pm }=\frac{2+2\rho \tau +2\kappa \tau \pm \sqrt{△}}{3+2{\mu }_{k,l}\tau +2\kappa \tau },$$

$$△={\left(2+2\rho \tau +2\kappa \tau \right)}^2-\left(3+2{\mu }_{k,l}\tau +2\kappa \tau \right)\left(1+2\rho \tau +2\kappa \tau \right).$$

$$-1<\frac{1+2\rho \tau +2\kappa \tau }{3+2{\mu }_{k,l}\tau +2\kappa \tau }<1⟹\kappa >-\frac{{\mu }_{k,l}+\rho }{2}-\frac{1}{\tau }.$$

$$4\kappa \tau +3\rho \tau +{\mu }_{k,l}\tau +4>0⟹\kappa >-\frac{{\mu }_{k,l}+3\rho }{4}-\frac{1}{\tau }.$$

$$\kappa >\frac{-{\mu }_{k,l}-3\rho }{4}-\frac{1}{\tau }.$$

$$\kappa >-\frac{3\rho }{4}.$$

$$\kappa >-\frac{3\rho }{4}-\frac{1}{{\tau }^{\ast }}.$$

$$e^n=\Vert U_N^n-u(\cdot ,t_n)\Vert .$$

$$\displaystyle\left\{\begin{array}[]{ll}{\partial_{t}}u=-K_{u}(-\Delta)^{\alpha/2}u+G_{1}(u,v)\quad(x,t)\in\Omega\times I,\\ {\partial_{t}}v=-K_{v}(-\Delta)^{\alpha/2}v+G_{2}(u,v),\quad(x,t)\in\Omega\times I.\\ \end{array}\right.$$

$$\{\begin{array}{rl}& {\mathcal{A}}_1^n(U_N,v)=(P_N(2G_1(U_N^{n-1})-G_1(U_N^{n-2})),v)-\tau \kappa (D_1{U}_N^n-D_1{U}_N^{n-1},v),\forall v\in X_N,\\ & {\mathcal{A}}_2^n(V_N,v)=(P_N(2G_2(V_N^{n-1})-G_2(V_N^{n-2})),v)-\tau \kappa (D_1{V}_N^n-D_1{V}_N^{n-1},v),\forall v\in X_N,\\ & U_N^0=P_N{u}^0,U_N^1=P_N(u_0+\tau {\partial }_{t}u(0)),\\ & V_N^0=P_N{v}^0,V_N^1=P_N(v_0+\tau {\partial }_{t}v(0)),\end{array}$$

$${\mathcal{A}}_1^n(u,v)=(D_2{u}^n,v)+K_u((-\Delta {)}^{\alpha /2}{u}^n,v),$$

$${\mathcal{A}}_2^n(u,v)=(D_2{u}^n,v)+K_v((-\Delta {)}^{\alpha /2}{u}^n,v).$$

$$G_1(u,v)=-u{v}^2+F(1-u),G_2(u,v)=u{v}^2-(F+\lambda )v,$$

$$U+2V$$

Full text

1 Introduction

We consider the following two-dimensional nonlinear space-fractional reaction-diffusion equation [2, 3, 4],

[TABLE]

where $\Omega=(0,2\pi)^{2}$, $I=(0,T]$, $x=(x_{1},x_{2})$, $K_{\alpha}>0$ is the diffusion coefficient, $G(u)$ is a nonlinear function that satisfies the local Lipschitz condition, and $(-\Delta)^{\alpha/2}$ $(1<\alpha\leq 2)$ is the fractional Laplacian operator defined by [5]

[TABLE]

in which $i^{2}=-1$ and $\hat{u}_{kl}$ are the Fourier coefficients,

[TABLE]

For any function $u(x)\in L_{per}^{2}(\Omega)$, we have [6, 7]

[TABLE]

The fractional reaction-diffusion equation has attracted much attention in the last few years due to its wide applications in many fields of science and engineering [4, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. For $G(u)=u-u^{3}$, (1) becomes the space-fractional Allen–Cahn equation [12], which describes the mixture of two incompressible fluids. The fractional order controls the sharpness of the interface, which is typically diffusive in integer-order phase-field models. As shown in [19], the space-fractional Gray–Scott model is a coupled counterpart of (1), which represents an autocatalytic reaction-diffusion process between two chemical species [18]. The space-fractional FitzHugh-Nagumo model describes the ionic fluxes as a function of the membrane potential [20].

Numerical methods for space-fractional reaction-diffusion equations similar to (1) have been extensively investigated. Ervin et al. [21] considered the numerical approximation for space-fractional diffusion equations containing a nonlocal quadratic nonlinearity, in which an optimal $H^{\alpha}$-error estimate was presented under the time step restriction $\tau\leq ch$, where $\tau$ is the temporal step size and $h$ is the spatial mesh size. An alternating direction implicit Galerkin-Legendre spectral method for the two-dimensional Riesz space fractional nonlinear reaction-diffusion equation was proposed in [22], which requires $\tau^{2}N\leq c$ to ensure stability or convergence, where $N$ is the degree of polynomial spatial basis. Similar results were obtained in [23, 24, 25, 26], in which numerical methods for integer-order partial differential equations were proposed and analyzed. The conditions as $\tau\leq ch$ and $\tau^{2}N\leq c$ were required to prove the stability and convergence of the aforementioned works due to the use of the inverse inequality in the numerical analysis. However, the conditions as $\tau\leq ch$ and $\tau^{2}N\leq c$ in [22, 23, 24, 25, 26] can be weakened by a temporal-spatial error splitting argument [27]. In this work, we extend the temporal-spatial error splitting argument to analyze the stabilized semi-implicit Fourier spectral method for the two-dimensional nonlinear space-fractional reaction-diffusion equations, which has not been fully studied. Stabilized time-stepping methods have been applied for solving nonlinear integer-order partial differential equations, such as Allen-Cahn and Cahn-Hilliard equations, see e.g. [28, 29, 30, 31].

The main contribution of this work is to prove that the stabilized semi-implicit method is stable without imposing the CFL condition by adopting the temporal-spatial error splitting argument [27, 32, 33], and present an optimal error estimate. We also analyze the linear stability of the stabilized semi-implicit method and obtain the explicit form of the stability condition criterion (see (20) and (21)), which is verified by numerical simulations. The temporal-spatial error splitting technique was originally developed for numerical analysis of integer-order equations. Up to now, this method has been applied to the numerical analysis of time-fractional partial differential equations [34, 35]. In this paper, we extend this method to the two-dimensional nonlinear space-fractional partial differential equations, which has not been fully explored. Our approach is illustrated by solving several problems of practical interest, including the fractional Allen-Cahn, Gray-Scott and FitzHugh-Nagumo models.

This paper is organized as follows. Section 2 develops a second-order stabilized semi-implicit time-stepping Fourier spectral method to solve the two-dimensional nonlinear space-fractional reaction-diffusion equation. In Section 3, we analyze the linear stability of the stabilized semi-implicit method and obtain a practical criteria to choose the time step size to guarantee the stability of the semi-implicit method. The semi-implicit method is extended to a system of space-fractional reaction-diffusion models in Section 4. Numerical examples are given in Section 5 to confirm the theoretical analysis. Finally, the conclusions are summarized in Section 6.

2 The stabilized semi-implicit Fourier method

Let $C_{per}^{\infty}(\Omega)$ be the set of all restrictions onto $\Omega$ of all complex-valued, $2\pi$-periodic, $C^{\infty}$-functions on $\mathbb{R}^{2}$. For a nonnegative real number $r$, let $H_{per}^{r}(\Omega)$ be the closure of $C_{per}^{\infty}(\Omega)$ with the norm $\|\cdot\|_{r}$ and semi-norm $|\cdot|_{r}$ defined by (see [5])

[TABLE]

For simplicity of the numerical analysis, we also assume that $\hat{u}_{00}=0$, which implies that the norm $\|\cdot\|_{r}$ and semi-norm $|\cdot|_{r}$ are equivalent.

For a positive integer $N$, the function space is denoted by

[TABLE]

Define the orthogonal projection operator $P_{N}$ as follows

[TABLE]

For the temporal discretization, we divide the interval $[0,T]$ into $K$ subintervals with a time step size $\tau=T/K$. Let $t_{n}=n\tau$, $u^{n}=u(x,t_{n})$, $0\leq n\leq K$, and denote

[TABLE]

The time derivative is discretized by the second-order backward difference formula, we obtain

[TABLE]

Dropping the truncation error and applying the Fourier spectral method (other methods can be applied, i.e., finite difference/element method) in space, we can obtain a fully implicit numerical scheme, which can be solved iteratively. Compared with the explicit method or the linearized method, the fully implicit method is costly and its implementation is a little complicated. The extrapolation has been widely used to obtain the semi-implicit method, which needs to solve a linear system. For example, The second-order extrapolation $G(u^{n})=2G(u^{n-1})-G(u^{n-2})+O(\tau^{2})$ can be applied, which preserves second-order accuracy in time, but it needs a small time step to ensure stability. In order to balance the stability and convergence, we follow [28] to develop the stabilized method for the time discretizetion, that is, we replace $G(u^{n})$ with $2G(u^{n-1})-G(u^{n-2})-\tau\kappa(D_{1}u^{n}-D_{1}u^{n-1})$ to obtain

[TABLE]

where $R^{n}=O(\tau^{2})$, $\kappa\geq 0$, and $\tau\kappa(D_{1}u^{n}-D_{1}u^{n-1})=\kappa(u^{n}-2u^{n-1}+u^{n-2})=O(\tau^{2})$ is a penalty term that balances the stability and accuracy.

From (8), the fully discrete Fourier spectral method for (1) is given by: Find $U_{N}^{n}\in X_{N}$ for $n\geq 2$ such that for all $v\in X_{N}$,

[TABLE]

where

[TABLE]

Remark 2.1**.**

In real applications, $U_{N}^{1}$ can be obtained by the following simple approach

[TABLE]

*Other approaches can be applied. For example, $U_{N}^{1}$ can be obtained by the fully implicit Crank–Nicolson Fourier spectral method. *

An optimal error estimate of the numerical scheme (9) is provided in the following theorem, the relevant proof can be seen in Appendix A.

Theorem 2.1**.**

Suppose that (1) admits a unique solution $u(\cdot,t)\in H^{r}(\Omega)$ satisfying (35). Let $U_{N}^{n}(2\leq n\leq K)$ be the solution of (9) and $G^{\prime}(z)$ is bounded for $|z|\leq M$, $M$ is a suitable positive constant. Then there exist two positive constants $\tau^{*}$ and $N^{*}$, when $\tau\leq\tau^{*}$ and $N\geq N^{*}$, it holds

[TABLE]

3 Linear stability

We investigate the linear stability of the method (9) in this section. Let $G(u)=\rho u$ and $\rho<0$, we have

[TABLE]

The above equation can be equivalently written as

[TABLE]

where $\mu_{k,l}=K_{\alpha}(k^{2}+l^{2})^{\alpha/2}$. The characteristic equation of the above recurrence formula is

[TABLE]

The two roots of (14) are

[TABLE]

where

[TABLE]

The recurrence relation (13) is stable if $|\eta^{\pm}|<1$. We discuss the following two cases.

Case I: For $\triangle<0$ and $|\eta^{\pm}|<1$, we have

[TABLE]

• 2.

Case II: For $\triangle>0$ and $|\eta^{\pm}|<1$, we have

[TABLE]

Combining (17)–(LABEL:e33.0356), we can obtain that the method (12) is stable if

[TABLE]

Since $\mu_{k,l}\geq 0$ and $\tau>0$, we have the following unconditional stability condition

[TABLE]

In real applications, we need to choose $\tau\leq\tau^{*}<1$ in the numerical simulations to obtain suitably accurate numerical solutions. From (19), for the step size $\tau\in(0,\tau^{*}]$, the method (12) is stable if

[TABLE]

In the following example, we verify the stability criterion (21).

Example 3.1**.**

We apply the method (9) to solve (1), where $G(u)=-8u+\exp(-t)(-\sin{(x_{1}+x_{2})}+8\sin{(x_{1}+x_{2})}+2^{\alpha/2}\ K_{\alpha}\sin{(x_{1}+x_{2})})$. The exact solution of (1) is $u=\exp(-t)\sin{(x_{1}+x_{2})}.$

We consider the following two cases:

Case I: $\kappa=2$, $\tau=2^{-j},1\leq j\leq 5,T=5$.

• 2.

Case II: $\tau=1/4$, $\kappa=1/2,1,2,4,8$, $T=5$.

The error is defined by

[TABLE]

From (21), we derive that the method is stable for this example if $\kappa>6-1/\tau^{*}$ and $\tau\leq\tau^{*}$. For Case I, we conclude that the method is stable when $\tau\leq\tau^{*}=1/4$ and is unstable when $\tau>\tau^{*}=1/4$, which is also verified in Fig. 1 (a). For Case II, we fix the step size $\tau=1/4$. It is easy to obtain that the method is stable if $\kappa>2$, which is also demonstrated from Fig. 1 (b). We have also verified that the stability condition (21) is independent of the fractional order $\alpha\in(1,2]$, but these results are not shown here.

4 Extension to a system of space-fractional reaction-diffusion models

We extend the semi-implicit method to the following coupled system,

[TABLE]

Each equation of system (22) is discretized by the second-order stabilized semi-implicit time-stepping Fourier spectral scheme as (9). The fully numerical approximation for (22) is developed as follows: Find $U_{N}^{n},V_{N}^{n}\in X_{N}$ for $n=2,3,\ldots,K$ such that

[TABLE]

where

[TABLE]

In (22), if we take

[TABLE]

then we obtain the space-fractional Gray–Scott model.

Reaction and diffusion of chemical species can produce a variety of patterns. A representative reaction and diffusion model is the Gray–Scott model, which is a variant of the autocatalytic Selkov model of glycolysis [36, 37]. Numerical simulations of this model were performed in [38] to find stationary lamellar patterns like those observed in earlier laboratory experiments on iodate-ferrocyanid-sulfite reactions [39]. The Gray–Scott model corresponds to the following two reactions [38],

[TABLE]

both reactions are irreversible, so $P$ is an inert product. A nonequilibrium constraint is represented by a feed term for $U$. Both $U$ and $V$ are removed by the feed process. Note that, in the aforementioned literature on pattern dynamics of the Gray–Scott model, the models are all with standard diffusion, i.e., the diffusion operator is the normal Laplacian. When we study the effects of the superdiffusion, we have to employ the fractional Laplacian $(-\Delta)^{\alpha/2}$ for $1<\alpha\leq 2$ on pattern formation of this model. This space-fractional Gray–Scott model describes an autocatalytic reaction-diffusion process between two chemical species with concentrations $u$ and $v$ [18].

Taking $K_{v}=0$ and

[TABLE]

we can get the space-fractional FitzHugh–Nagumo model.

Spiral waves in excitable media provide an important example of self-organization phenomena in spatially distributed biological, physical, and chemical systems [40], and can be modeled by the FitzHugh–Nagumo equations [41, 42]. The fractional FitzHugh–Nagumo model consists of a space-fractional nonlinear reaction-diffusion model and a system of ordinary differential equations, describing the ionic fluxes as a function of the membrane potential [18]. Mathematical models of electrical activity in cardiac tissue are becoming increasingly powerful tools in the study of cardiac arrhythmias. In this context, the fractional model presented here represents a new approach to dealing with the propagation of the electrical potential in heterogeneous cardiac tissue. In [18], the authors proposed a fundamental rethinking of the homogenization approach via the use of a fractional Fick’s law [43, 44] and, in particular, a fractional FitzHugh–Nagumo model is introduced to capture the spatial heterogeneities and spatial connectivities in the extracellular domain through the use of fractional derivatives. A semi-alternating direction method for a 2-D fractional FitzHugh-Nagumo monodomain model on an approximate irregular domain was considered in [20].

5 Numerical results

In this section, we provide some numerical results to confirm our theoretical analysis. The error at $t=t_{n}$ is denoted by

[TABLE]

where $u_{ref}$ is the reference solution obtained using a smaller time step size $\tau=0.0005$ for a fixed $N$, or a larger $N=1024$ for a fixed time step size if the exact solution is not explicitly given. The $L^{2}$-norm convergence orders in time and space are defined as

[TABLE]

where $\tau_{1}\neq\tau_{2}$ and $N_{1}\neq N_{2}$.

Example 5.1**.**

The Allen-Cahn equation with a quartic double well potential is a simple nonlinear reaction-diffusion model that arises in the study of formation and motion of phase boundaries [45, 46]. The space-fractional equation takes the form [12, 47]

[TABLE]

where $K_{\alpha}$ is a small positive constant.

Case I: Choose the following initial condition

[TABLE]

We choose the values of the parameters $K_{\alpha}=0.01$. For Case I, Table 1 presents the $L^{2}$-errors for different values of $\kappa$ at $t=2$ when $N=32$. It is obvious that second-order accuracy is observed, which is in line with the theoretical analysis. Next, we fix the step size $\tau=0.01$, $\kappa=1$, the $L^{2}$-errors at $t=2$ against $N$ are shown in Fig. 2 (a), spectral accuracy is displayed, which agrees with the theoretical analysis.

Case II: Choose the following initial condition

[TABLE]

Considering $K_{\alpha}=0.01$, $N=32$, we present the $L^{2}$-error for different values of $\kappa$ at $t=2$ and its convergence order in time in Table 2. Clearly, second-order accuracy is still observed for Case II. Fixing $\tau=0.01$, Fig. 2 (b) shows spectral accuracy in space for a fixed $\tau=0.01$, $\kappa=1$.

Example 5.2** (Fractional Gray–Scott model).**

For the fractional Gray–Scott model (see (22) and (LABEL:eq4.0)), $u$ and $v$ are the concentrations of $U$ and $V$ (see (25)), respectively, $K_{u}$ and $K_{v}$ are the diffusion rates in the process, $F$ is the dimensionless feed rate, and $\lambda$ is the dimensionless rate constant of the second reaction. In order to illustrate the effect of varying $\alpha$ in the model, we take the spatially initial condition as

[TABLE]

where $\Omega=(-1,2)^{2}$ and $O_{c}=\{(x_{1},x_{2})|(x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}\leq 0.04^{2}\}$.

Take $\kappa=2$, $N=1024$, $\tau=0.1$, $K_{u}=2\times 10^{-5}$, $K_{v}=K_{u}/2$, and $F=0.03$ and vary $\lambda$ in a range in which the standard diffusion model is known to exhibit interesting dynamics [38]. For a ratio of diffusion coefficients $K_{u}/K_{v}>1$, the model is known to generate different mechanisms of pattern formation depending on the values of the feed rate $F$ and decay rate $\lambda$.

Figs. 3-5 show the evolutions of the numerical solution $v$ and summarize the effects of the super-diffusion for the fractional Gray–Scott model. The domain of interest is taken to be $(0,1)^{2}$. All the plots shown in Figs. 3-5 are snapshots from the numerical solutions $v$ in the above domain. The speeds of pattern formation are different for different $\alpha$ due to the influence of fractional powers.

For $\lambda=0.063$ (Fig. 3), the Gray–Scott model exhibits patterns of mitosis under conditions of the standard diffusion ($\alpha=2$). When $\alpha=1.7$, the replication pattern has completely changed as the fractional order of the model is decreased. As is shown for $\alpha=1.5$, the patterns present different behavior and exert dynamical states where solitons and filaments may coexist. We observe that the structure of patterns and the size of the spots are different. In Fig. 4 with $\lambda=0.061$, the original model produces a wavefront propagation partially driven by curvature. When the order $\alpha$ of the fractional Laplacian operator decreases, we can find that the area generated by diffusion propagation becomes larger. The final field is made up of much thinner filaments and transforms from point shape to line shape. For $\lambda=0.055$ (Fig. 5), the model with standard diffusion $(\alpha=2)$ is well known to organize in a steady state field of negative solitons. Moreover, the reduction of the fractional order $(\alpha=1.7)$ generates a decrease in the velocity of propagation of the initial perturbation, and affects the size of patterns with smaller spots. For smaller values of the fractional power $(\alpha=1.5)$, we observe a new process of nucleation of structures in the centre and in the boundaries of the domain. The process propagates outward until the entire area reaches the final steady state.

Example 5.3** (Fractional FitzHugh–Nagumo model).**

Here $\Omega=(0,2.5)\times(0,2.5)$, $u$ is the excitation variable which can be identified with the transmembrane potential, $v$ is the recovery variable, and $K_{u}$ is a small positive diffusion coefficient. In order to illustrate the effect of varying $\alpha$ in the model. The initial conditions are taken as

[TABLE]

[TABLE]

We choose the parameters $\kappa=2$, $\mu=0.1$, $\varepsilon=0.01$, $\beta=0.5$, $\gamma=1$, $\delta=0$, which are known to generate stable patterns in the system in the form of spiral waves. Stable rotating solutions at different time are presented in Fig. 6 to illustrate the effect of fractional order in the FitzHugh–Nagumo model with $K_{u}=1\times 10^{-4}$. We can observe that the width of the excitation wavefront is markedly reduced for decreasing $\alpha$, so is the wavelength of the system, with the domain being able to accommodate a larger number of wavefronts for smaller $\alpha$. By the comparison between Fig. 6 and Fig. 7 $(K_{u}=1\times 10^{-5})$, it is obvious that the role of reducing the fractional power $\alpha$ is not equivalent to the influence of a decreased diffusion coefficient $K_{u}$ in the pure diffusion case. For approximately the same width of the excitation wavefront, the wavelength of the model is larger in the fractional diffusion case, due to the long-tailed mechanisms of the fractional Laplacian operator. These results are consistent with those of Engler [48] for the fractional reaction-diffusion Fisher equation, showing distinct effects of fractional diffusion to those of reduced conductivity for a family of traveling wave solutions. They also illustrate the use of fractional diffusion as a modeling tool to characterize intermediate dynamic states not solely described by pure diffusion mechanisms.

6 Conclusions

Fractional differential equations have been proved to be useful tools for modeling diffusive processes associated with anomalous diffusion or spatial heterogeneity. In this paper, we study a second-order stabilized semi-implicit time-stepping Fourier spectral method for the reaction-diffusion systems of equations with space described by the fractional Laplacian in bounded rectangular domains. An optimal error estimate of the numerical scheme is obtained without imposing the CFL condition by adopting the temporal-spatial error splitting argument. We also analyze the linear stability of the stabilized semi-implicit method and obtain a practical criteria to choose the time step size to guarantee the stability of the semi-implicit method in real applications. Our approach is illustrated by solving several problems of practical interest, including the fractional Allen-Cahn, FitzHugh-Nagumo and Gray-Scott models show that such systems can have different dynamics to pure diffusion, and drastically different evolving patterns compared to integer-order models.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grants Nos. 11102102, 11672163), the Natural Science Foundation of Shandong Province (Grant ZR2017MA030) and the MURI/ARO on “Fractional PDEs for Conservation Laws and Beyond: Theory, Numerics and Applications (W911NF-15-1-0562)”.

Appendix A. Proof of Theorem 2.1

In the following, we denote $C,C_{q},\widetilde{C}_{q}$, and $\widehat{C}_{q}$ as positive constants that are independent of the time step size $\tau$ and the degree of approximation space $N$, where $q=1,2,\cdots$.

We first introduce some Lemmas that will be used in the numerical analysis later.

Lemma A.1** ([5]).**

Let $s,r\geq 0$, then for any $u,v\in H_{per}^{s+r}(\Omega)$, it holds that

[TABLE]

Lemma A.2** ([49, 50]).**

Suppose that $u\in H_{per}^{r}(\Omega)$. Then the following estimate holds for all $0\leq s\leq r$,

[TABLE]

Lemma A.3** ([51]).**

Let $\Omega$ be a domain in $\mathbb{R}^{n}$ satisfying the cone condition. If $mp>n$, let $p\leq q\leq\infty$; if $mp=n$, let $p\leq q<\infty$; if $mp<n$, let $p\leq q\leq q^{*}=np/(n-mp)$. Then there exists a constant $C_{2}$ depending on $m$, $n$, $p$, $q$, and the dimensions of $\Omega$ such that for all $u\in W^{m,p}(\Omega)$,

[TABLE]

where $\varepsilon=(n/mp)-(n/mq)$.

Lemma A.4** ([5]).**

If $s\leq r$ and $u\in H_{per}^{r}(\Omega)$, then

[TABLE]

Lemma A.5** ([5]).**

If $u,v\in H_{per}^{2s}(\Omega)$ for all $s\geq 0$, then

[TABLE]

where $C_{s}=\max{\{1,2^{2s-1}\}}$.

Lemma A.6** ([22]).**

For any $u\in X_{N}$, the following inverse inequality holds,

[TABLE]

We assume that the solution $u$ of (1) satisfies

[TABLE]

where $I=(0,T]$, $M$ is a constant independent of $n$, $N$, and $\tau$. The above regularity assumption is used in the convergence analysis.

We introduce a time-discrete system,

[TABLE]

Thus, we split the error $u^{n}-U_{N}^{n}$ into two parts, that is, $u^{n}-U_{N}^{n}=(u^{n}-U^{n})+(U^{n}-U_{N}^{n})$. The boundedness of $u^{n}-U_{N}^{n}$ can be obtained by estimating $u^{n}-U^{n}$ and $U^{n}-U_{N}^{n}$ separately.

Denote by $e^{n}=u^{n}-U^{n}$. Combining (8) and (36) yields

[TABLE]

where

[TABLE]

Define

[TABLE]

According to (35) and (39), we can find that $M_{1}$ is a positive constant independent of $n$, $\tau$, and $N$.

We have the following theorem for the time-discrete system (36).

Theorem A.1**.**

Suppose that (1) has a unique solution $u(\cdot,t)\in H^{\alpha}_{per}(\Omega)$ satisfying (35). Then, the time-discrete system (36) has a unique solution $U^{n}$. If $U^{n}\in H^{\alpha}_{per}(\Omega)$, then there exists a positive constant $\tau_{1}^{*}>0$, when $\tau\leq\tau_{1}^{*}$, it holds,

[TABLE]

[TABLE]

Proof.

When $n=1$, $e^{1}=0$. The proof is obvious.

Assume that (40) holds for $2\leq n\leq k-1$. By Lemma A.3, we have

[TABLE]

where $\tau\leq C_{2}^{-1}$.

Next, we prove (40) for $n\leq k$. Multiplying (37) by $2\tau(-\Delta)^{\alpha/2}D_{1}e^{n}$ and integrating the resulting system over $\Omega$, we obtain

[TABLE]

where Lemma A.1 and the Cauchy–Schwartz inequality is used. Summing $n$ in (43) from $2$ to $\hat{k}\,(2\leq\hat{k}\leq k)$ and using $e^{0}=e^{1}=0$ gives

[TABLE]

Next, we estimate $\|(-\Delta)^{\alpha/4}F_{1}^{n-1}\|^{2}$. From Lemma A.5, (42), and the mean value theorem, we derive

[TABLE]

where $\xi^{l}\in\left[\min\{u^{l},U^{l}\},\max\{u^{l},U^{l}\}\right]$, $|\xi^{l}|\leq\max{\{M,M_{1}\}}$, $l=n-2,n-1$.

Combining (44) and (LABEL:e3.30) leads to

[TABLE]

By the Gronwall inequality, one has $|e^{\hat{k}}|_{\alpha}\leq C_{6}\tau^{2}$, which implies

[TABLE]

when $\tau\leq(2C_{6})^{-1}$.

If $\tau\leq C_{2}^{-1}$, then

[TABLE]

Taking $\tau_{1}^{*}=\min\{C_{2}^{-1},(2C_{6})^{-1}\}$, we complete the proof. ∎

From Lemma A.2 and Lemma A.3 ($m=\alpha,p=n=2,q=\infty$), we have $\|P_{N}v\|_{\infty}\leq\widetilde{C}_{1}\|v\|_{\alpha}$ for any $v\in H^{\alpha}_{per}(\Omega)$. Therefore, we can obtain the boundedness of $\|P_{N}U^{n}\|_{\infty}$ for $n=1,2,\ldots,K$. Define the following constant

[TABLE]

Let $\eta_{N}^{n}=P_{N}U^{n}-U^{n}_{N}$. From (9) and (36), we obtain

[TABLE]

where

[TABLE]

The following Theorem gives a bound of $U_{N}^{n}$.

Theorem A.2**.**

Suppose that (36) has a unique solution $U^{n}\in H^{\alpha}_{per}(\Omega)$ and $U^{n}_{N}$ is the solution of (9), $0\leq n\leq K$. Then there exists two positive constants $\widetilde{C}_{4}$ and $N_{1}^{*}$, when $N\geq N_{1}^{*}$, it holds

[TABLE]

Proof.

For $n=0,1$, we have $\eta_{N}^{0}=\eta_{N}^{1}=0$. The proof is obvious.

Suppose that (49) holds for $2\leq n\leq k-1$. By Lemma A.6, we have

[TABLE]

where $N\geq(C_{3}\widetilde{C}_{4})^{1/(\alpha-1)}$.

Now, we prove (49) for $n\leq k$. Letting $v=2\tau D_{1}\eta^{n}_{N}$ in (47), it can be seen that

[TABLE]

Similar to (LABEL:e3.30), we obtain

[TABLE]

According to (52) and (LABEL:eqq3.58), one has

[TABLE]

Summing $n$ from $2$ to $\hat{k}$ $(2\leq\hat{k}\leq k)$, it arrives at

[TABLE]

where $\eta_{N}^{1}=0$ has been used. Combining (LABEL:eq3.55) and $\|\eta^{\hat{k}}_{N}\|^{2}=\tau^{2}\|\sum_{n=2}^{\hat{k}}D_{1}\eta^{n}_{N}\|^{2}\leq T\tau\sum_{n=1}^{\hat{k}}$ $\|D_{1}\eta^{n}_{N}\|^{2}$ presents

[TABLE]

Applying the Gronwall inequality,

[TABLE]

It is obvious that

[TABLE]

provided that $N\geq(C_{3}\widetilde{C}_{4})^{1/(\alpha-1)}$. Taking $N_{1}^{*}=(C_{3}\widetilde{C}_{4})^{1/(\alpha-1)}$, this proof is completed. ∎

Next, we prove the error estimate of of fully discrete scheme (9). Denote

[TABLE]

We can obtain the error equation of (9) as

[TABLE]

where

[TABLE]

We are now in a position to prove Theorem 2.1.

Proof.

It is obvious that (11) holds for $n=0,1$. Taking $v=2\tau D_{1}\theta_{N}^{n}$ in (60) leads to

[TABLE]

From Theorem A.2, we derive that $U_{N}^{n}$ is bounded, which indicates

[TABLE]

The following bound is easily obtained

[TABLE]

Combining (62), (LABEL:eq3.58), and (64) yields

[TABLE]

Summing $n$ from $2$ to $\hat{k}$ $(2\leq\hat{k}\leq k)$ presents

[TABLE]

where $\theta_{N}^{0}=\theta_{N}^{1}=0$ was used. Combining (LABEL:e3.55) and $\|\theta_{N}^{\hat{k}}\|^{2}=\tau^{2}\|\sum_{n=2}^{\hat{k}}D_{1}\theta^{n}_{N}\|^{2}\leq T\tau\sum_{n=1}^{\hat{k}}$ $\|D_{1}\theta^{n}_{N}\|^{2}$ gives

[TABLE]

By the Gronwall inequality, we obtain

[TABLE]

The above inequality and Lemma A.2 imply

[TABLE]

which ends the proof. ∎