# Local Discontinuous Galerkin Methods Coupled with Implicit Integration Factor Methods for Solving Reaction-Cross-Diffusion Systems

Authors: Na An; Xijun Yu; Chaobao Huang; Maochang Duan
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## 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 (Δ
t
=
O
(
h
m
i
n
2
)), the implicit integration factor (IIF) method is employed for the temporal discretization so that the time step can be relaxed as Δ
t
=
O
(
h
m
i
n
). 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.


---

## Body

## 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:(1)
∂
u
∂
t
=
Δ
b
11
u
+
b
12
v
+
f
u
,
v
,
∂
v
∂
t
=
Δ
b
21
u
+
b
22
v
+
g
u
,
v
,where u
,
v are the two biological or physical species or even two chemical concentrations, f
,
g describe the reaction kinetics, and B
=
b
11
b
12
b
21
b
22 is the diffusion constant matrix. Here, the diagonal elements b
11
,
b
22 are called self-diffusion coefficients; the nondiagonals b
12
,
b
21 are called cross-diffusion coefficients. The term Δ
(
b
11
u
)
=
∇
·
(
∇
(
b
11
u
)
) takes into account the flux of u, ∇
(
b
11
u
), induced by the gradient of species u, 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 Δ
t
=
O
(
h
m
i
n
2
) necessary for explicit schemes to Δ
t
=
O
(
h
m
i
n
). 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 Section2, we present the LDG formulation for spatial discretization, eliminate the auxiliary variables q
h
,
p
h 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 u
h
,
v
h 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 Ω
×
[
0
,
T
], together with no-flux boundary conditions,(2)
∂
u
x
,
y
,
t
∂
n
=
0
,
∂
v
x
,
y
,
t
∂
n
=
0
,and appropriate initial conditions,(3)
u
x
,
y
,
0
=
u
0
x
,
y
,
v
x
,
y
,
0
=
v
0
x
,
y
,where Ω is an open, bounded domain Ω
⊂
R
2 and n is the outward unit normal to ∂
Ω.LetT
h
=
{
K
} be regular triangulation of Ω, and h is the mesh size. E
h denotes the collection of all edges in T
h. E
h
∘ and E
h
b are the sets of interior edges and boundary edges, respectively.


### 2.1. The LDG Method for Spatial Discretization

Firstly, by introducing two auxiliary variablesq
=
b
11
∇
u and p
=
b
22
∇
v, we rewrite (1) as the following first-order differential equations:(4)
q
-
b
11
∇
u
=
0
,
p
-
b
22
∇
v
=
0
,
∂
u
∂
t
-
∇
·
q
-
b
12
b
22
∇
·
p
=
f
u
,
v
,
∂
v
∂
t
-
∇
·
p
-
b
21
b
11
∇
·
q
=
g
u
,
v
.Define the finite element space as follows:(5)
V
h
=
r
h
∈
L
2
Ω
:
r
h
K
∈
P
1
K
,
∀
K
∈
T
h
,
W
h
=
w
h
∈
L
2
Ω
2
:
w
h
K
∈
P
1
K
2
,
∀
K
∈
T
h
,where P
1
(
K
) denotes the linear polynomials defined on element K.Then the semidiscrete LDG formulation can be defined as follows. Fort
∈
(
0
,
T
], find q
h
(
t
)
,
p
h
(
t
)
∈
W
h and u
h
(
t
)
,
v
h
(
t
)
∈
V
h such that, for ∀
w
h
u
,
w
h
v
∈
W
h and ∀
r
h
u
,
r
h
v
∈
V
h,(6)
∫
K
q
h
·
w
h
u
d
x
+
∫
K
u
h
∇
·
b
11
w
h
u
d
x
-
∫
∂
K
u
h
^
b
11
w
h
u
·
n
K
d
s
=
0
,
∫
K
p
h
·
w
h
v
d
x
+
∫
K
v
h
∇
·
b
22
w
h
v
d
x
-
∫
∂
K
v
h
^
b
22
w
h
v
·
n
K
d
s
=
0
,
∫
K
∂
u
h
∂
t
r
h
u
d
x
+
∫
K
q
h
·
∇
r
h
u
d
x
-
∫
∂
K
q
h
^
·
n
K
r
h
u
d
s
+
b
12
b
22
∫
K
p
h
·
∇
r
h
u
d
x
-
∫
∂
K
p
h
^
·
n
K
r
h
u
d
s
=
∫
K
f
u
h
,
v
h
r
h
u
d
x
,
∫
K
∂
v
h
∂
t
r
h
v
d
x
+
∫
K
p
h
·
∇
r
h
v
d
x
-
∫
∂
K
p
h
^
·
n
K
r
h
v
d
s
+
b
21
b
11
∫
K
q
h
·
∇
r
h
v
d
x
-
∫
∂
K
q
h
^
·
n
K
r
h
v
d
s
=
∫
K
g
u
h
,
v
h
r
h
v
d
x
,where n
K is the outward unit normal vector to ∂
K. The quantities u
h
^ and q
h
^ are the so-called numerical fluxes and are chosen as [35] (7)
u
h
^
e
=
u
h
+
C
12
·
u
h
,
q
h
^
e
=
q
h
-
C
11
u
h
-
C
12
q
h
,
∀
e
∈
E
h
∘
.The stability parameter C
11
>
0 is taken to be O
(
h
e
-
1
) to enhance the accuracy of the LDG method. The auxiliary vector parameter C
12 is generally chosen as C
12
·
n
e
=
O
(
1
) on each edge e.The boundary conditions (2) are imposed through the following definition of the numerical fluxes:(8)
u
h
^
e
=
u
h
,
q
h
^
e
=
0
,
∀
e
∈
E
h
b
.The numerical fluxes v
h
^ and p
h
^ can be defined similar to u
h
^ and q
h
^, respectively.By use of basis functions, we expressu
h
,
v
h and q
h
=
(
q
h
,
x
,
q
h
,
y
), p
h
=
(
p
h
,
x
,
p
h
,
y
) in element K as(9)
u
h
=
∑
i
=
1
3
u
i
t
ϕ
i
x
,
y
=
Φ
T
u
¯
,
v
h
=
∑
i
=
1
3
v
i
t
ϕ
i
x
,
y
=
Φ
T
v
¯
,
q
h
,
l
=
∑
i
=
1
3
q
l
,
i
t
ϕ
i
x
,
y
=
Φ
T
q
l
¯
,
p
h
,
l
=
∑
i
=
1
3
p
l
,
i
t
ϕ
i
x
,
y
=
Φ
T
p
l
¯
,
l
=
x
,
y
,where Φ denotes the basis functions: (10)
Φ
=
ϕ
1
x
,
y
ϕ
2
x
,
y
ϕ
3
x
,
y
,and u
¯, v
¯, q
l
¯, p
l
¯, l
=
x
,
y, are the corresponding degrees of freedom: (11)
u
¯
=
u
1
t
u
2
t
u
3
t
,
v
¯
=
v
1
t
v
2
t
v
3
t
,
q
l
¯
=
q
l
,
1
t
q
l
,
2
t
q
l
,
3
t
,
p
l
¯
=
p
l
,
1
t
p
l
,
2
t
p
l
,
3
t
.For elementK, let K
i, i
=
1,2
,
3, denote its three adjacent elements. And we employ the subscript (
N
B
) to mark the quantities belonging to the adjacent elements (see Figure 1).Figure 1
The sketch of triangular elementK and its neighbor elements.Then, substituting (7) into (6), and applying expressions (9), we can obtain the matrix form for interior element K, where we do by the same way with that in [36]: (12)
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
M
0
0
0
0
0
0
M
q
x
¯
t
q
y
¯
t
p
x
¯
t
p
y
¯
t
u
¯
t
v
¯
t
+
M
0
0
0
b
11
H
x
0
0
M
0
0
b
11
H
y
0
0
0
M
0
0
b
22
H
x
0
0
0
M
0
b
22
H
y
J
x
J
y
b
12
b
22
J
x
b
12
b
22
J
y
0
0
b
21
b
11
J
x
b
21
b
11
J
y
J
x
J
y
0
0
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
+
∑
i
=
1
3
0
0
0
0
b
11
H
x
,
i
0
0
0
0
0
b
11
H
y
,
i
0
0
0
0
0
0
b
22
H
x
,
i
0
0
0
0
0
b
22
H
y
,
i
J
x
,
i
J
y
,
i
b
12
b
22
J
x
,
i
b
12
b
22
J
y
,
i
G
u
,
i
b
12
b
22
G
u
,
i
b
21
b
11
J
x
,
i
b
21
b
11
J
y
,
i
J
x
,
i
J
y
,
i
b
21
b
11
G
u
,
i
G
u
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
+
∑
i
=
1
3
0
0
0
0
b
11
H
x
,
B
,
i
0
0
0
0
0
b
11
H
y
,
B
,
i
0
0
0
0
0
0
b
22
H
x
,
B
,
i
0
0
0
0
0
b
22
H
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
b
12
b
22
J
x
,
B
,
i
b
12
b
22
J
y
,
B
,
i
G
u
,
B
,
i
b
12
b
22
G
u
,
B
,
i
b
21
b
11
J
x
,
B
,
i
b
21
b
11
J
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
b
21
b
11
G
u
,
B
,
i
G
u
,
B
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
N
B
,
i
=
0
0
0
0
S
u
S
v
,where (
u
¯
)
t
=
∂
u
¯
/
∂
t and (
v
¯
)
t
=
∂
v
¯
/
∂
t denote the time derivatives and the matrices are calculated as follows: for m
,
n
=
1,2
,
3, (13)
M
m
n
=
∫
K
ϕ
m
ϕ
n
d
x
,
S
u
,
m
=
∫
K
f
u
h
,
v
h
ϕ
m
d
x
,
S
v
,
m
=
∫
K
g
u
h
,
v
h
ϕ
m
d
x
,
J
l
,
m
n
=
∫
K
∂
ϕ
m
∂
l
ϕ
n
d
x
,
H
l
,
m
n
=
J
l
,
m
n
=
∫
K
∂
ϕ
m
∂
l
ϕ
n
d
x
,
l
=
x
,
y
,
G
u
,
i
,
m
n
=
∫
∂
K
i
C
11
ϕ
m
ϕ
n
d
x
,
G
u
,
B
,
i
,
m
n
=
-
∫
∂
K
i
C
11
ϕ
m
ϕ
n
N
B
,
i
d
x
,
i
=
1,2
,
3
,
H
l
,
i
,
m
n
=
-
∫
∂
K
i
1
2
+
C
12
ϕ
m
ϕ
n
n
l
d
s
,
H
l
,
B
,
i
,
m
n
=
-
∫
∂
K
i
1
2
-
C
12
ϕ
m
ϕ
n
N
B
,
i
n
l
d
s
,
J
l
,
i
,
m
n
=
-
∫
∂
K
i
1
2
-
C
12
ϕ
m
ϕ
n
n
l
d
s
,
J
l
,
B
,
i
,
m
n
=
-
∫
∂
K
i
1
2
+
C
12
ϕ
m
ϕ
n
N
B
,
i
n
l
d
s
,where ∂
K
=
∪
i
=
1
3
(
∂
K
)
i, with (
∂
K
)
i denoting the common edge between element K and its adjacent element K
i. We also use the relations that n
=
n
K
=
(
n
x
,
n
y
), C
12
=
C
12
·
n, and n
(
N
B
)
=
-
n.Remark 1.
If the edge(
∂
K
)
i shared by K and K
i is a boundary edge, that is to say, K
i does not exist, then the above matrices related to the edge (
∂
K
)
i 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 G
u
,
B
,
i, H
l
,
B
,
i, and J
l
,
B
,
i
,
l
=
x
,
y, are not needed and should be gotten rid of. In addition, (14)
J
l
,
i
,
m
n
=
0
,
G
u
,
i
,
m
n
=
0
,
H
l
,
i
,
m
n
=
-
∫
∂
K
i
ϕ
m
ϕ
n
n
l
d
s
,
l
=
x
,
y
.To facilitate computations, we rewrite the above matrix form into two separate matrix equations:(15)
M
0
0
0
0
M
0
0
0
0
M
0
0
0
0
M
q
x
¯
q
y
¯
p
x
¯
p
y
¯
+
b
11
H
H
x
0
b
11
H
H
y
0
0
b
22
H
H
x
0
b
22
H
H
y
u
¯
v
¯
+
∑
i
=
1
3
b
11
H
x
,
B
,
i
0
b
11
H
y
,
B
,
i
0
0
b
22
H
x
,
B
,
i
0
b
22
H
y
,
B
,
i
u
¯
v
¯
N
B
,
i
=
0
0
0
0
,
(16)
M
0
0
M
u
¯
t
v
¯
t
+
J
J
x
J
J
y
b
12
b
22
J
J
x
b
12
b
22
J
J
y
b
21
b
11
J
J
x
b
21
b
11
J
J
y
J
J
x
J
J
y
q
x
¯
q
y
¯
p
x
¯
p
y
¯
+
G
U
b
12
b
22
G
U
b
21
b
11
G
U
G
U
u
¯
v
¯
+
∑
i
=
1
3
J
x
,
B
,
i
J
y
,
B
,
i
b
12
b
22
J
x
,
B
,
i
b
12
b
22
J
y
,
B
,
i
b
21
b
11
J
x
,
B
,
i
b
21
b
11
J
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
N
B
,
i
+
∑
i
=
1
3
G
u
,
B
,
i
b
12
b
22
G
u
,
B
,
i
b
21
b
11
G
u
,
B
,
i
G
u
,
B
,
i
u
¯
v
¯
N
B
,
i
=
S
u
S
v
,where we use the notifications (17)
H
H
x
=
H
x
+
∑
i
=
1
3
H
x
,
i
,
H
H
y
=
H
y
+
∑
i
=
1
3
H
y
,
i
,
J
J
x
=
J
x
+
∑
i
=
1
3
J
x
,
i
,
J
J
y
=
J
y
+
∑
i
=
1
3
J
y
,
i
,
G
U
=
∑
i
=
1
3
G
u
,
i
.In this paper, we take the degrees of freedom as the values of midpoints at three edges ofK; then (18)
M
≈
K
3
I
,
S
u
≈
K
3
f
u
¯
,
v
¯
,
S
v
≈
K
3
g
u
¯
,
v
¯
,where |
K
| is the area of element K, I is the unit matrix, and (19)
f
u
¯
,
v
¯
=
f
u
1
,
v
1
f
u
2
,
v
2
f
u
3
,
v
3
,
g
u
¯
,
v
¯
=
g
u
1
,
v
1
g
u
2
,
v
2
f
u
3
,
v
3
.At this moment, (15) can be rewritten as ∀
K
∈
T
h; l
=
x
,
y:(20)
q
l
¯
=
-
3
K
b
11
H
H
l
u
¯
+
∑
i
=
1
3
H
l
,
B
,
i
u
¯
N
B
,
i
,
p
l
¯
=
-
3
K
b
22
H
H
l
v
¯
+
∑
i
=
1
3
H
l
,
B
,
i
v
¯
N
B
,
i
,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 K
i, i
=
1,2
,
3, we have (21)
q
l
¯
N
B
,
i
=
-
3
K
i
b
11
H
H
l
N
B
,
i
u
¯
N
B
,
i
+
∑
j
=
1
3
H
l
,
B
,
j
N
B
,
i
u
¯
N
B
,
i
N
B
,
j
,
p
l
¯
N
B
,
i
=
-
3
K
i
b
22
H
H
l
N
B
,
i
v
¯
N
B
,
i
+
∑
j
=
1
3
H
l
,
B
,
j
N
B
,
i
v
¯
N
B
,
i
N
B
,
j
,where l
=
x
,
y and (
N
B
,
i
), i
=
1,2
,
3, denote the quantities belonging to K
′
s adjacent elements K
i (see Figure 1). Similarly, (
N
B
,
i
)
(
N
B
,
j
), j
=
1,2
,
3, are used to mark the quantities belonging to the adjacent elements K
i
j of K
i, respectively. In addition, u
¯
(
N
B
,
i
)
(
N
B
,
j
)
=
u
¯ and v
¯
(
N
B
,
i
)
(
N
B
,
j
)
=
v
¯ when i
=
j.Substituting the above equations and (20) into (16), we derive a system including original variables only: (22)
d
w
d
t
=
-
3
K
N
w
+
∑
i
=
1
3
N
¯
i
w
N
B
,
i
+
∑
i
=
1
3
∑
j
=
1
3
N
~
i
j
w
N
B
,
i
N
B
,
j
+
F
w
,where (23)
w
=
u
¯
v
¯
,
F
w
=
f
u
¯
,
v
¯
g
u
¯
,
v
¯
,
N
=
N
u
b
12
b
22
N
v
b
21
b
11
N
u
N
v
,
N
¯
i
=
N
¯
i
,
u
b
12
b
22
N
¯
i
,
v
b
21
b
11
N
¯
i
,
u
N
¯
i
,
v
,
N
~
i
j
=
N
~
i
j
,
u
b
12
b
22
N
~
i
j
,
v
b
21
b
11
N
~
i
j
,
u
N
~
i
j
,
v
,
i
,
j
=
1,2
,
3
,
N
u
=
G
U
-
3
K
b
11
J
J
x
H
H
x
+
J
J
y
H
H
y
,
N
v
=
G
U
-
3
K
b
22
J
J
x
H
H
x
+
J
J
y
H
H
y
,
N
¯
i
,
u
=
G
u
,
B
,
i
-
3
K
b
11
J
J
x
H
x
,
B
,
i
+
J
J
y
H
y
,
B
,
i
-
3
K
i
b
11
J
x
,
B
,
i
H
H
x
N
B
,
i
+
J
y
,
B
,
i
H
H
y
N
B
,
i
,
N
¯
i
,
v
=
G
u
,
B
,
i
-
3
K
b
22
J
J
x
H
x
,
B
,
i
+
J
J
y
H
y
,
B
,
i
-
3
K
i
b
22
J
x
,
B
,
i
H
H
x
N
B
,
i
+
J
y
,
B
,
i
H
H
y
N
B
,
i
,
N
~
i
,
j
,
u
=
-
3
K
i
b
11
J
x
,
B
,
i
H
x
,
B
,
j
N
B
,
i
+
J
y
,
B
,
i
H
y
,
B
,
j
N
B
,
i
,
N
~
i
,
j
,
v
=
-
3
K
i
b
22
J
x
,
B
,
i
H
x
,
B
,
j
N
B
,
i
+
J
y
,
B
,
i
H
y
,
B
,
j
N
B
,
i
.


### 2.2. The IIF Methods Based on Krylov Subspace Approximation for Temporal Discretization

Assembling (22) over all of the elements in T
h, we derive a global system of ODEs:(24)
d
W
d
t
=
A
W
+
F
W
,where W
=
(
w
1
T
,
w
2
T
,
…
,
w
N
e
T
)
T, F
(
W
)
=
(
F
(
w
1
)
T
,
F
(
w
2
)
T
,
…
,
F
(
w
N
e
)
T
)
T, w
j is the degrees of freedom on K
j
,
j
=
1,2
,
…
,
N
e, and here N
e denotes the number of triangular elements. The 6
N
e
×
6
N
e global matrix A is sparse and formulated element by element according to (22). Each element K
∈
T
h contributes to the global matrix A with no more than ten 6
×
6 block matrices at corresponding six rows of A.Then we apply the second-order IIF scheme to system (24): (25)
W
n
+
1
=
e
A
Δ
t
W
n
+
Δ
t
2
F
W
n
+
Δ
t
2
F
W
n
+
1
,where n is the time level, t
n
+
1
=
t
n
+
Δ
t, and W
n
=
W
(
t
n
).When we computeW
n
+
1, the vector Q
=
e
A
Δ
t
(
W
n
+
(
Δ
t
/
2
)
F
(
W
n
)
) is a known quantity related to the earlier time level and can be computed by the Krylov subspace approximation shown in [30]. The nonlinear system at t
n
+
1 is decoupled from the diffusion term with a simple form:(26)
W
n
+
1
=
Q
+
Δ
t
2
F
W
n
+
1
,which can be solved element by element. And the N
e
6
×
6 systems are independent of each other with every system of the same structure. The local algebraic system on every element K
j
,
j
=
1,2
,
…
,
N
e, is of the following form: (27)
R
w
j
n
+
1
=
0with (28)
w
j
n
+
1
=
u
¯
j
n
+
1
v
¯
j
n
+
1
=
u
j
,
1
n
+
1
u
j
,
2
n
+
1
u
j
,
3
n
+
1
v
j
,
1
n
+
1
v
j
,
2
n
+
1
v
j
,
3
n
+
1
,
R
w
j
n
+
1
=
u
j
,
1
n
+
1
-
Q
6
j
-
1
+
1
-
Δ
t
2
f
u
j
,
1
n
+
1
,
v
j
,
1
n
+
1
u
j
,
2
n
+
1
-
Q
6
j
-
1
+
2
-
Δ
t
2
f
u
j
,
2
n
+
1
,
v
j
,
2
n
+
1
u
j
,
3
n
+
1
-
Q
6
j
-
1
+
3
-
Δ
t
2
f
u
j
,
3
n
+
1
,
v
j
,
3
n
+
1
v
j
,
1
n
+
1
-
Q
6
j
-
1
+
4
-
Δ
t
2
g
u
j
,
1
n
+
1
,
v
j
,
1
n
+
1
v
j
,
2
n
+
1
-
Q
6
j
-
1
+
5
-
Δ
t
2
g
u
j
,
2
n
+
1
,
v
j
,
2
n
+
1
v
j
,
3
n
+
1
-
Q
6
j
-
1
+
6
-
Δ
t
2
g
u
j
,
3
n
+
1
,
v
j
,
3
n
+
1
.The Newton iterative method can be applied to solve the above system. In the iterations to compute w
j
n
+
1, we use the numerical value w
j
n at time t
n as the initial guess. The threshold value for judging Newton iteration can be set small enough and is taken as 1
0
-
13 in the numerical examples.


## 2.1. The LDG Method for Spatial Discretization

Firstly, by introducing two auxiliary variablesq
=
b
11
∇
u and p
=
b
22
∇
v, we rewrite (1) as the following first-order differential equations:(4)
q
-
b
11
∇
u
=
0
,
p
-
b
22
∇
v
=
0
,
∂
u
∂
t
-
∇
·
q
-
b
12
b
22
∇
·
p
=
f
u
,
v
,
∂
v
∂
t
-
∇
·
p
-
b
21
b
11
∇
·
q
=
g
u
,
v
.Define the finite element space as follows:(5)
V
h
=
r
h
∈
L
2
Ω
:
r
h
K
∈
P
1
K
,
∀
K
∈
T
h
,
W
h
=
w
h
∈
L
2
Ω
2
:
w
h
K
∈
P
1
K
2
,
∀
K
∈
T
h
,where P
1
(
K
) denotes the linear polynomials defined on element K.Then the semidiscrete LDG formulation can be defined as follows. Fort
∈
(
0
,
T
], find q
h
(
t
)
,
p
h
(
t
)
∈
W
h and u
h
(
t
)
,
v
h
(
t
)
∈
V
h such that, for ∀
w
h
u
,
w
h
v
∈
W
h and ∀
r
h
u
,
r
h
v
∈
V
h,(6)
∫
K
q
h
·
w
h
u
d
x
+
∫
K
u
h
∇
·
b
11
w
h
u
d
x
-
∫
∂
K
u
h
^
b
11
w
h
u
·
n
K
d
s
=
0
,
∫
K
p
h
·
w
h
v
d
x
+
∫
K
v
h
∇
·
b
22
w
h
v
d
x
-
∫
∂
K
v
h
^
b
22
w
h
v
·
n
K
d
s
=
0
,
∫
K
∂
u
h
∂
t
r
h
u
d
x
+
∫
K
q
h
·
∇
r
h
u
d
x
-
∫
∂
K
q
h
^
·
n
K
r
h
u
d
s
+
b
12
b
22
∫
K
p
h
·
∇
r
h
u
d
x
-
∫
∂
K
p
h
^
·
n
K
r
h
u
d
s
=
∫
K
f
u
h
,
v
h
r
h
u
d
x
,
∫
K
∂
v
h
∂
t
r
h
v
d
x
+
∫
K
p
h
·
∇
r
h
v
d
x
-
∫
∂
K
p
h
^
·
n
K
r
h
v
d
s
+
b
21
b
11
∫
K
q
h
·
∇
r
h
v
d
x
-
∫
∂
K
q
h
^
·
n
K
r
h
v
d
s
=
∫
K
g
u
h
,
v
h
r
h
v
d
x
,where n
K is the outward unit normal vector to ∂
K. The quantities u
h
^ and q
h
^ are the so-called numerical fluxes and are chosen as [35] (7)
u
h
^
e
=
u
h
+
C
12
·
u
h
,
q
h
^
e
=
q
h
-
C
11
u
h
-
C
12
q
h
,
∀
e
∈
E
h
∘
.The stability parameter C
11
>
0 is taken to be O
(
h
e
-
1
) to enhance the accuracy of the LDG method. The auxiliary vector parameter C
12 is generally chosen as C
12
·
n
e
=
O
(
1
) on each edge e.The boundary conditions (2) are imposed through the following definition of the numerical fluxes:(8)
u
h
^
e
=
u
h
,
q
h
^
e
=
0
,
∀
e
∈
E
h
b
.The numerical fluxes v
h
^ and p
h
^ can be defined similar to u
h
^ and q
h
^, respectively.By use of basis functions, we expressu
h
,
v
h and q
h
=
(
q
h
,
x
,
q
h
,
y
), p
h
=
(
p
h
,
x
,
p
h
,
y
) in element K as(9)
u
h
=
∑
i
=
1
3
u
i
t
ϕ
i
x
,
y
=
Φ
T
u
¯
,
v
h
=
∑
i
=
1
3
v
i
t
ϕ
i
x
,
y
=
Φ
T
v
¯
,
q
h
,
l
=
∑
i
=
1
3
q
l
,
i
t
ϕ
i
x
,
y
=
Φ
T
q
l
¯
,
p
h
,
l
=
∑
i
=
1
3
p
l
,
i
t
ϕ
i
x
,
y
=
Φ
T
p
l
¯
,
l
=
x
,
y
,where Φ denotes the basis functions: (10)
Φ
=
ϕ
1
x
,
y
ϕ
2
x
,
y
ϕ
3
x
,
y
,and u
¯, v
¯, q
l
¯, p
l
¯, l
=
x
,
y, are the corresponding degrees of freedom: (11)
u
¯
=
u
1
t
u
2
t
u
3
t
,
v
¯
=
v
1
t
v
2
t
v
3
t
,
q
l
¯
=
q
l
,
1
t
q
l
,
2
t
q
l
,
3
t
,
p
l
¯
=
p
l
,
1
t
p
l
,
2
t
p
l
,
3
t
.For elementK, let K
i, i
=
1,2
,
3, denote its three adjacent elements. And we employ the subscript (
N
B
) to mark the quantities belonging to the adjacent elements (see Figure 1).Figure 1
The sketch of triangular elementK and its neighbor elements.Then, substituting (7) into (6), and applying expressions (9), we can obtain the matrix form for interior element K, where we do by the same way with that in [36]: (12)
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
M
0
0
0
0
0
0
M
q
x
¯
t
q
y
¯
t
p
x
¯
t
p
y
¯
t
u
¯
t
v
¯
t
+
M
0
0
0
b
11
H
x
0
0
M
0
0
b
11
H
y
0
0
0
M
0
0
b
22
H
x
0
0
0
M
0
b
22
H
y
J
x
J
y
b
12
b
22
J
x
b
12
b
22
J
y
0
0
b
21
b
11
J
x
b
21
b
11
J
y
J
x
J
y
0
0
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
+
∑
i
=
1
3
0
0
0
0
b
11
H
x
,
i
0
0
0
0
0
b
11
H
y
,
i
0
0
0
0
0
0
b
22
H
x
,
i
0
0
0
0
0
b
22
H
y
,
i
J
x
,
i
J
y
,
i
b
12
b
22
J
x
,
i
b
12
b
22
J
y
,
i
G
u
,
i
b
12
b
22
G
u
,
i
b
21
b
11
J
x
,
i
b
21
b
11
J
y
,
i
J
x
,
i
J
y
,
i
b
21
b
11
G
u
,
i
G
u
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
+
∑
i
=
1
3
0
0
0
0
b
11
H
x
,
B
,
i
0
0
0
0
0
b
11
H
y
,
B
,
i
0
0
0
0
0
0
b
22
H
x
,
B
,
i
0
0
0
0
0
b
22
H
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
b
12
b
22
J
x
,
B
,
i
b
12
b
22
J
y
,
B
,
i
G
u
,
B
,
i
b
12
b
22
G
u
,
B
,
i
b
21
b
11
J
x
,
B
,
i
b
21
b
11
J
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
b
21
b
11
G
u
,
B
,
i
G
u
,
B
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
u
¯
v
¯
N
B
,
i
=
0
0
0
0
S
u
S
v
,where (
u
¯
)
t
=
∂
u
¯
/
∂
t and (
v
¯
)
t
=
∂
v
¯
/
∂
t denote the time derivatives and the matrices are calculated as follows: for m
,
n
=
1,2
,
3, (13)
M
m
n
=
∫
K
ϕ
m
ϕ
n
d
x
,
S
u
,
m
=
∫
K
f
u
h
,
v
h
ϕ
m
d
x
,
S
v
,
m
=
∫
K
g
u
h
,
v
h
ϕ
m
d
x
,
J
l
,
m
n
=
∫
K
∂
ϕ
m
∂
l
ϕ
n
d
x
,
H
l
,
m
n
=
J
l
,
m
n
=
∫
K
∂
ϕ
m
∂
l
ϕ
n
d
x
,
l
=
x
,
y
,
G
u
,
i
,
m
n
=
∫
∂
K
i
C
11
ϕ
m
ϕ
n
d
x
,
G
u
,
B
,
i
,
m
n
=
-
∫
∂
K
i
C
11
ϕ
m
ϕ
n
N
B
,
i
d
x
,
i
=
1,2
,
3
,
H
l
,
i
,
m
n
=
-
∫
∂
K
i
1
2
+
C
12
ϕ
m
ϕ
n
n
l
d
s
,
H
l
,
B
,
i
,
m
n
=
-
∫
∂
K
i
1
2
-
C
12
ϕ
m
ϕ
n
N
B
,
i
n
l
d
s
,
J
l
,
i
,
m
n
=
-
∫
∂
K
i
1
2
-
C
12
ϕ
m
ϕ
n
n
l
d
s
,
J
l
,
B
,
i
,
m
n
=
-
∫
∂
K
i
1
2
+
C
12
ϕ
m
ϕ
n
N
B
,
i
n
l
d
s
,where ∂
K
=
∪
i
=
1
3
(
∂
K
)
i, with (
∂
K
)
i denoting the common edge between element K and its adjacent element K
i. We also use the relations that n
=
n
K
=
(
n
x
,
n
y
), C
12
=
C
12
·
n, and n
(
N
B
)
=
-
n.Remark 1.
If the edge(
∂
K
)
i shared by K and K
i is a boundary edge, that is to say, K
i does not exist, then the above matrices related to the edge (
∂
K
)
i 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 G
u
,
B
,
i, H
l
,
B
,
i, and J
l
,
B
,
i
,
l
=
x
,
y, are not needed and should be gotten rid of. In addition, (14)
J
l
,
i
,
m
n
=
0
,
G
u
,
i
,
m
n
=
0
,
H
l
,
i
,
m
n
=
-
∫
∂
K
i
ϕ
m
ϕ
n
n
l
d
s
,
l
=
x
,
y
.To facilitate computations, we rewrite the above matrix form into two separate matrix equations:(15)
M
0
0
0
0
M
0
0
0
0
M
0
0
0
0
M
q
x
¯
q
y
¯
p
x
¯
p
y
¯
+
b
11
H
H
x
0
b
11
H
H
y
0
0
b
22
H
H
x
0
b
22
H
H
y
u
¯
v
¯
+
∑
i
=
1
3
b
11
H
x
,
B
,
i
0
b
11
H
y
,
B
,
i
0
0
b
22
H
x
,
B
,
i
0
b
22
H
y
,
B
,
i
u
¯
v
¯
N
B
,
i
=
0
0
0
0
,
(16)
M
0
0
M
u
¯
t
v
¯
t
+
J
J
x
J
J
y
b
12
b
22
J
J
x
b
12
b
22
J
J
y
b
21
b
11
J
J
x
b
21
b
11
J
J
y
J
J
x
J
J
y
q
x
¯
q
y
¯
p
x
¯
p
y
¯
+
G
U
b
12
b
22
G
U
b
21
b
11
G
U
G
U
u
¯
v
¯
+
∑
i
=
1
3
J
x
,
B
,
i
J
y
,
B
,
i
b
12
b
22
J
x
,
B
,
i
b
12
b
22
J
y
,
B
,
i
b
21
b
11
J
x
,
B
,
i
b
21
b
11
J
y
,
B
,
i
J
x
,
B
,
i
J
y
,
B
,
i
q
x
¯
q
y
¯
p
x
¯
p
y
¯
N
B
,
i
+
∑
i
=
1
3
G
u
,
B
,
i
b
12
b
22
G
u
,
B
,
i
b
21
b
11
G
u
,
B
,
i
G
u
,
B
,
i
u
¯
v
¯
N
B
,
i
=
S
u
S
v
,where we use the notifications (17)
H
H
x
=
H
x
+
∑
i
=
1
3
H
x
,
i
,
H
H
y
=
H
y
+
∑
i
=
1
3
H
y
,
i
,
J
J
x
=
J
x
+
∑
i
=
1
3
J
x
,
i
,
J
J
y
=
J
y
+
∑
i
=
1
3
J
y
,
i
,
G
U
=
∑
i
=
1
3
G
u
,
i
.In this paper, we take the degrees of freedom as the values of midpoints at three edges ofK; then (18)
M
≈
K
3
I
,
S
u
≈
K
3
f
u
¯
,
v
¯
,
S
v
≈
K
3
g
u
¯
,
v
¯
,where |
K
| is the area of element K, I is the unit matrix, and (19)
f
u
¯
,
v
¯
=
f
u
1
,
v
1
f
u
2
,
v
2
f
u
3
,
v
3
,
g
u
¯
,
v
¯
=
g
u
1
,
v
1
g
u
2
,
v
2
f
u
3
,
v
3
.At this moment, (15) can be rewritten as ∀
K
∈
T
h; l
=
x
,
y:(20)
q
l
¯
=
-
3
K
b
11
H
H
l
u
¯
+
∑
i
=
1
3
H
l
,
B
,
i
u
¯
N
B
,
i
,
p
l
¯
=
-
3
K
b
22
H
H
l
v
¯
+
∑
i
=
1
3
H
l
,
B
,
i
v
¯
N
B
,
i
,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 K
i, i
=
1,2
,
3, we have (21)
q
l
¯
N
B
,
i
=
-
3
K
i
b
11
H
H
l
N
B
,
i
u
¯
N
B
,
i
+
∑
j
=
1
3
H
l
,
B
,
j
N
B
,
i
u
¯
N
B
,
i
N
B
,
j
,
p
l
¯
N
B
,
i
=
-
3
K
i
b
22
H
H
l
N
B
,
i
v
¯
N
B
,
i
+
∑
j
=
1
3
H
l
,
B
,
j
N
B
,
i
v
¯
N
B
,
i
N
B
,
j
,where l
=
x
,
y and (
N
B
,
i
), i
=
1,2
,
3, denote the quantities belonging to K
′
s adjacent elements K
i (see Figure 1). Similarly, (
N
B
,
i
)
(
N
B
,
j
), j
=
1,2
,
3, are used to mark the quantities belonging to the adjacent elements K
i
j of K
i, respectively. In addition, u
¯
(
N
B
,
i
)
(
N
B
,
j
)
=
u
¯ and v
¯
(
N
B
,
i
)
(
N
B
,
j
)
=
v
¯ when i
=
j.Substituting the above equations and (20) into (16), we derive a system including original variables only: (22)
d
w
d
t
=
-
3
K
N
w
+
∑
i
=
1
3
N
¯
i
w
N
B
,
i
+
∑
i
=
1
3
∑
j
=
1
3
N
~
i
j
w
N
B
,
i
N
B
,
j
+
F
w
,where (23)
w
=
u
¯
v
¯
,
F
w
=
f
u
¯
,
v
¯
g
u
¯
,
v
¯
,
N
=
N
u
b
12
b
22
N
v
b
21
b
11
N
u
N
v
,
N
¯
i
=
N
¯
i
,
u
b
12
b
22
N
¯
i
,
v
b
21
b
11
N
¯
i
,
u
N
¯
i
,
v
,
N
~
i
j
=
N
~
i
j
,
u
b
12
b
22
N
~
i
j
,
v
b
21
b
11
N
~
i
j
,
u
N
~
i
j
,
v
,
i
,
j
=
1,2
,
3
,
N
u
=
G
U
-
3
K
b
11
J
J
x
H
H
x
+
J
J
y
H
H
y
,
N
v
=
G
U
-
3
K
b
22
J
J
x
H
H
x
+
J
J
y
H
H
y
,
N
¯
i
,
u
=
G
u
,
B
,
i
-
3
K
b
11
J
J
x
H
x
,
B
,
i
+
J
J
y
H
y
,
B
,
i
-
3
K
i
b
11
J
x
,
B
,
i
H
H
x
N
B
,
i
+
J
y
,
B
,
i
H
H
y
N
B
,
i
,
N
¯
i
,
v
=
G
u
,
B
,
i
-
3
K
b
22
J
J
x
H
x
,
B
,
i
+
J
J
y
H
y
,
B
,
i
-
3
K
i
b
22
J
x
,
B
,
i
H
H
x
N
B
,
i
+
J
y
,
B
,
i
H
H
y
N
B
,
i
,
N
~
i
,
j
,
u
=
-
3
K
i
b
11
J
x
,
B
,
i
H
x
,
B
,
j
N
B
,
i
+
J
y
,
B
,
i
H
y
,
B
,
j
N
B
,
i
,
N
~
i
,
j
,
v
=
-
3
K
i
b
22
J
x
,
B
,
i
H
x
,
B
,
j
N
B
,
i
+
J
y
,
B
,
i
H
y
,
B
,
j
N
B
,
i
.


## 2.2. The IIF Methods Based on Krylov Subspace Approximation for Temporal Discretization

Assembling (22) over all of the elements in T
h, we derive a global system of ODEs:(24)
d
W
d
t
=
A
W
+
F
W
,where W
=
(
w
1
T
,
w
2
T
,
…
,
w
N
e
T
)
T, F
(
W
)
=
(
F
(
w
1
)
T
,
F
(
w
2
)
T
,
…
,
F
(
w
N
e
)
T
)
T, w
j is the degrees of freedom on K
j
,
j
=
1,2
,
…
,
N
e, and here N
e denotes the number of triangular elements. The 6
N
e
×
6
N
e global matrix A is sparse and formulated element by element according to (22). Each element K
∈
T
h contributes to the global matrix A with no more than ten 6
×
6 block matrices at corresponding six rows of A.Then we apply the second-order IIF scheme to system (24): (25)
W
n
+
1
=
e
A
Δ
t
W
n
+
Δ
t
2
F
W
n
+
Δ
t
2
F
W
n
+
1
,where n is the time level, t
n
+
1
=
t
n
+
Δ
t, and W
n
=
W
(
t
n
).When we computeW
n
+
1, the vector Q
=
e
A
Δ
t
(
W
n
+
(
Δ
t
/
2
)
F
(
W
n
)
) is a known quantity related to the earlier time level and can be computed by the Krylov subspace approximation shown in [30]. The nonlinear system at t
n
+
1 is decoupled from the diffusion term with a simple form:(26)
W
n
+
1
=
Q
+
Δ
t
2
F
W
n
+
1
,which can be solved element by element. And the N
e
6
×
6 systems are independent of each other with every system of the same structure. The local algebraic system on every element K
j
,
j
=
1,2
,
…
,
N
e, is of the following form: (27)
R
w
j
n
+
1
=
0with (28)
w
j
n
+
1
=
u
¯
j
n
+
1
v
¯
j
n
+
1
=
u
j
,
1
n
+
1
u
j
,
2
n
+
1
u
j
,
3
n
+
1
v
j
,
1
n
+
1
v
j
,
2
n
+
1
v
j
,
3
n
+
1
,
R
w
j
n
+
1
=
u
j
,
1
n
+
1
-
Q
6
j
-
1
+
1
-
Δ
t
2
f
u
j
,
1
n
+
1
,
v
j
,
1
n
+
1
u
j
,
2
n
+
1
-
Q
6
j
-
1
+
2
-
Δ
t
2
f
u
j
,
2
n
+
1
,
v
j
,
2
n
+
1
u
j
,
3
n
+
1
-
Q
6
j
-
1
+
3
-
Δ
t
2
f
u
j
,
3
n
+
1
,
v
j
,
3
n
+
1
v
j
,
1
n
+
1
-
Q
6
j
-
1
+
4
-
Δ
t
2
g
u
j
,
1
n
+
1
,
v
j
,
1
n
+
1
v
j
,
2
n
+
1
-
Q
6
j
-
1
+
5
-
Δ
t
2
g
u
j
,
2
n
+
1
,
v
j
,
2
n
+
1
v
j
,
3
n
+
1
-
Q
6
j
-
1
+
6
-
Δ
t
2
g
u
j
,
3
n
+
1
,
v
j
,
3
n
+
1
.The Newton iterative method can be applied to solve the above system. In the iterations to compute w
j
n
+
1, we use the numerical value w
j
n at time t
n as the initial guess. The threshold value for judging Newton iteration can be set small enough and is taken as 1
0
-
13 in the numerical examples.


## 3. Numerical Experiments

In this section, numerical experiments are presented to demonstrate the validity and accuracy of the LDG method with the IIF scheme for solving the reaction-cross-diffusion system on two-dimensional triangular meshes. Firstly, we give a test example with exact solutions to manifest the spatial accuracy of the method. Then we apply the method to the Brusselator model with cross-diffusion. And numerical results agree well with those in other references [23, 26]. In addition, our Krylov IIF method for temporal discretization reduces the computational cost greatly, and LDG method for spatial discreetization derives the numerical approximations not only for solutions but also for fluxes at the same time.All of the numerical examples considered in this section are subject to no-flux boundary conditions (2). The triangular partitions used here are Delaunay partitions gotten from EasyMesh (see Figure 2). And the time step size is taken as(29)
Δ
t
=
c
f
l
∗
h
m
i
n
,where h
m
i
n is the length of the minimum edge in the triangular partition.Figure 2
Delaunay partitions for square domain in Example2 (a) and circular domain in Example 4 (b).
(a)
Partition for[
0,2
π
]
2 with h
=
0.1261
(b)
Partition for circular domian withh
=
0.4962In addition, the auxiliary parameters in the numerical fluxes (7) and (8) are taken as (30)
C
11
=
C
~
h
e
,
C
12
·
n
e
=
1
2
sign
⁡
n
e
1
+
n
e
2
,
∀
e
∈
∂
K
∩
E
h
,where h
e is the length of edge e; C
~
>
0 is the penalization parameter and is set to be C
~
=
1 in the following computation. n
e
=
(
n
e
1
,
n
e
2
) is the outward unit normal vector of K on e.


### 3.1. A Test Problem with Exact Solutions

To validate the spatial accuracy of our method numerically, we firstly consider a simple auxiliary reaction-diffusion system given in [23].Example 2.
The test reaction-diffusion system with cross-diffusion is of the following form:(31)
∂
u
∂
t
-
Δ
b
11
u
+
b
12
v
=
b
12
v
,
∂
v
∂
t
-
Δ
b
21
u
+
b
22
v
=
4
b
21
u
,defined on square Ω
=
[
0,2
π
]
2 and here b
11
=
b
22
=
1, b
12
=
1.5, and b
21
=
0.5. The exact solution is (32)
u
x
,
y
,
t
=
exp
⁡
-
4
t
cos
⁡
2
x
+
cos
⁡
2
y
,
v
x
,
y
,
t
=
exp
⁡
-
t
cos
⁡
x
+
cos
⁡
y
,from which we can derive the initial conditions.The simulation for Example2 is carried up to T
=
0.5 with c
f
l
=
0.05. The errors in L
2-norm and L
∞-norm are measured for both solutions and the fluxes on various triangular meshes. And the results are displayed in Table 1, where N
n denotes the number of triangular vertices. We can observe optimal convergence rates of O
(
h
2
) for solutions u
,
v in both L
2-norm and L
∞-norm and suboptimal convergence rates of O
(
h
) for fluxes q
,
p. These rates are consistent with theoretical results for the LDG method [35].Table 1
Numerical results for Example2 at T
=
0.5.
N
e
N
n
u
-
u
h
0
Order
u
-
u
h
∞
Order
v
-
v
h
0
Order
v
-
v
h
∞
Order
42
30
2.9159E − 1
—
1.0756E − 1
—
3.9750E − 1
—
1.1964E − 1
—
150
92
1.0714E − 1
1.5449
3.7922E − 2
1.6085
1.1495E − 1
1.9144
4.6267E − 2
1.4659
596
331
2.5577E − 2
2.0136
1.1788E − 2
1.6425
2.8564E − 2
1.9572
1.3543E − 2
1.7271
2392
1261
6.1548E − 3
2.2018
3.2021E − 3
2.0145
6.7407E − 3
2.2320
3.6577E − 3
2.0233
9224
4739
1.5429E − 3
2.0212
8.5707E − 4
1.9255
1.6793E − 3
2.0304
9.7541E − 4
1.9309
36838
18672
3.7359E − 4
1.9233
2.2556E − 4
1.8104
4.0320E − 4
1.9348
2.5540E − 4
1.8173
h
CPU (s)
q
-
q
h
0
Order
q
-
q
h
∞
Order
p
-
p
h
0
Order
p
-
p
h
∞
Order
1.8596
0.22
5.5622E − 1
—
1.6628E − 1
—
3.9209E − 1
—
1.2379E − 1
—
0.9727
0.28
1.8244E − 1
1.7200
5.4618E − 2
1.7178
1.9275E − 1
1.0957
7.6123E − 2
0.7502
0.4775
1.34
7.1312E − 2
1.3205
2.6911E − 2
0.9950
7.8795E − 2
1.2575
3.2105E − 2
1.2136
0.2501
9.75
3.1589E − 2
1.2586
1.3518E − 2
1.0643
3.5300E − 2
1.2412
1.5165E − 2
1.1593
0.1261
86.75
1.5204E − 2
1.0683
6.9372E − 3
0.9746
1.7065E − 2
1.0618
7.5730E − 3
1.0145
0.0603
473.50
7.4107E − 3
0.9746
3.3396E − 3
0.9914
8.2970E − 3
0.9780
3.8574E − 3
0.9148Moreover, Figure3 demonstrates graphs of numerical solutions u
h (left) and v
h (right) for this example, which are derived under the Delaunay partition shown in Figure 2(a). In addition, the numerical approximates for fluxes q
h (top) and p
h (bottom) are plotted in Figure 4.Figure 3
Graphs of numerical solutionsu
h and v
h at T
=
0.5 under the partition shown in Figure 2(a).
(a)
u
h
(b)
v
hFigure 4
Numerical approximates for fluxesq (top) and p (bottom) at T
=
0.5 under the partition shown in Figure 2(a).
(a)
x-direction component of q
h
(b)
y-direction component of q
h
(c)
x-direction component of p
h
(d)
y-direction component of p
h


### 3.2. Application to the Brusselator Model

In this subsection, we apply the LDG method coupled with the IIF scheme to solve the Brusselator model with cross-diffusion:(33)
∂
u
∂
t
-
Δ
b
11
u
+
b
12
v
=
-
β
+
1
u
+
u
2
v
+
α
,
∂
v
∂
t
-
Δ
b
21
u
+
b
22
v
=
β
u
-
u
2
v
,with the initial conditions chosen as small random perturbations of the equilibrium:(34)
u
x
,
0
=
5.8
+
1
3
rand
⁡
x
,
v
x
,
0
=
0.13
+
1
10
rand
⁡
x
,where x
=
(
x
,
y
)
∈
Ω, and rand: Ω
→
[
0,1
] is a random function in Fortran. In the following simulations, we take Ω as the square domain Ω
=
[
0,20
]
2 and circular domain Ω
=
{
(
x
,
y
)
∣
x
2
+
y
2
≤
1
0
2
}, respectively.Similar to [23], we set the parameters in the Brusselator model as(35)
b
11
=
0.4
,
b
22
=
2
,
b
21
=
0.02
,
α
=
6
,
β
=
1
,with b
12
=
24 and 32, respectively, based on which patterns are expected to appear. Actually, according to Theorem  2.3 in [23], the choice (35) is sufficient for the positive equilibrium point being linearly unstable and b
12 is taken as the Turing bifurcation parameter. Then the threshold b
12
=
22.2665 is obtained from (2.5) in [23]. Furthermore, it is observed that the increase of b
12 over the threshold with b
12
=
24 and b
12
=
32 yields the formation of spotted and labyrinthine patterns.Example 3.
We solve the Brusselator model (33)-(34) on the square domain Ω
=
[
0,20
]
2. In the computation, the square is divided into 16370 triangles and 8354 vertices, and now h
m
i
n
=
0.1492; mesh size h
=
0.2930. The time step size is taken as (29) with c
f
l
=
0.05.Different patterns will be obtained by selecting two sets of values for parametersb
12. The first set (b
12
=
24) leads to a spotted pattern as shown in Figure 5. It presents numerical approximations for solutions u (left) and v (right) at different values of final time T
=
1, T
=
10, and T
=
25, respectively. In addition, our method derives the approximations for fluxes q
=
b
11
∇
u and p
=
b
22
∇
v at the same time, and the graphs for q
h are plotted in Figure 6. The second set (b
12
=
32) generates a labyrinthine pattern (see Figure 7). The numerical approximations for fluxes q and p can also be obtained, and here we give only graphs of q
h in Figure 8. We observe a larger amplitude of the patterns (higher gradients) for both species.Figure 5
The time evolution of numerical solutionsu
h (left) and v
h (right) on the square domain with b
12
=
24.
(a)
u
h at T
=
1
(b)
v
h at T
=
1
(c)
u
h at T
=
10
(d)
v
h at T
=
10
(e)
u
h at T
=
25
(f)
v
h at T
=
25Figure 6
The time evolution of numerical approximates for fluxesq on the square domain with b
12
=
24.
(a)
x-direction component of q
h at T
=
1
(b)
y-direction component of q
h at T
=
1
(c)
x-direction component of q
h at T
=
10
(d)
y-direction component of q
h at T
=
10
(e)
x-direction component of q
h at T
=
25
(f)
y-direction component of q
h at T
=
25Figure 7
The time evolution of numerical solutionsu
h (left) and v
h (right) on the square domain with b
12
=
32.
(a)
u
h at T
=
1
(b)
v
h at T
=
1
(c)
u
h at T
=
10
(d)
v
h at T
=
10
(e)
u
h at T
=
25
(f)
v
h at T
=
25Figure 8
The time evolution of numerical approximates for fluxesq on the square domain with b
12
=
32.
(a)
x-direction component of q
h at T
=
1
(b)
y-direction component of q
h at T
=
1
(c)
x-direction component of q
h at T
=
10
(d)
y-direction component of q
h at T
=
10
(e)
x-direction component of q
h at T
=
25
(f)
y-direction component of q
h at T
=
25Simulations for these two sets of parameters have also been conducted in [23] with finite volume methods and in [26] with two kinds of finite element methods. We observe that the patterns obtained by our method agree well with those in [23, 26]. Meanwhile our method possesses its own advantages. By using the LDG method for spatial discretization, our method derives not only numerical solutions as the references do but also numerical approximations for fluxes. And it is easy to derive high-order spatial approximations which is difficult for finite volume schemes. Furthermore, our method reduces the computational cost greatly and CPU time for this simulation carried up to T
=
25 is 2632 s for the first set (b
12
=
24) and 2665 s for the second set (b
12
=
32). The reason is that, by employing the IIF scheme for temporal discretization, our method relaxes the strict time step restriction that is necessary for explicit schemes including the fourth-order exponential time differencing scheme used in [26] and allows us to compute element by element and avoids solving a global system of nonlinear algebraic equations as the backward difference advancing scheme applied in [23] does.Example 4.
We compute the Brusselator model (33)-(34) on the circular domain Ω
=
x
,
y
∣
x
2
+
y
2
≤
1
0
2. In the computation, the circular domain is divided into 4794 triangles and 2478 vertices, and now h
m
i
n
=
0.2405; mesh size h
=
0.4962 (see Figure 2(b)). The time step size is taken as (29) with c
f
l
=
0.03.From the numerical simulations in Example3, we have found that the patterns of numerical solutions u
h and v
h are always of the same type. Consequently, we can restrict our analysis of pattern formation to u
h.Figure9 exhibits graphs of numerical solutions u
h at T
=
1, T
=
10, and T
=
50 for b
12
=
24 (left) and b
12
=
32 (right), respectively, which also agree well with those in [23]. And the x-direction components of q
h at T
=
1, T
=
10, and T
=
50 for b
12
=
24 (left) and b
12
=
32 (right) are plotted in Figure 10. Compared with Example 3, we observe that the same patterns can be obtained on more complex geometry with larger final time.Figure 9
The time evolution ofu
h for b
12
=
24 (left) and b
12
=
32 (right) on the circular domain.
(a)
u
h at T
=
1 (b
12
=
24)
(b)
u
h at T
=
1 (b
12
=
32)
(c)
u
h at T
=
10 (b
12
=
24)
(d)
u
h at T
=
10 (b
12
=
32)
(e)
u
h at T
=
50 (b
12
=
24)
(f)
u
h at T
=
50 (b
12
=
32)Figure 10
Graphs ofx-direction component of q
h for b
12
=
24 (left) and b
12
=
32 (right) on the circular domain.
(a)
T
=
1 (b
12
=
24)
(b)
T
=
1 (b
12
=
32)
(c)
T
=
10 (b
12
=
24)
(d)
T
=
10 (b
12
=
32)
(e)
T
=
50 (b
12
=
24)
(f)
T
=
50 (b
12
=
32)And it is also worthy to point out that our method reduces the computational cost greatly and CPU time for this case carried up toT
=
50 is 1653 s for the first set (b
12
=
24) and 1604 s for the second set (b
12
=
32).


## 3.1. A Test Problem with Exact Solutions

To validate the spatial accuracy of our method numerically, we firstly consider a simple auxiliary reaction-diffusion system given in [23].Example 2.
The test reaction-diffusion system with cross-diffusion is of the following form:(31)
∂
u
∂
t
-
Δ
b
11
u
+
b
12
v
=
b
12
v
,
∂
v
∂
t
-
Δ
b
21
u
+
b
22
v
=
4
b
21
u
,defined on square Ω
=
[
0,2
π
]
2 and here b
11
=
b
22
=
1, b
12
=
1.5, and b
21
=
0.5. The exact solution is (32)
u
x
,
y
,
t
=
exp
⁡
-
4
t
cos
⁡
2
x
+
cos
⁡
2
y
,
v
x
,
y
,
t
=
exp
⁡
-
t
cos
⁡
x
+
cos
⁡
y
,from which we can derive the initial conditions.The simulation for Example2 is carried up to T
=
0.5 with c
f
l
=
0.05. The errors in L
2-norm and L
∞-norm are measured for both solutions and the fluxes on various triangular meshes. And the results are displayed in Table 1, where N
n denotes the number of triangular vertices. We can observe optimal convergence rates of O
(
h
2
) for solutions u
,
v in both L
2-norm and L
∞-norm and suboptimal convergence rates of O
(
h
) for fluxes q
,
p. These rates are consistent with theoretical results for the LDG method [35].Table 1
Numerical results for Example2 at T
=
0.5.
N
e
N
n
u
-
u
h
0
Order
u
-
u
h
∞
Order
v
-
v
h
0
Order
v
-
v
h
∞
Order
42
30
2.9159E − 1
—
1.0756E − 1
—
3.9750E − 1
—
1.1964E − 1
—
150
92
1.0714E − 1
1.5449
3.7922E − 2
1.6085
1.1495E − 1
1.9144
4.6267E − 2
1.4659
596
331
2.5577E − 2
2.0136
1.1788E − 2
1.6425
2.8564E − 2
1.9572
1.3543E − 2
1.7271
2392
1261
6.1548E − 3
2.2018
3.2021E − 3
2.0145
6.7407E − 3
2.2320
3.6577E − 3
2.0233
9224
4739
1.5429E − 3
2.0212
8.5707E − 4
1.9255
1.6793E − 3
2.0304
9.7541E − 4
1.9309
36838
18672
3.7359E − 4
1.9233
2.2556E − 4
1.8104
4.0320E − 4
1.9348
2.5540E − 4
1.8173
h
CPU (s)
q
-
q
h
0
Order
q
-
q
h
∞
Order
p
-
p
h
0
Order
p
-
p
h
∞
Order
1.8596
0.22
5.5622E − 1
—
1.6628E − 1
—
3.9209E − 1
—
1.2379E − 1
—
0.9727
0.28
1.8244E − 1
1.7200
5.4618E − 2
1.7178
1.9275E − 1
1.0957
7.6123E − 2
0.7502
0.4775
1.34
7.1312E − 2
1.3205
2.6911E − 2
0.9950
7.8795E − 2
1.2575
3.2105E − 2
1.2136
0.2501
9.75
3.1589E − 2
1.2586
1.3518E − 2
1.0643
3.5300E − 2
1.2412
1.5165E − 2
1.1593
0.1261
86.75
1.5204E − 2
1.0683
6.9372E − 3
0.9746
1.7065E − 2
1.0618
7.5730E − 3
1.0145
0.0603
473.50
7.4107E − 3
0.9746
3.3396E − 3
0.9914
8.2970E − 3
0.9780
3.8574E − 3
0.9148Moreover, Figure3 demonstrates graphs of numerical solutions u
h (left) and v
h (right) for this example, which are derived under the Delaunay partition shown in Figure 2(a). In addition, the numerical approximates for fluxes q
h (top) and p
h (bottom) are plotted in Figure 4.Figure 3
Graphs of numerical solutionsu
h and v
h at T
=
0.5 under the partition shown in Figure 2(a).
(a)
u
h
(b)
v
hFigure 4
Numerical approximates for fluxesq (top) and p (bottom) at T
=
0.5 under the partition shown in Figure 2(a).
(a)
x-direction component of q
h
(b)
y-direction component of q
h
(c)
x-direction component of p
h
(d)
y-direction component of p
h


## 3.2. Application to the Brusselator Model

In this subsection, we apply the LDG method coupled with the IIF scheme to solve the Brusselator model with cross-diffusion:(33)
∂
u
∂
t
-
Δ
b
11
u
+
b
12
v
=
-
β
+
1
u
+
u
2
v
+
α
,
∂
v
∂
t
-
Δ
b
21
u
+
b
22
v
=
β
u
-
u
2
v
,with the initial conditions chosen as small random perturbations of the equilibrium:(34)
u
x
,
0
=
5.8
+
1
3
rand
⁡
x
,
v
x
,
0
=
0.13
+
1
10
rand
⁡
x
,where x
=
(
x
,
y
)
∈
Ω, and rand: Ω
→
[
0,1
] is a random function in Fortran. In the following simulations, we take Ω as the square domain Ω
=
[
0,20
]
2 and circular domain Ω
=
{
(
x
,
y
)
∣
x
2
+
y
2
≤
1
0
2
}, respectively.Similar to [23], we set the parameters in the Brusselator model as(35)
b
11
=
0.4
,
b
22
=
2
,
b
21
=
0.02
,
α
=
6
,
β
=
1
,with b
12
=
24 and 32, respectively, based on which patterns are expected to appear. Actually, according to Theorem  2.3 in [23], the choice (35) is sufficient for the positive equilibrium point being linearly unstable and b
12 is taken as the Turing bifurcation parameter. Then the threshold b
12
=
22.2665 is obtained from (2.5) in [23]. Furthermore, it is observed that the increase of b
12 over the threshold with b
12
=
24 and b
12
=
32 yields the formation of spotted and labyrinthine patterns.Example 3.
We solve the Brusselator model (33)-(34) on the square domain Ω
=
[
0,20
]
2. In the computation, the square is divided into 16370 triangles and 8354 vertices, and now h
m
i
n
=
0.1492; mesh size h
=
0.2930. The time step size is taken as (29) with c
f
l
=
0.05.Different patterns will be obtained by selecting two sets of values for parametersb
12. The first set (b
12
=
24) leads to a spotted pattern as shown in Figure 5. It presents numerical approximations for solutions u (left) and v (right) at different values of final time T
=
1, T
=
10, and T
=
25, respectively. In addition, our method derives the approximations for fluxes q
=
b
11
∇
u and p
=
b
22
∇
v at the same time, and the graphs for q
h are plotted in Figure 6. The second set (b
12
=
32) generates a labyrinthine pattern (see Figure 7). The numerical approximations for fluxes q and p can also be obtained, and here we give only graphs of q
h in Figure 8. We observe a larger amplitude of the patterns (higher gradients) for both species.Figure 5
The time evolution of numerical solutionsu
h (left) and v
h (right) on the square domain with b
12
=
24.
(a)
u
h at T
=
1
(b)
v
h at T
=
1
(c)
u
h at T
=
10
(d)
v
h at T
=
10
(e)
u
h at T
=
25
(f)
v
h at T
=
25Figure 6
The time evolution of numerical approximates for fluxesq on the square domain with b
12
=
24.
(a)
x-direction component of q
h at T
=
1
(b)
y-direction component of q
h at T
=
1
(c)
x-direction component of q
h at T
=
10
(d)
y-direction component of q
h at T
=
10
(e)
x-direction component of q
h at T
=
25
(f)
y-direction component of q
h at T
=
25Figure 7
The time evolution of numerical solutionsu
h (left) and v
h (right) on the square domain with b
12
=
32.
(a)
u
h at T
=
1
(b)
v
h at T
=
1
(c)
u
h at T
=
10
(d)
v
h at T
=
10
(e)
u
h at T
=
25
(f)
v
h at T
=
25Figure 8
The time evolution of numerical approximates for fluxesq on the square domain with b
12
=
32.
(a)
x-direction component of q
h at T
=
1
(b)
y-direction component of q
h at T
=
1
(c)
x-direction component of q
h at T
=
10
(d)
y-direction component of q
h at T
=
10
(e)
x-direction component of q
h at T
=
25
(f)
y-direction component of q
h at T
=
25Simulations for these two sets of parameters have also been conducted in [23] with finite volume methods and in [26] with two kinds of finite element methods. We observe that the patterns obtained by our method agree well with those in [23, 26]. Meanwhile our method possesses its own advantages. By using the LDG method for spatial discretization, our method derives not only numerical solutions as the references do but also numerical approximations for fluxes. And it is easy to derive high-order spatial approximations which is difficult for finite volume schemes. Furthermore, our method reduces the computational cost greatly and CPU time for this simulation carried up to T
=
25 is 2632 s for the first set (b
12
=
24) and 2665 s for the second set (b
12
=
32). The reason is that, by employing the IIF scheme for temporal discretization, our method relaxes the strict time step restriction that is necessary for explicit schemes including the fourth-order exponential time differencing scheme used in [26] and allows us to compute element by element and avoids solving a global system of nonlinear algebraic equations as the backward difference advancing scheme applied in [23] does.Example 4.
We compute the Brusselator model (33)-(34) on the circular domain Ω
=
x
,
y
∣
x
2
+
y
2
≤
1
0
2. In the computation, the circular domain is divided into 4794 triangles and 2478 vertices, and now h
m
i
n
=
0.2405; mesh size h
=
0.4962 (see Figure 2(b)). The time step size is taken as (29) with c
f
l
=
0.03.From the numerical simulations in Example3, we have found that the patterns of numerical solutions u
h and v
h are always of the same type. Consequently, we can restrict our analysis of pattern formation to u
h.Figure9 exhibits graphs of numerical solutions u
h at T
=
1, T
=
10, and T
=
50 for b
12
=
24 (left) and b
12
=
32 (right), respectively, which also agree well with those in [23]. And the x-direction components of q
h at T
=
1, T
=
10, and T
=
50 for b
12
=
24 (left) and b
12
=
32 (right) are plotted in Figure 10. Compared with Example 3, we observe that the same patterns can be obtained on more complex geometry with larger final time.Figure 9
The time evolution ofu
h for b
12
=
24 (left) and b
12
=
32 (right) on the circular domain.
(a)
u
h at T
=
1 (b
12
=
24)
(b)
u
h at T
=
1 (b
12
=
32)
(c)
u
h at T
=
10 (b
12
=
24)
(d)
u
h at T
=
10 (b
12
=
32)
(e)
u
h at T
=
50 (b
12
=
24)
(f)
u
h at T
=
50 (b
12
=
32)Figure 10
Graphs ofx-direction component of q
h for b
12
=
24 (left) and b
12
=
32 (right) on the circular domain.
(a)
T
=
1 (b
12
=
24)
(b)
T
=
1 (b
12
=
32)
(c)
T
=
10 (b
12
=
24)
(d)
T
=
10 (b
12
=
32)
(e)
T
=
50 (b
12
=
24)
(f)
T
=
50 (b
12
=
32)And it is also worthy to point out that our method reduces the computational cost greatly and CPU time for this case carried up toT
=
50 is 1653 s for the first set (b
12
=
24) and 1604 s for the second set (b
12
=
32).


## 4. Conclusions

In this paper, we have developed the LDG method, coupled with the Krylov IIF time discretization, for solving reaction-cross-diffusion systems. By using LDG methods, we can obtain not only numerical solutions but also approximations for fluxes at the same time. Furthermore, another important property of LDG method is that the computation can proceed element by element and can also be remained, which benefits from applying the IIF method for temporal discretization. And it also relaxes the strict time step restriction that is necessary for explicit schemes.Experimental convergence rates were obtained for the LDG approximations by a test system with exact solutions, which shows optimal order for the solutions and suboptimal order for the fluxes in bothL
2-norm and L
∞-norm, just as the theoretical analysis in [35]. Then we conduct numerical simulations for the Brusselator system modelling an autocatalytic chemical reaction to verify the expected results in terms of behavior of the generated patterns.



---

*Source: 5345032-2016-06-29.xml*