Abstract

We present a new numerical method for solving nonlinear reaction-diffusion systems with cross-diffusion which are often taken as mathematical models for many applications in the biological, physical, and chemical sciences. The two-dimensional system is discretized by the local discontinuous Galerkin (LDG) method on unstructured triangular meshes associated with the piecewise linear finite element spaces, which can derive not only numerical solutions but also approximations for fluxes at the same time comparing with most of study work up to now which has derived numerical solutions only. Considering the stability requirement for the explicit scheme with strict time step restriction (), the implicit integration factor (IIF) method is employed for the temporal discretization so that the time step can be relaxed as . Moreover, the method allows us to compute element by element and avoids solving a global system of nonlinear algebraic equations as the standard implicit schemes do, which can reduce the computational cost greatly. Numerical simulations about the system with exact solution and the Brusselator model, which is a theoretical model for a type of autocatalytic chemical reaction, are conducted to confirm the expected accuracy, efficiency, and advantages of the proposed schemes.

1. Introduction

In 1952, Turing proposed the reaction-diffusion systems in the seminal paper [1], which constitute an essential basis to describe morphogenetic mechanisms. It was suggested that, in a reaction-diffusion system describing the interaction between two species, different diffusion rates can lead to the destabilization of a constant steady state, followed by the transition to a nonhomogeneous steady state. According to this result, the equilibrium of the nonlinear system is asymptotically stable in the absence of diffusion but unstable in the presence of diffusion, which is called Turing unstable [2, 3]. This mechanism, known as diffusion driven instability, leads to the pattern appearance. In [4], Levin and Segel added diffusion to a planktonic system and demonstrated that diffusion plays an important role in generating spatial patterns. And these phenomena of spatial patterns have also been reported in [5–7] (see [8] for an extensive review).

Most of the reaction-diffusion systems used to predict the formation of patterns assumed that the diffusion of each species depends only on the gradient of the density of the species itself. However, cross-diffusion terms should be introduced when the gradient of the density of one species induces not only a flux of the species itself but also a flux of another species, which was originally introduced in the context of population dynamics [9] and has now gained a renewed interest in diverse contexts like ecology [10–12], social systems [13], drift-diffusion in semiconductors [14–16], granular materials [17], and other references [18–21].

In this paper, we study the reaction-diffusion system with both self-diffusion and cross-diffusion:where are the two biological or physical species or even two chemical concentrations, describe the reaction kinetics, and is the diffusion constant matrix. Here, the diagonal elements are called self-diffusion coefficients; the nondiagonals are called cross-diffusion coefficients. The term takes into account the flux of , , induced by the gradient of species , and the other three terms are the same.

In addition to the above theoretical aspects, an important interest lies in the behavior of numerical approximations exhibiting spatial pattern. And there are many numerical schemes to simulate system (1) including the finite difference methods, spectral methods, finite volume methods, and finite element methods. However, finite difference methods and spectral methods are constrained with the complicated and irregular domain geometries. At this point, they are not so popular as finite volume methods and finite element methods. In [22–24], the finite volume method proposed by Andreianov et al. [25] was adopted for the numerical treatment of the reaction-cross-diffusion system, and the formation and identification of spatial patterns were studied. And in [26], two kinds of finite element methods, which contain variational multiscale element-free Galerkin (VMEFG) and local discontinuous Galerkin (LDG) methods proposed by Cockburn and Shu [27], were applied to discretize the space derivative of the system. And fourth-order exponential time differencing Runge-Kutta method [28, 29] has been employed for temporal discretization.

In this paper, we choose to pursue LDG method, where more general numerical fluxes than those in [26] are used, coupled with Krylov implicit integration factor (IIF) methods [30, 31] for temporal discretization which is based on the IIF method [32]. By applying this method, we can derive the numerical approximations not only for solutions but also for fluxes at the same time. What is more, we can relax the time step necessary for explicit schemes to . Moreover, the method allows us to compute element by element and avoids solving a global system of nonlinear algebraic equations as the standard implicit schemes do, which can reduce the computational cost greatly.

The rest of this paper is organized as follows. In Section 2, we present the LDG formulation for spatial discretization, eliminate the auxiliary variables at the element level, and then apply the second-order IIF methods to discretize the resulting ordinary differential equations (ODEs) which have only original variables as unknowns. In Section 3, some numerical experiments including the test system with exact solutions and the Brusselator model (see, e.g., [33, 34]) are conducted to show that the results obtained by our method agree well with those in [23, 26] and our method possesses its own advantages. Finally, some conclusions are drawn in Section 4.

2. Construction of the Fully Discrete Scheme

In this section, we present the fully discrete scheme, which was obtained by combining the LDG method with the IIF method, to solve the nonlinear reaction-diffusion system (1). Here we consider the system defined on , together with no-flux boundary conditions,and appropriate initial conditions,where is an open, bounded domain and is the outward unit normal to .

Let be regular triangulation of , and is the mesh size. denotes the collection of all edges in . and are the sets of interior edges and boundary edges, respectively.

2.1. The LDG Method for Spatial Discretization

Firstly, by introducing two auxiliary variables and , we rewrite (1) as the following first-order differential equations:

Define the finite element space as follows: where denotes the linear polynomials defined on element .

Then the semidiscrete LDG formulation can be defined as follows. For , find and such that, for and ,where is the outward unit normal vector to . The quantities and are the so-called numerical fluxes and are chosen as [35] The stability parameter is taken to be to enhance the accuracy of the LDG method. The auxiliary vector parameter is generally chosen as on each edge .

The boundary conditions (2) are imposed through the following definition of the numerical fluxes:The numerical fluxes and can be defined similar to and , respectively.

By use of basis functions, we express and , in element aswhere denotes the basis functions: and , , , , , are the corresponding degrees of freedom:

For element , let , , denote its three adjacent elements. And we employ the subscript to mark the quantities belonging to the adjacent elements (see Figure 1).

Figure 1 

The sketch of triangular element and its neighbor elements.

Then, substituting (7) into (6), and applying expressions (9), we can obtain the matrix form for interior element , where we do by the same way with that in [36]: where and denote the time derivatives and the matrices are calculated as follows: for , where , with denoting the common edge between element and its adjacent element . We also use the relations that , , and .

Remark 1. If the edge shared by and is a boundary edge, that is to say, does not exist, then the above matrices related to the edge can have some differences. We only need to insert the numerical fluxes on boundary edges (8) into (6). Then, by the same way, we can derive the matrices for the boundary edges. The quantities , , and , are not needed and should be gotten rid of. In addition,

To facilitate computations, we rewrite the above matrix form into two separate matrix equations: where we use the notifications

In this paper, we take the degrees of freedom as the values of midpoints at three edges of ; then where is the area of element , is the unit matrix, and At this moment, (15) can be rewritten as ; :which is an advantage of LDG method where the auxiliary variables can be expressed by original variables locally.

As a special case of (20), for the adjacent elements , , we have where and , , denote the quantities belonging to adjacent elements (see Figure 1). Similarly, , , are used to mark the quantities belonging to the adjacent elements of , respectively. In addition, and when .

Substituting the above equations and (20) into (16), we derive a system including original variables only: where