Splitting spectral element method for fractional reaction-diffusion equations

Abstract

In this paper, we propose a second-order operator splitting spectral element method for solving fractional reaction-diffusion equations. In order to achieve a fast second-order scheme in time, we decompose the original equation into linear and nonlinear sub-equations, and combine a quarter-time nonlinear solver and a half-time linear solver followed by final quarter-time nonlinear solver. The spatial discretization is eigen-decomposition based on spectral element method. Since this method gives a full diagonal representation of the fractional operator and gets an exponential convergence in space. We have an accurate and efficient approach for solving spacial fractional reaction-diffusion equations. Some numerical experiments are carried out to demonstrate the accuracy and efficiency of this method. Finally, we apply the proposed method to investigate the effect of the fractional order in the fractional reaction-diffusion equations.

Keywords

Fractional reaction-diffusion equations operator splitting spectral element method eigen-decomposition exponential convergence

Introduction

In recent years, the fractional partial differential equations are wining more and more scientific applications across a variety of fields including control theory,¹ biology,^2,3 electrochemical processes,⁴ finance,⁵ etc. There have been many numerical methods for discretizing fractional operators and resolving fractional differential equations^6–11 etc. Splitting methods have been fruitfully used to solve large systems of partial differential equations.^12–14 Usually, it is very difficult to find the exact (analytical) solution of a given problem in practice. Thus, we can use numerical methods to obtain an approximate solution of the equations. There are many major benefits of operator splitting, including dimension reduction, problem simplification, preservation of any order accuracy in time, and computational speed-up for some complex problems (cf.^15–17 and the references therein). In this paper, we try to extend the splitting method for spacial fractional reaction-diffusion equations as follows

\frac{\partial u}{\partial t} = μ L^{α} (u) + N (u), (x, t) \in Ω \times (0, T)

(1)

u |_{\partial Ω} = b (x)

(2)

with

0 < α < 1

, where Ω is a bounded domain in

R^{2}

, μ is viscosity coefficient, the fractional operator is defined as

L^{α} = - {(- Δ)}^{α}

and the nonlinear term

N (u) = f (u)

(e.g. in the Fisher-Kolmogorov equation¹⁸

f = r u (1 - \frac{u}{K})

, here r and K > 0 are variable coefficients depending on the space). The initial condition is given as

u_{0} (x) = g (x)

The fractional Laplacian ${(- Δ)}^{α}$ in this paper is defined as follows

{(- Δ)}^{α} u (x, t) = \sum_{i = 1}^{\infty} λ_{i}^{α} u_{i} (t) ϕ_{i} (x)

(3)

where

{(λ_{i}, ϕ_{i})}_{i = 1}^{\infty}

are the eigenpairs of the standard Laplacian

- Δ

with certainty boundary conditions and u has the expansion

u (x, t) = \sum_{i = 1}^{\infty} u_{i} (t) ϕ_{i} (x)

. Alternatively, the fractional operator

{(- Δ)}^{α}

can be defined¹⁹ through the Fourier transform

\overset{}{{(- Δ)}^{α}} u (ω) = ω^{2 α} \hat{u} (ω)

(4)

where

\hat{u}

is the Fourier transform of u.²⁰

The remainder of this paper is organized as follows. In next section, we describe the numerical method for the fractional reaction-diffusion equations, and provide the error estimation for the splitting scheme. Several numerical examples are given in ‘Numerical results’ section. We will apply our method to solve the fractional Allen-Cahn, Fisher-Kolmogorov, Gray-Scott and FitzHugh–Nagumo equations.^2,11,21 The numerical results show that our method is both efficient and highly accurate for solving fractional reaction-diffusion equations. We provide a short summary in the final section.

Numerical method

We introduce a second-order splitting scheme for solving the above reaction-diffusion equation (1) and derive the error estimation for the scheme in this section. Let L be the number of the time step to integrate up to final time T, then the time step $Δ t = T / L, t_{n} = n Δ t, t_{n + \frac{1}{4}} = t_{n} + \frac{Δ t}{4}$ and $t_{n + \frac{3}{4}} = t_{n} + \frac{3 Δ t}{4}$ for $n = 0, 1, \dots, L - 1$ . We denote the time levels by superscripts $u^{n + 1} = u (t_{n + 1})$ and set initial condition $u^{0} = u_{0} (x)$ . We look for solution of $u^{n + 1}, n = 0, 1, \dots, L - 1$

{\begin{array}{l} \frac{\partial u^{(1)}}{\partial t} = 2 N (u^{(1)}), t \in (t_{n}, t_{n + \frac{1}{4}}) \\ \frac{\partial u^{(2)}}{\partial t} = 2 μ L^{α} (u^{(2)}), t \in (t_{n + \frac{1}{4}}, t_{n + \frac{3}{4}}) \\ \frac{\partial u}{\partial t} = 2 N (u), t \in (t_{n + \frac{3}{4}}, t_{n + 1}) \end{array}

(5)

with the initial conditions in each sub-step

u^{(1)} (t_{n}) = u^{n}, u^{(2)} (t_{n + \frac{1}{4}}) = u^{(1)} (t_{n + \frac{1}{4}})

and

u (t_{n + \frac{3}{4}}) = u^{(2)} (t_{n + \frac{3}{4}})

respectively.

Remark 1. In particular case we can find the analytic solutions of the initial value problem with a given nonlinear term f(u), since the ordinary differential equations in the first sub-step and the third sub-step of the above scheme (5) are solvable. Alternatively, we solve a diffusion equation for the second sub-step with an efficiency way since the fractional operator is defined as (3).

Next, we will derive the error estimation for the above scheme (5). Here we have second-order convergence in time $O (Δ t^{2})$ which is the same as the splitting method used for standard reaction-diffusion equation (i.e. α = 1) in Harwood.²²

Theorem 1. Given a semi-linear parabolic PDE, which can be written as $\frac{\partial u}{\partial t} = μ L^{α} (u) + N (u),$ where $L^{α}$ and $N$ are linear and nonlinear operator upon u, respectively, the solutions to the set of split equations

{\begin{array}{l} \frac{\partial u^{(1)}}{\partial t} = 2 N (u^{(1)}) \\ \frac{\partial u^{(2)}}{\partial t} = 2 μ L^{α} (u^{(2)}) \\ \frac{\partial u}{\partial t} = 2 N (u) \end{array}

(6)

provide an

O (Δ t^{2})

splitting accuracy when solved sequentially over the set of subintervals

{(t_{n}, t_{n + \frac{1}{4}}], (t_{n + \frac{1}{4}}, t_{n + \frac{3}{4}}], (t_{n + \frac{3}{4}}, t_{n + 1}]}

Proof. Using Taylor expansion at time step $t_{n + 1}$ , the solution to the full problem (1)

\begin{array}{l} u_{T}^{n + 1} = u^{n} + Δ t (μ L^{α} (u^{n}) + N (u^{n})) \\ + \frac{Δ t^{2}}{2} (μ^{2} L^{2 α} (u^{n}) + μ L^{α} (N (u^{n})) \\ + \frac{\partial N (u^{n})}{\partial u^{n}} μ L^{α} (u^{n}) + \frac{\partial N (u^{n})}{\partial u^{n}} N (u^{n})) \\ + O (Δ t^{3}) \end{array}

(7)

Combining the splitting equations (19) we have the following approximation, here we denote $u^{n + \frac{3}{4}} = u^{(2)} (t_{n + \frac{3}{4}}), u^{n + \frac{1}{4}} = u^{(1)} (t_{n + \frac{1}{4}})$ and I represents the identity operator for convenience. By using Taylor expansion

\begin{array}{l} u^{n + 1} = u^{n + \frac{3}{4}} + \frac{Δ t}{4} \frac{\partial u^{n + \frac{3}{4}}}{\partial t} + \frac{Δ t^{2}}{32} \frac{\partial^{2} u^{n + \frac{3}{4}}}{\partial t^{2}} + O (Δ t^{3}) \\ = u^{n + \frac{3}{4}} + \frac{Δ t}{2} N (u^{n + \frac{3}{4}}) + \frac{Δ t^{2}}{16} \frac{\partial N (u^{n + \frac{3}{4}})}{\partial t} + O (Δ t^{3}) \\ = u^{n + \frac{1}{4}} + μ Δ t L^{α} (u^{n + \frac{1}{4}}) + \frac{Δ t^{2} μ^{2}}{2} L^{2 α} (u^{n + \frac{1}{4}}) \\ + \frac{Δ t}{2} N (u^{n + \frac{1}{4}} + μ Δ t L^{α} (u^{n + \frac{1}{4}}) + O (Δ t^{2})) \\ + \frac{Δ t^{2}}{16} (\frac{\partial N (u^{n + \frac{1}{4}})}{\partial t} + O (Δ t)) + O (Δ t^{3}) \\ = u^{n + \frac{1}{4}} + μ Δ t L^{α} (u^{n + \frac{1}{4}}) + \frac{Δ t^{2} μ^{2}}{2} L^{2 α} (u^{n + \frac{1}{4}}) \\ + \frac{Δ t}{2} (N (u^{n + \frac{1}{4}}) + μ Δ t \frac{\partial N (u^{n + \frac{1}{4}})}{\partial u^{n + \frac{1}{4}}} L^{α} (u^{n + \frac{1}{4}}) + O (Δ t^{2})) \\ + \frac{Δ t^{2}}{8} \frac{\partial N (u^{n + \frac{1}{4}})}{\partial u^{n + \frac{1}{4}}} N (u^{n + \frac{1}{4}}) + O (Δ t^{3}) \end{array}

\begin{array}{l} = (I + μ Δ t L^{α} + \frac{Δ t^{2} μ^{2}}{2} L^{2 α}) (u^{n} + \frac{Δ t}{4} \frac{\partial u^{n}}{\partial t} + \frac{Δ t^{2}}{32} \frac{\partial^{2} u^{n}}{\partial t^{2}} + O (Δ t^{3})) \\ + \frac{Δ t}{2} (N (u^{n} + \frac{Δ t}{4} \frac{\partial u^{n}}{\partial t} + O (Δ t^{2})) \\ + Δ t (\frac{\partial N (u^{n})}{\partial u^{n}} + O (Δ t)) (μ L^{α} (u^{n}) + O (Δ t))) \\ + \frac{Δ t^{2}}{8} (\frac{\partial N (u^{n})}{\partial u^{n}} + O (Δ t)) (N (u^{n}) + O (Δ t)) + O (Δ t^{3}) \\ = u^{n} + Δ t (μ L^{α} (u^{n}) + \frac{1}{2} N (u^{n})) \\ + \frac{Δ t^{2}}{2} (μ^{2} L^{2 α} (u^{n}) + μ L^{α} (N (u^{n})) + \frac{1}{4} \frac{\partial N (u^{n})}{\partial u^{n}} N (u^{n})) \\ + \frac{Δ t}{2} N (u^{n}) + Δ t^{2} (\frac{1}{4} \frac{\partial N (u^{n})}{\partial u^{n}} N (u^{n}) + \frac{1}{2} \frac{\partial N (u^{n})}{\partial u^{n}} μ L^{α} (u^{n})) \\ + \frac{Δ t^{2}}{8} \frac{\partial N (u^{n})}{\partial u^{n}} N (u^{n}) + O (Δ t^{3}) \\ = u^{n} + Δ t (μ L^{α} (u^{n}) + N (u^{n})) \\ + \frac{Δ t^{2}}{2} (μ^{2} L^{2 α} (u^{n}) + μ L^{α} (N (u^{n})) \\ + \frac{\partial N (u^{n})}{\partial u^{n}} μ L^{α} (u^{n}) + \frac{\partial N (u^{n})}{\partial u^{n}} N (u^{n})) + O (Δ t^{3}) \end{array}

Where $u^{n + \frac{3}{4}} = u^{n + \frac{1}{4}} + μ Δ t L^{α} (u^{n + \frac{1}{4}}) + \frac{Δ t^{2} μ^{2}}{2} L^{2 α} (u^{n + \frac{1}{4}}) + O (Δ t^{3})$ is derived from the second sub-step of the scheme (6). The $O (Δ t^{2})$ splitting error can be obtained by computing the $O (Δ t^{3})$ solution difference from this ordering as

| u_{T}^{n + 1} - u^{n + 1} | = O (Δ t^{3})

Thus, the symmetric splitting scheme results in second-order accuracy in time. The proof is complete. □

In previous work, we have developed a numerical method for computing fractional Laplacians on complex-geometry domains.²³ By using the numerical eigenpairs $(λ_{i}, ϕ_{i})$ of the Laplace operator $- Δ$ from the SEM solution, we can approximate the fractional Laplace operator as

{(- Δ)}^{α} u \approx \sum_{i = 1}^{N_{d}} u_{i} λ_{i}^{α} ϕ_{i}

(8)

where

u_{i} = (ϕ_{i}, u), (\cdot, \cdot)

represents inner product and

N_{d}

denotes the degree of freedom (DoF) which depends on the element number El and the polynomial degree N of the problem. The variational statement for the second sub-step of the scheme (5) reads as: Find

u_{N_{d}} \in V_{N_{d}} (Ω), \forall v \in V_{N_{d}} (Ω)

so that

(\frac{\partial u_{N_{d}}}{\partial t}, v) = - 2 μ ({(- Δ_{N_{d}})}^{α} u_{N_{d}}, v)

(9)

where the time interval is

t \in (t_{n + \frac{1}{4}}, t_{n + \frac{3}{4}}]

and

V_{N_{d}} (Ω) = span {ϕ_{i}, i = 1, \dots, N_{d}}

and

u^{(2)} (t)

in equation (5) is replaced by u(t) for convenience. The numerical solution is expanded as

u_{N_{d}} = \sum_{i = 1}^{N_{d}} u_{i} (t) ϕ_{i}

and the fractional operator

{(- Δ_{N_{d}})}^{α}

is defined as follows

{(- Δ_{N_{d}})}^{α} u = \sum_{i = 1}^{N_{d}} λ_{i}^{α} u_{i} (t) ϕ_{i}

(10)

By the orthogonality of the eigenfunctions we obtain

\partial_{t} u_{i} (t) = - 2 μ u_{i} (t) λ_{i}^{α}

(11)

The solution of the above equations with initial condition $u_{i} (t_{n + \frac{1}{4}})$ is given as follows

u_{i} (t) = e^{- 2 μ λ_{i}^{α} (t - t_{n + \frac{1}{4}})} u_{i} (t_{n + \frac{1}{4}})

where

u_{i} (t_{n + \frac{1}{4}}) = \int_{Ω} u (x, t_{n + \frac{1}{4}}) ϕ_{i} d x

. Thus we obtain the analytic solution of (9)

u_{N_{d}} (x, t) = \sum_{i = 1}^{N_{d}} e^{- 2 μ λ_{i}^{α} (t - t_{n + \frac{1}{4}})} u_{i} (t_{n + \frac{1}{4}}) ϕ_{i}

(12)

Numerical results

Fractional Allen-Cahn equation (verification of convergence results)

In order to test the convergence rates of the splitting scheme (5) for reaction-diffusion equations, we consider the fractional Allen-Cahn equation

\frac{\partial u}{\partial t} = - μ {(- Δ)}^{α} u + u - u^{3}, (x, t) \in Ω \times (0, T), 0 < α \leq 1

(13)

with the Neumann boundary condition

\frac{\partial u (x, t)}{\partial n} = 0, x \in \partial Ω

(14)

where the bounded domain

Ω = {[- 1, 1]}^{2}

and μ is the viscosity. The initial condition is given as

u (x, 0) = \sin (π x) \cos (2 π y)

The corresponding nonlinear initial value problem for the Allen-Cahn equation is given as follows

\frac{\partial u}{\partial t} = 2 (u - u^{3}), s \in (t_{ξ}, t]

(15)

here the initial condition is given as

u (t_{ξ})

, here ξ is equal to n or

n + \frac{3}{4}

at the n-th time step. The equation (15) can be integrated to obtain^24,25

u (s) = \frac{u (t_{ξ})}{\sqrt{e^{- 4 (s - t_{ξ})} + u {(t_{ξ})}^{2} (1 - e^{- 4 (s - t_{ξ})})}}, s \in (t_{ξ}, t)

(16)

We set $μ = 0.01, T = 2$ , the element number El = 64 and the polynomial degree N = 20 in each element such that the spatial discretization errors are negligible compared to the time discretization errors. Since the analytical solution is not available in this test, we choose the result obtained from the finest time step $τ_{0} = 10^{- 3}$ as our reference solution and compute the L² errors between the reference and the numerical solutions from several time steps $Δ t$ , denote $e_{1} = | | u_{N_{d}, τ_{0}} - u_{N_{d}, Δ t} | |_{L^{2} (Ω)}$ . We also use the numerical solution $u_{N_{d}, Δ t}^{f}$ as another reference solution, where $u_{N_{d}, Δ t}^{f}$ is computed by using the method we developed in Song et al.,²³ and denote $e_{2} = | | u_{N_{d}, Δ t}^{f} - u_{N_{d}, Δ t} | |_{L^{2} (Ω)}$ . In Figure 1, we plot the L² errors e₁ and e₂ with different time steps $Δ t$ for several fractional orders α. It is clear from Figure 1 that the scheme is second-order accurate method in time for any α. Moreover, the splitting scheme (5) is more accuracy than the scheme in Song et al.²³ for $α < 1$ , since we solve the nonlinear equation (15) exactly in each time step. In Figure 2, we also show the spatial errors in semi-log scale for several α. The time step is $Δ = 10^{- 3}$ and final time T = 1 in this test. It is observed that this method has exponential convergence rate for all α. Since the analytical solutions are not available in this test either, we choose the result obtained from the finest mesh N = 25 in each element as the reference solution and compute the L² errors between the reference and the numerical solutions from several polynomial degree N, denote $e = | | u_{25} - u_{N} | |_{L^{2} (Ω)}$ .

Figure 1.

L² error as a function of $Δ t$ in log-log scale with different fractional order α of the Allen-Cahn equation (13).

Figure 2.

L² error as a function of the polynomial degree N in semi-log scale with different fractional order α of the Allen-Cahn equation (13).

Fisher-Kolmogorov model

In this subsection, we consider the Fisher-Kolmogorov equation

\frac{\partial u}{\partial t} = - μ {(- Δ)}^{α} u + r u (1 - \frac{u}{K}), (x, t) \in Ω \times (0, T)

(17)

with the Dirichlet boundary condition

u (x, t) = 0, x \in \partial Ω

(18)

where the bounded domain $Ω = [- 3, 10] \times [- 2, 2], μ = 0.02, r = 0.2$ , K(x, y) is a variable coefficient

K (x, y) = {\begin{array}{l} 10^{- 7}, x \in Ω_{b} = [1, 2] \times ([- 2, 0.2] \cup [0.4, 2]) \\ 1, otherwise \end{array}

(19)

In applications to biology, Ω _b represents a region where populations cannot grow, due to unfavorable environmental conditions. The geometry is a slitted barrier, through which the solution will eventually penetrate. K smoothly interpolate on the interface between Ω _b and outside region for our simulation. The initial condition is given as

u_{0} (x) = 0.8 \exp (- 25 (x^{2} + y^{2}))

We consider the corresponding nonlinear initial value problem for the Fisher-Kolmogorov model as follows

\frac{\partial u}{\partial t} = 2 r u (1 - \frac{u}{K}), s \in (t_{ξ}, t)

(20)

here the initial condition is given as $u (t_{ξ})$ . The equation (20) can be integrated to obtain

u (s) = \frac{K u (t_{ξ}) e^{2 r (s - t_{ξ})}}{K + u (t_{ξ}) (e^{2 r (s - t_{ξ})} - 1)}, s \in (t_{ξ}, t)

(21)

We set $Δ t = 0.0002, T = 90$ , the element number El = 117 and the polynomial degree N = 8 in each element to this simulation. Figure 3 shows the evolution of the solution for different fractional order $α = 1.0, 0.85$ . These results conform to the results presented in Baeumer et al.²⁶ Figure 3(a) shows the solution in the case of α = 1 (i.e. the classical reaction-diffusion equation) at t = 50, 70, 90. We can observe that the solution is very slow to penetrate the barrier. Next we change $α = 0.85$ to represent anomalous diffusion, and repeat the numerical simulation. Figure 3(b) shows the solution has penetrated significantly at time t = 70. Due to the strongly non-local character of the fractional operator, it is much easier for members of the population to cross (jump) over the barrier. Figure 4 shows the population masses versus time t for the cases $α = 1, 0.85$ in two different regions $Ω_{1} = [- 3, 1] \times [- 2, 2]$ and $Ω_{2} = [2, 10] \times [- 2, 2]$ . The population mass of the fractional case (i.e. $α = 0.85$ ) in $Ω_{2}$ grows significantly faster than the classical case (i.e. α = 1) after time t = 40.

Figure 3.

Evolution of the solution corresponding to Fisher-Kolmogorov model. u contours are shown at times t = 50, 70, 90 from top to bottom. (a) $α = 1.0$ , (b) $α = 0.85$ .

Figure 4.

Population masses of the Fisher-Kolmogorov model versus t: $M_{i} = \int_{Ω_{i}} u d x, i = 1, 2$ , here $Ω_{1} = [- 3, 1] \times [- 2, 2]$ and $Ω_{2} = [2, 10] \times [- 2, 2]$ .

Gray-Scott model

In this subsection, we consider the Gray-Scott model. The Gray-Scott equations describe an autocatalytic reaction-diffusion process between two chemical constituents with concentrations u and v

\begin{array}{l} \frac{\partial u}{\partial t} = - μ_{u} {(- Δ)}^{α} u - u v^{2} + c (1 - u) \\ \frac{\partial v}{\partial t} = - μ_{v} {(- Δ)}^{α} v + u v^{2} - (c + k) v \end{array}

(22)

with periodic boundary condition, the corresponding chemical reaction is as follows

\begin{array}{l} U + 2 V \to 3 V \\ V \to P \end{array}

(23)

where the bounded domain

Ω = {[- 1, 1]}^{2}, U, V

and P are chemical species, u and v represent their concentrations, μ_u and μ_v are the diffusion rates, k represents the rate of conversion of V to P and c represents the rate of the process that feeds U and drains

U, V

and P.

By splitting the linear terms of the reactive contribution into the linear system

\begin{array}{l} \frac{\partial u}{\partial t} = - μ_{u} {(- Δ)}^{α} u + c (1 - u) \\ \frac{\partial v}{\partial t} = - μ_{v} {(- Δ)}^{α} v + (c + k) v \end{array}

(24)

the corresponding nonlinear system for $t \in (t_{ξ}, t)$

\begin{array}{l} \frac{\partial u}{\partial t} = - 2 u v^{2} \\ \frac{\partial v}{\partial t} = 2 u v^{2} \end{array}

(25)

which has the conservation relation

u + v = u_{0} + v_{0}

, where u₀ and v₀ are initial conditions. This allows the nonlinear system to be uncoupled to obtain

v_{t} = (u_{0} + v_{0} - v) v^{2}

, which can be integrated to obtain

{\begin{array}{l} u = u_{0} + v_{0} - v \\ v = \frac{u_{0} + v_{0}}{1 + W_{0} (\frac{u_{0}}{v_{0}} \exp (\frac{u_{0}}{v_{0}} - 2 {(u_{0} + v_{0})}^{2} (s - t_{0})))}, \end{array} s \in (t_{0}, t)

(26)

where W₀ is the principal branch of the Lambert W function.^27,28 Writing equation (26) in terms of W₀ is only possible when

v_{0} \neq 0

. In this simulation, the initial condition

v_{0} = 0

in the some regions of the domain. The Lambert W function can not be used in this case. Alternatively, we use fix-point iteration method for solving the nonlinear equation as follows

\frac{\partial v}{\partial t} = 2 (u_{0} + v_{0} - v) v^{2}

(27)

from the conservation relation

u_{0} + v_{0} = u + v

We set $Δ t = 0.1, Ω = {[0, 1]}^{2}$ the element number El = 1600 and the polynomial degree N = 8 in each element to this simulation. The coefficients are given as $μ_{u} = 2 \times 10^{- 5}, μ_{v} = μ_{u} / 2$ , c = 0.03 and $k = k_{1}, k_{2}$ , where $k_{1} = 0.055, k_{2} = 0.063$ . For the ratio of diffusion coefficients $μ_{u} / μ_{v} > 1$ , the model is known to generate different mechanisms of pattern formation depending on the values of the feed c and depletion k rates. We select k in a range in which the standard diffusion model is known to exhibit interesting dynamics.³ The initial condition of the entire system is placed in the state $(u_{0}, v_{0}) = (1, 0)$ , and the center plate ${(x - 0.5)}^{2} + {(y - 0.5)}^{2} \leq {0.04}^{2}$ of the domain is perturbed to $(u_{0}, v_{0}) = (1 / 2, 1 / 4)$ . Figure 5 shows the evolutions of the numerical solution v and summarizes the effects of the super-diffusion for the fractional Gray-Scott model. Since the evolutions of the numerical solution u is similar to that the v, we will not elaborate it in this paper. In Figure 5(a) to (c) with $κ = 0.063$ , under the influence of the standard diffusion (α = 1), the Gray-Scott model exhibits patterns of mitosis.³ However, in the super-diffusion case $α = 0.85, 0.75$ , the way of the replication has a slight change and the dissociated spots become smaller. In Figure 5(d) to (f) with $κ = 0.055$ , the Gray-Scott model with the normal diffusion ( $α = 1.0$ ) produces a circular wave propagating outward and forms a structured pattern shown in picture (d). Moreover, the reduction of the fractional order $α = 0.85$ effects the size of patterns with smaller spots. For smaller values of the fractional order $α = 0.75$ , a new process of nucleation of structures in the center and in the boundaries of the domain is observed, that then grow outward until the entire domain reaches its final steady state configuration. These results show the same phenomenon as the results in literature.^2,29

Figure 5.

Evolution of the solution v in Gray-Scott equation (22): (a)-(c) v contours are shown at times t = 2000, 4600, 5200, 6200, 7200, 9200, 12000, 18000 from top to bottom for three fractional orders $α = 1.0, 0.85, 0.75$ with k = 0.063; (d)-(f) v contours are shown at times t = 2000, 4600, 5200, 6200, 7200, 9200, 12000, 18000 from top to bottom for three fractional orders $α = 1.0, 0.85, 0.75$ with k = 0.055. (a) $α = 1.0$ , (b) $α = 0.85$ , (c) $α = 0.75$ , (d) $α = 1.0$ , (e) $α = 0.85$ , (f) $α = 0.75$ .

FitzHugh–Nagumo model

In the last subsection, we consider the fractional FitzHugh-Nagumo model as follows

\frac{\partial u}{\partial t} = - K_{u} {(- Δ)}^{α} u + u (1 - u) (u - σ) - v

(28)

\frac{\partial v}{\partial t} = ϵ (β u - r v - δ)

(29)

with homogeneous Neumann boundary condition on the bounded domain

Ω = {[- 1.25, 1.25]}^{2}

, where u is a normalized trans-membrane potential and v is a dimensionless time-dependent recovery variable, K_u is the diffusion coefficient. We set the parameters of the model:

K_{u} = 10^{- 4}, σ = 0.1, ϵ = 0.01, β = 0.5, r = 1

and δ = 0, which is known to generate stable patterns in a system in the form of reentrant spiral waves. In order to illustrate the effect of varying fractional order α in the FitzHugh–Nagumo model, we take an initial condition as

\begin{array}{l} u (x, 0) = {\begin{array}{l} 1, x, y < 0 \\ 0, otherwise \end{array} \\ v (x, 0) = {\begin{array}{l} 0.1, y > 0, \forall x \\ 0, otherwise \end{array} \end{array}

(30)

The time step in the numerical experiments is $Δ t = 0.1$ .

Figure 6 presents the stable rotation solutions of the FitzHugh-Nagumo model at t = 2000 for $α = 1, 0.85, 0.75$ . We find that as expected, the wave travels more slowly as fractional order decreases. The width of the excitation wavefront is markedly reduced for decreasing α, as is the wavelength of the system, with the domain being able to accommodate a larger number of wavefronts for smaller fractional order α. These visualizations are consistent with results presented in Bueno-Orovio et al.² (Figure 3) for time up to t = 2000 for all α.

Figure 6.

Spiral waves in the FitzHugh-Nagumo model for varying α. (a) a = 1.0, (b) a = 0.85, (c) a = 0.75.

Summary

We develop a second-order operator splitting algorithm for solving reaction-diffusion equations with fractional operator. We analyze the discretization errors with benchmark solutions that show second-order convergence in time. We employ the spectral element method in space to solve the fractional diffusion equation. And we observe an exponential convergence in space. We investigate the effect of the fractional order for the different models, including the Allen-Cahn, Fisher-Kolmogorov, Gray-Scott and FitzHugh–Nagumo equations. Moreover, our method can be applied to the fractional reaction-diffusion models on complex-geometry domains, since the fractional Laplacian is defined as (8). We can use the spectral element method for solving the integer-order eigenvalue problem on complex-geometry domains. After we obtain the numerical eigenpairs of the eigenvalue problem, we can approach the fractional Laplace based on the eigenvalue decomposition on these domains. The approach is efficient as it does not require the solution of any linear systems since the eigenfunction decomposition leads naturally to a system of ODEs, which can be readily solved with standard methods.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the NSF of China (Grant No. 11901100), the Natural Science Foundation of Fujian Province of China (Grant No. 2020J02001) and Funds of Education Department of Fujian Province (Grant No. 510881/GXRC-20046).

ORCID iD

Fangying Song

References

A space-time spectral method for the time fractional diffusion optimal control problems. Adv Differ Equ 2015; 2015: 156.

Bueno-Orovio

Kay

Burrage

Fourier spectral methods for fractional-in-space reaction-diffusion equations. BIT Numer Math 2014; 54: 937–954.

Pearson

JE.

Complex patterns in a simple system. Science 1993; 261: 189–192.

Yuste

Acedo

Lindenberg

Reaction front in an a+B →C reaction-subdiffusion process. Phys Rev E 2004; 69: 36126.

Scalas

Gorenflo

Mainardi

Fractional calculus and continuous-time finance. Phys A – Stat Mech Appl 2000; 284: 376–384.

Lin

Finite difference/spectral approximations for the time-fractional diffusion equation. J Comput Phys 2007; 225: 1533–1552.

Zhang

Meerschaert

Baeumer

Particle tracking for time-fractional diffusion. Phys. Rev. E 2008; 78: 036705.

Wang

Tian

A fast and faithful collocation method with efficient matrix assembly for a two-dimensional nonlocal diffusion model. Comput Meth Appl Mech Eng 2014; 273: 19–36.

Liu

Zhuang

Anh

, et al. Stability and convergence of the difference methods for the space–time fractional advection–diffusion equation. Appl Math Comput 2007; 191: 12–20.

10.

Deng

Finite element method for the space and time fractional Fokker–Planck equation. SIAM J Numer Anal 2009; 47: 204–226.

11.

Zeng

Liu

, et al. Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation. SIAM J Numer Anal 2014; 52: 2599–2622.

12.

Ladics

Application of operator splitting to solve reaction-diffusion equations. Proc 9th Coll QTDE 2012; 9: 1–20.

13.

Feng

N-N

Zhou

G-R

Huang

W-P.

An efficient split-step time-domain beam-propagation method for modeling of optical waveguide devices. J Lightwave Technol 2005; 23: 2186.

14.

Chen

Shen

Spectral direction splitting schemes for the incompressible Navier-Stokes equations. E Asian J Appl Math 2011; 1: 215–234.

15.

Guermond

J-L

Minev

Salgado

Convergence analysis of a class of massively parallel direction splitting algorithms for the Navier-Stokes equations in simple domains. Math Comp 2012; 81: 1951–1977.

16.

Lee

HG.

A second-order operator splitting Fourier spectral method for fractional-in-space reaction–diffusion equations. J Comput Appl Math 2018; 333: 395–403.

17.

Peaceman

Rachford

Jr.

The numerical solution of parabolic and elliptic differential equations. J Soc Ind Appl Math 1955; 3: 28–41.

18.

Adomian

Fisher-Kolmogorov equation. Appl Math Lett 1995; 8: 51–52.

19.

Constantin

Behavior of solutions of 2D quasi-geostrophic equations. SIAM J Math Anal 1999; 30: 937–948.

20.

Taylor

ME.

Tools for PDE: pseudo-differential operators, paradifferential operators, and layer potentials. Volume 81. Cambridge: Cambridge Univ Press, 2000.

21.

Liu

Zhuang

Turner

, et al. A semi-alternating direction method for a 2-D fractional FitzHugh–Nagumo monodomain model on an approximate irregular domain. J Comput Phys 2015; 293: 252–263.

22.

Harwood

RC.

Operator splitting method and applications for semilinear parabolic partial differential equations. PhD thesis, Washington State University, USA, 2011.

23.

Song

Karniadakis

GE.

Computing fractional laplacians on complex-geometry domains: algorithms and simulations. SIAM J Sci Comput 2017; 39: A1320–A1344.

24.

Lee

Jeong

, et al. An unconditionally stable hybrid numerical method for solving the allen–cahn equation. Comput Math Appl 2010; 60: 1591–1606.

25.

Lee

Kim

A fast, robust, and accurate operator splitting method for phase-field simulations of crystal growth. J Crystal Growth 2011; 321: 176–182.

26.

Baeumer

Kovács

Meerschaert

MM.

Numerical solutions for fractional reaction–diffusion equations. Comput Math Appl 2008; 55: 2212–2226.

27.

Corless

Gonnet

Hare

, et al. On the lambert W function. Adv Comput Math 1996; 5: 329–359.

28.

Veberič

Lambert W function for applications in physics. Comput Phys Commun 2012; 183: 2622–2628.

29.

Wang

Song

Wang

, et al. Fractional Gray–Scott model: well-posedness, discretization, and simulations. Comput Meth Appl Mech Eng 2019; 347: 1030–1049.