Fractional Gray-Scott Model: Well-posedness, Discretization, and Simulations1footnote 11footnote 1This work was supported by the OSD/ARO/MURI on “Fractional PDEs for Conservation Laws and Beyond: Theory, Numerics and Applications (W911NF-15-1-0562)" and the National Science Foundation under Grant DMS-1620194. The first author was supported by the China Scholarship Council under 201706220157.

Fractional Gray-Scott Model: Well-posedness, Discretization, and Simulations111This work was supported by the OSD/ARO/MURI on “Fractional PDEs for Conservation Laws and Beyond: Theory, Numerics and Applications (W911NF-15-1-0562)" and the National Science Foundation under Grant DMS-1620194. The first author was supported by the China Scholarship Council under 201706220157.

Tingting Wang222School of Mathematics, Shandong University, Jinan 250100, Shandong, China. 333Division of Applied Mathematics, Brown University, Providence, RI 02912, USA (,,    Fangying Song333Division of Applied Mathematics, Brown University, Providence, RI 02912, USA (,,    Hong Wang444Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA (    George Em Karniadakis333Division of Applied Mathematics, Brown University, Providence, RI 02912, USA (,,

The Gray-Scott (GS) model represents the dynamics and steady state pattern formation in reaction-diffusion systems and has been extensively studied in the past. In this paper, we consider the effects of anomalous diffusion on pattern formation by introducing the fractional Laplacian into the GS model. First, we prove that the continuous solutions of the fractional GS model are unique. We then introduce the Crank-Nicolson (C-N) scheme for time discretization and weighted shifted Grünwald difference operator for spatial discretization. We perform stability analysis for the time semi-discrete numerical scheme, and furthermore, we analyze numerically the errors with benchmark solutions that show second-order convergence both in time and space. We also employ the spectral collocation method in space and C-N scheme in time to solve the GS model in order to verify the accuracy of our numerical solutions. We observe the formation of different patterns at different values of the fractional order, which are quite different than the patterns of the corresponding integer-order GS model, and quantify them by using the radial distribution function (RDF). Finally, we discover the scaling law for steady patterns of the RDFs in terms of the fractional order .

Key words. pattern formation, ADI algorithm, anomalous transport, finite difference, spectral collocation, radial distribution function

1 Introduction

In the past several decades, the formation of spatial and temporal patterns has become a very active area of research. There are many diverse patterns formed in physical, biological and chemical systems [12, 21, 27]. Among various systems, the reaction and diffusion systems attract much attention, since they create a variety of patterns that could be found in nature, for instance, spots, spot replication, stripes, and travelling waves (e.g. [3, 22]). A representative reaction and diffusion model is the GS model, which is a variant of the autocatalytic Selkov model of glycolysis [9, 24]. This model includes the following two reactions


where and represent the chemical species. The two reactions take place in an open flow reactor, in which is continuously supplied and the final product is removed [17]. The first reaction shows a process of autocatalysis and the second reaction describes the decay of into . In addition, note that there is a non-equilibrium constraint on by constantly feeding it and removing , which leads to a variety of unstable phenomena.

The classical (integer-order) GS model is expressed in the form


where and are the concentrations of the two chemical components. is the feed rate and is the decay rate of the second reaction; and are the diffusion coefficients.

The response of one dimensional GS model was studied by Vastano and co-works previously [27]. Doelman et al. [5, 6] investigated the asymptotic scaling of parameters and variables necessary for the analysis of the patterns in one dimension. Pearson [22] studied this system in two dimensions and presented very complicated spatio-temporal patterns. Pearson also did a thorough numerical study for this system. Many complex structures were observed in the numerical solutions. Hale studied the exact homoclinic and heteroclinic solutions of the GS model for autocatalysis [10]. McGough and Riley produced the bifurcation analysis to support the existing numerical evidence for patterns and derived the bifurcation results for nonuniform steady states [18].

Figure 1.1: The difference in pattern formation and the corresponding RDFs between normal diffusion and anomalous diffusion described by the GS model for .

Note that, in the aforementioned literature on pattern dynamics of the GS model, the models are all with standard diffusion, i.e. the diffusion operator is the normal Laplacian. In addition, different pattern formation was observed by changing the values of parameters and . Hence, here we aim to study the effects of the super-diffusion (with fractional Laplacian for ) on pattern formation of this model and perform some theoretical analysis for the GS model. Fig. LABEL:Fig._0 shows the different pattern formation and the corresponding radial distribution functions (RDFs) between normal diffusion and anomalous diffusion for the GS model. The standard diffusion systems correspond to Brownian motion, while the fractional diffusion systems correspond to Lévy process. The relationship between Lévy process and certain types of space fractional models has been established in [20]. In the paper by Wu et al. [15], the Turing instability and pattern formation of the Lengyel-Epstein model with super-diffusion were studied. In that paper, it was emphasized that more complex dynamics will appear under the super-diffusion. Bueno-Orovio and co-workers used the Fourier spectral method to solve several types of fractional reaction diffusion equations in their paper [2]. Lee [13] also introduced a second-order operator splitting Fourier spectral method to approximate the fractional-in-space reaction diffusion equations. All these papers described investigations of the fractional reaction-diffusion models. So in this paper, we will consider the fractional GS model with space fractional Laplacian (super-diffusion effects) in two-dimensions.

Since the exact solutions for the fractional GS model are not known, we will develop second-order numerical methods to solve the GS equations. The standard approach for solving space fractional diffusion equations is to use finite difference method, finite element method to discrete fractional derivatives and then use Euler formulation for the evolution of time [4, 11, 16, 19]. However, these methods require the solution of a linear system at each time step, which correspond to a large, dense matrix due to the nonlocal nature of the fractional operator. Roop et al. [7] analyzed a fully discrete finite element approximation to a time dependent fractional order diffusion equation, which contains a nonlocal quadratic nonlinearity. In the paper by Meerschaert [19, 25], the approach based on the C-N method combined with spatial extrapolation was used to derive temporally and spatially second-order accurate numerical estimates. Deng et al. [26] introduced the weighted shifted Grünwald difference operator to approximate the Riemann-Liouville fractional derivative and obtained a second-order accuracy in space. To overcome the bottleneck of expensive computations, Wang et al. introduced a fast numerical algorithm, which tackled the problems for two or three dimensions effectively [30, 31]. In addition to these classical numerical methods, spectral methods also have been used for space fractional equations [2, 13].

In this paper we shall use the finite difference method to derive a second-order numerical scheme for the fractional GS model. In this fractional model, the classical Laplacian operator is replaced by the fractional Laplacian operator with . We apply the weighted shifted Grünwald difference discretization method, which leads to well structured, relatively sparse and positive definite coefficient matrix [26]. For time direction, the C-N scheme is employed to obtain a temporally second-order estimate. Since this fractional GS model is a nonlinear system, we use the second-order implicit-explicit methods [1] to handle this problem, i.e. an implicit scheme is used for the linear terms and an explicit scheme is used for the nonlinear terms. Moreover, we carry out the linear stability analysis for the steady states of the GS model. We also derive the well-posedness of the fractional GS model. In addition, we provide the stability analysis for the time semi-discrete scheme. Several numerical experiments have been conducted to verify the accuracy in time and space of this numerical scheme. In the simulations of the fractional GS model, a small perturbation has been added to the initial states. We observe the formation of patterns under the condition of different parameter values. Moreover, the spectral collocation method is used in space discretization to simulate this fractional model. We compare the steady patterns obtained by the two different numerical methods and calculate the corresponding RDFs. Finally, we estimate the scaling law between the fractional orders and the distance between all spot pairs in the steady spot patterns.

This paper is organized as follows. In section 2, we introduce the fractional GS model and perform stability analysis for steady states. The well-posedness of this model is presented in section 3. In section 4, we propose a second-order accurate both in time and space numerical scheme for the discretization of the fractional GS model. We also provide the stability analysis for the time semi-discrete scheme. In section 5, we present numerical simulations of the fractional GS model, including convergence results and the scaling law for steady patterns of the RDFs. We conclude in section 6 with a summary.

2 The fractional GS model

The fractional GS model that describes an autocatalytic reaction-diffusion process between two chemical species with concentrations and is written as:


where . The diffusion coefficients satisfy . The parameters are positive constants representing feed rate and decay rate, respectively. In this paper we define the fractional Laplacian operator by Riesz fractional derivatives as follows


with and being the Riemann-Liouville fractional operators.

The fractional GS system is an activation-substrate depletion system [17]. The chemical specie grows auto-catalytically on the the specie i.e., the continuously fed substrate. For instance, the concentrations vary opposite to each other. In other words, the existence of in the second equation of (LABEL:s1:e1) will prompt the production of and reduce the concentration of . In practical experiments, the formation of patterns is affected by the parameter values and anomalous diffusion.

2.1 Steady states and linear stability analysis

In this subsection, we consider the spatially uniform steady states of the fractional GS model with :


This system has an unique trivial steady state for all the values of and . In addition, there also exists two steady states and when . Namely, we have


As a result, we can get the saddle-node bifurcation


In order to analyze the stability of these spatially uniform steady states, we introduce perturbations to the fractional system by and and obtain


By neglecting terms , we can get the normal mode solution


with amplitudes and wave number , .

We first investigate the eigenvalues of their characteristic equations to examine the stability of the steady states. By substituting into , the corresponding characteristic equation is derived as follows


Then we obtain the following dispersion relation



From , we get the eigenvalues corresponding to the trivial steady state,


which depend on the fractional order . These are extensive models of the integral case [18]. For , it is obvious that this trivial steady state is stable for all values of and . In terms of the other two alternative steady states, first we insert the point given by (LABEL:ss:e2) to (LABEL:ss:ee6). The computation shows that for . Since the steady states satisfy the equations in (LABEL:ss:e1), making use of , we find that and for . Therefore, the steady state is always unstable. However, the steady state may have stable and unstable structure.

We consider the case that the diffusion constants are negligibly small with ratio of order 1. In this case, when in (LABEL:ss:ee6) is purely imaginary, the system undergoes a Hopf bifurcation which closes to the lower branch of the saddle-node curve (see [17]). When , the eigenvalues of equation (LABEL:ss:ee6) are complex conjugate. Then the Hopf bifurcation can be obtained under the condition that for . Inserting the steady state , the critical feed rate is given by

Figure 2.1: Phase diagram of the inviscid dynamics. In region III, there is only one spatially uniform steady state , which is stable for all and . In region I and II , there are three steady states. In region II, only two steady states and are stable. However, the steady state is stable in region I, while the uniform state loses stability when is decreased through the Hopf bifurcation curve (the red dashed-line). The nontrivial fixed point is always unstable. The bifurcating periodic orbit is stable for and unstable for .

In the plane, we plot the phase diagram for the homogeneous fractional GS system. In Fig. LABEL:Fig._1, outside the region bounded by the saddle-node curve (region III), there only exists a single spatially uniform state (the trivial state ) that is stable for all , while in the remainder region (region I and II) of this plane, there are three spatially uniform steady states. In region II, the system is bistable with the steady states and . Since the coefficient is less than zero in region I, the steady state loses stability when is decreased through the Hopf bifurcation curve (the red dashed-line). The third steady state is always unstable for arbitrary parameter. In addition, the intersection point of the two curves is .

3 Well-posedness

The well-posedness of the GS model with classical diffusion has been analyzed in [32]. In this section, we consider the well-posedness of the fractional GS model. Firstly, we give some useful property and lemmas (see [8, 14, 23]).

Property 3.1

If , then

Lemma 3.1

For real , , if , then

Lemma 3.2

For real , then

We can extend the fractional derivatives to by zero outside easily. Using the above property, lemmas and the definition of the fractional Laplacian (LABEL:s2:ee1), we derive the following inner product formula i.e., for ,


where is the completely space of in . We omit the subscripts of the norm in this paper whenever it is clear from the context. Then the solution of the fractional GS model can be bounded as in the following theorem.

Theorem 3.3

Suppose the solutions of the fractional GS model belong to for any initial data . Then the following estimates hold


where is the solution of the following equation


The proof is presented in Appendix LABEL:adx1.

4 Numerical discretization and stability analysis

In this section, we use the C-N difference scheme for time discretization and use the weighted shifted Grünwald difference operator introduced in [26] to approximate the spatial fractional Laplacian operator. In addition, we apply the second-order implicit-explicit method to handle the nonlinear terms.

4.1 Numerical discretization

Let and be positive integers. We define the space steps and time step as and . Then we partition the space domain and time interval into the uniform mesh with and for and . We let , for , and introduce the following notations:


We consider the first equation of in problem . In time direction, we use the Taylor expansion as the first step ,


and for ,


where .

In space, we use the second-order Grünwald difference operators and to approximate the fractional diffusion operators and , i.e.,

where the above Grünwald difference operators are derived in [26] with ,

The coefficients are defined as follows

where .

Then, multiplying equation (LABEL:s1:e3) with , we obtain


where is the truncation error. We define

Therefore, the equation can be rewritten as


Using the Taylor expansion, we have


Adding to , we obtain the alternative direction iteration (ADI) scheme as follows


Similarly, we can obtain the discretization scheme for


where .

We replace and by the numerical approximations and , and obtain an ADI finite difference scheme for the model for ,


and for ,


We define the following matrices

Therefore, for , we can rewrite the numerical scheme into the matrix form as follows