An effective operator splitting method based on spectral deferred correction for the fractional Gray–Scott model

Highlights

• We propose an effective operator splitting method to simulate the fractional Gray–Scott (GS) model.

• The second-order operator splitting scheme, which is to split the original problem into linear and nonlinear parts.

• The linear subproblem is numerically solved using the Fourier spectral method.

• The nonlinear one is solved via the Crank–Nicolson formula and Rubin-Graves linearization technique.

• The stability and convergence of this method are analyzed.

Abstract

This paper presents a method by combining the semi-implicit spectral deferred correction (SDC) method with the operator splitting scheme to simulate the fractional Gray-Scott (GS) model. We start with the second-order operator splitting scheme, which is to split the original problem into linear and nonlinear parts. The linear subproblem is numerically solved using the Fourier spectral method, which is based on the exact solution and thus has no stability restriction on the time-step size. The nonlinear subproblem is solved via the Crank-Nicolson formula and Rubin-Graves linearization technique, which can be solved effectively. The stability and convergence of this method are analyzed in $L^2$-norm. Moreover, the scheme also takes advantage of the semi-implicit SDC method to improve the temporal accuracy. Numerical results are given to illustrate that the proposed method is a practical, accurate and efficient simulation tool for solving fractional GS problems.

Introduction

The Gray–Scott (GS) model of chemical kinetics which describes the spatiotemporal changes of the concentrations in the reactor is given by [1]

$$\left\{\begin{array}{c}{u}_t={\mu }_u\Delta u-u{v}^2+F(1-u),\hfill \\ v_t={\mu }_v\Delta v+u{v}^2-(F+\kappa )v,\hfill \end{array}\right.$$

where $u$ and $v$ are the concentrations, ${\mu }_u>0$ and ${\mu }_v>0$ are the diffusion rates in the process. $F\ge 0$ is the feed rate and $\kappa \ge 0$ is the decay rate of the second reaction. The theoretical property of this model has been investigated in [2].

The GS model is highly non-linear due to the presence of the coupling terms other than the time derivative and the diffusion terms. Many authors have proposed different numerical methods for solving the GS model, such as [1], [3], [4], [5], [6] Note that in these works on pattern dynamics of the GS model are all concerned with standard diffusion. However, recent studies show that non-homogeneities of the medium may fundamentally alter the laws of Markov diffusion, leading to long range fluxes, and non-Gaussian, heavy tailed profiles, and these motions may no longer obey the classical Fick’s law. This phenomenon is called anomalous diffusion. Fractional diffusion equations provide an adequate and accurate description of transport processes that exhibit anomalous diffusion, which cannot be modeled properly by integer-order diffusion equations [7], [8], [9], [10], [11]. In this context, the fractional GS model  [12], [13], [14] has attracted much attention. A clear physical interpretation of this model can be found in [13] and the references cited therein.

In this work, we consider the following fractional GS model [12]

$$\left\{\begin{array}{c}{u}_t=-{\mu }_u{(-\Delta )}^{\alpha /2}u-u{v}^2+F(1-u),(\mathbf{x},t)\in \Omega \times (0,T],\hfill \\ v_t=-{\mu }_v{(-\Delta )}^{\alpha /2}v+u{v}^2-(F+\kappa )v,(\mathbf{x},t)\in \Omega \times (0,T]\hfill \end{array}\right.$$

subject to the periodic boundary conditions and the initial conditions

$$u(\mathbf{x},0)=u_0\left(\mathbf{x}\right),v(\mathbf{x},0)=v_0\left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

Here, $\Omega ={[a,b]}^d$ (d = 1, 2, 3). The operator ${(-\Delta )}^{\alpha /2}$ with $\alpha \in (1,2]$ is the fractional Laplacian operator, which can be defined in various ways, such as the Riesz definition, the spectral definition and the directional definition. A review of these existing definitions can be seen in [15], [16]. In this paper, we utilize the following commonly used definition in spectral decomposition form [17], [18], [19]. For example, in the 2D case

$${(-\Delta )}^{\alpha /2}u=\sum _{p=-\infty }^{\infty }\sum _{q=-\infty }^{\infty }{\lambda }_{pq}^{\alpha }{\hat{u}}_{pq}{\phi }_{pq},$$

where

$${\lambda }_{pq}^{\alpha }={\left({\left(\frac{2p\pi }{b-a}\right)}^2+{\left(\frac{2q\pi }{b-a}\right)}^2\right)}^{\frac{\alpha }{2}},{\phi }_{pq}=exp\left(\text{i}\frac{2p\pi (x-a)}{b-a}+\text{i}\frac{2q\pi (y-a)}{b-a}\right),{\hat{u}}_{pq}=\frac{1}{{(b-a)}^2}{\int }_{\Omega }u{\bar{\phi }}_{pq}d\mathbf{x}.$$

Analytical solutions for the fractional GS model are not always feasible and thus discrete numerical techniques are desired. Bueno-Orovio et al. [9] used the Fourier spectral method to solve several types of fractional reaction–diffusion equations. The method in [9] required iteration and no theoretical analysis was provided. Wang et al. [13] proposed a second-order finite difference scheme for the fractional GS model and provided the stability analysis for the semi-discrete scheme in time. Liu et al. [14] introduced a temporal two-mesh finite element method and derived the stability and error estimates for the fully discrete scheme. Zhang et al. [12] proposed a stabilized semi-implicit Fourier spectral method for nonlinear space-fractional reaction–diffusion equations. The algorithm was successfully applied to solve the fractional GS model. Lee et al. [10] presented an efficient operator splitting method for solving the model (1.2). However, the method in [10] requires the Newton-type iteration and no theoretical analysis was provided. In addition, all the above schemes are not higher than second-order accuracy in time.

Our goal of this paper is to propose a numerical algorithm with high-order accuracy in both time and space for the fractional GS model (1.2). We will mainly follow the operator splitting method [20], [21] and deferred correction method [22], [23] to construct our time stepping scheme, and the space variables are discretized by using the Fourier-spectral method [24], [25], [26], [27], [28]. Operator splitting is a powerful method for the numerical investigation of complex time-dependent models and has been used to solve many complicated problems, such as [29], [30], [31], [32], [33], [34]. However, for the operator splitting method to be applied to nonlinear partial differential equations, the theoretical analysis turns out to be very challenging due to its multi-stage nature. The readers are referred to the works [35], [36], [37], [38], [39], [40] for more detailed information on the theoretical analysis of operator splitting method. In this paper, we obtain the desired convergence estimate of the operator splitting method for the fractional GS model. This study provides a general framework for the convergence estimate of the operator splitting method. The main contribution consists of four aspects:

$•$ A new linearizing algorithm is proposed based on the operator splitting method for the fractional GS model.

$•$ The stability and convergence of this method in $L^2$-norm are studied.

$•$ A semi-implicit SDC method is further used to improve time accuracy.

$•$ The results from three numerical experiments are presented to compare the dynamics of different $\alpha$.

The paper is organized as follows. In Section 2, we present some notations used for discretization and discuss some preliminary results concerning the fractional Laplacian operator. In Section 3, a new linearizing second-order operator splitting spectral scheme is proposed for the fractional GS model. In Section 4, we give the stability and convergence of the proposed method. In Section 5, a high-order semi-implicit SDC method combined with the operator splitting spectral scheme is introduced. Numerical experiments are presented to test the performance of the proposed numerical algorithm in Section 6, and some concluding comments are given in the final section.