The main purpose of this paper is to introduce the gradient discretisation method (GDM) to a system of reaction diffusion equations subject to non-homogeneous Dirichlet boundary conditions. Then, the authors show that the GDM provides a comprehensive convergence analysis of several numerical methods for the considered model. The convergence is established without non-physical regularity assumptions on the solutions.
In this paper, the authors use the GDM to discretise a system of reaction diffusion equations with non-homogeneous Dirichlet boundary conditions.
The authors provide a generic convergence analysis of a system of reaction diffusion equations. The authors introduce a specific example of numerical scheme that fits in the gradient discretisation method. The authors conduct a numerical test to measure the efficiency of the proposed method.
This work provides a unified convergence analysis of several numerical methods for a system of reaction diffusion equations. The generic convergence is proved under the classical assumptions on the solutions.
1. Introduction
In this paper, we study a system of reaction diffusion equations:
where μ1 and μ2 are the diffusion coefficients corresponding to the chemical concentrations and , respectively. The functions F and G describe the governing kinetics of the chemical reactions. A preprint has been posted on arXive [1]. Some biochemical phenomena have been expressed in the literature based on the choice of these reaction terms. For an example, in 2014, Yung–Rong Lee and Sy-Sang Lia [2] proposed a reaction–diffusion model that stimulated the relationship between concentrations of oxygen and lactic acid to simulate oxidative phosphorylation and glycolysis reactions in tissues. They concluded that a reaction–diffusion model can generate and maintain the ideal micro-environment for stem cells.
Another example is the Gray-Scott model, a very well-known reaction–diffusion system. The model describes the chemical reaction between two substances, an activator and substrate , where both of which diffuse over time and so represents what is called cubic autocatalysis [3].
The Brusselator model is also an example of such chemical reaction systems, which is used to describe mechanism of chemical reaction–diffusion with non-linear oscillations [4]. The model occurs in a large number of physical problems such as the formation of ozone by atomic oxygen, in enzymatic reactions, and arises in laser and plasma physics from multiple coupling between modes [5, 6].
The gradient discretisation method (GDM) is a generic framework to design numerical schemes together with their convergence analysis for different models, which are based on partial differential equations. It covers a variety of numerical schemes, such as finite volumes, finite elements, discontinuous Galerkin, etc. We refer the reader to Refs. [7–14] and the monograph [15] for a complete presentation. The main purpose of this paper is to introduce the GDM to a system of reaction–diffusion equations subject to non-homogeneous Dirichlet boundary conditions. Then, we show that the GDM provides a comprehensive convergence analysis of several numerical methods for the considered model. The convergence is established without non-physical regularity assumptions on the solutions since it is based on discrete compactness techniques detailed in Ref. [16].
The paper is organised as follows. Section 2 introduces the continuous model and its weak formulation. Section 3 describes the GDM for the model and states four required properties. Section 4 states the theorem corresponding to the convergence results. In Section 5, we include numerical test by employing a finite volume scheme, namely the hybrid mimetic mixed (HMM) method, to study and analyse the behaviour of the solutions of the Brusselator model as an example of the biochemical systems. The resultant relative errors with respect to the mesh size are also studied.
2. Continuous model
We consider the following biochemical system of partial differential equations:
Our analysis focuses on the weak formulation of the above reaction diffusion model. Let us assume the following properties on the data of the model.
The assumptions on the data in Problem (2.1)–(2.6) are the following:
Ω is an open bounded connected subset of , T > 0, ,
(uini, vini) are in L∞(Ω) × L∞(Ω),
g and h are traces of functions in L2(0, T; H1(Ω)) whose time derivatives are in L2(0, T; H−1(Ω)) and
the functions are Lipschitz continuous with Lipschitz constants LF and LG, respectively.
Under Assumptions 2.1, the weak solution of (2.1)–(2.6) is seeking
3. Discrete problem
The analysis of numerical schemes for the approximation of solutions to our model is performed using the GDM. The first step to reach this analysis is the reconstruction of a set of discrete spaces and operators, which is called gradient discretisation.
(Gradient discretisation for time-dependent problems with non-homogeneous Dirichlet boundary conditions). Let Ω be an open subset of (with d > 1) and T > 0. A gradient discretisation for time-dependent problems with non-homogeneous Dirichlet boundary conditions is ), where
the set of discrete unknowns is the direct sum of two finite dimensional spaces on , corresponding, respectively, to the interior unknowns and to the boundary unknowns,
the linear mapping is an interpolation operator for the trace,
the function reconstruction is a linear operator,
the gradient reconstruction is a linear operator and must be defined so that defines a norm on ,
is a linear and continuous interpolation operator for the initial conditions and
t(0) = 0 < t(1) < …. < t(N) = T are discrete times.
Let us introduce some notations to define the space–time reconstructions and and the discrete time derivative for .
For a.e x ∈ Ω, for all n ∈ {0, …, N − 1} and for all t ∈ (t(n), t(n+1)], let
Set and to define
Setting the gradient discretisation defined previously in the place of the continuous space and operators in the weak formulation of the model leads to a numerical scheme, called a gradient scheme.
(Gradient scheme). The gradient scheme for the continuous problem (2.7) is to find families of pair such that , and for all n = 0, …, N − 1, u(n+1) and v(n+1) satisfy the following equalities:
In order to establish the stability and convergence of the above gradient scheme, sequences of gradient discretisations described in Definition 3 are required to satisfy four properties: coercivity, consistency, limit–conformity and compactness.
(Coercivity). Let be a gradient discretisation and let be defined by
A sequence of a gradient discretisation is coercive if is bounded.
(Consistency). If is a gradient discretisation, define the function by, for φ ∈ H1(Ω),
A sequence of gradient discretisations is consistent if, as m → ∞
for all φ ∈ H1(Ω), ,
for all w ∈ L2(Ω), in L2(Ω) and
.
(Limit-conformity). If is a gradient discretisation and the space Hdiv = {ψ ∈ L2(Ω)d: divψ ∈ L2(Ω)}, define the function by, for all ψ ∈ Hdiv,
A sequence of gradient discretisations is limit-conforming if for all ψ ∈ Hdiv, , as m → ∞.
(Compactness). A sequence of gradient discretisation in the sense of Definition 3.1 is compact if for any sequence , such that is bounded, the sequence is relatively compact in L2(Ω).
(Dual norm on ). If be a gradient discretisation, the dual norm on is defined by, for all ,
4. Convergence results
Our convergence results are stated in the following theorem.
(Convergence of the gradient scheme). Assume (2.1) and let be a sequence of gradient discretisations, that is coercive, consistent, limit-conforming and compact. For , let be a solution to the gradient scheme (3.1) with . Then there exists a weak solution of (2.7) and a subsequence of gradient discretisations, still denoted by , such that, as m → ∞,
converges strongly to in L∞(0, T; L2(Ω)),
converges strongly to in L∞(0, T; L2(Ω)),
converges strongly to in L2(Ω × (0,T))d and
converges strongly to in L2(Ω × (0,T))d.
The proof relies on the compactness arguments as in Ref. [15] and is divided into four stages.
Step 1: Take liftings of g and of h such that and . Thanks to the density of space–time tensorial functions in the space L2(0, T; H1(Ω)) established in [[17], Corollary 1.3.1], we can express the liftings and in the following way: let , (ϕ)i = 1,…,ℓ, (ξ)i = 1,…,ℓ ⊂ C∞([0, T]) and such that
Let be defined by and for 0, …, Nm, where
From the consistency property, as m → ∞, we have.
strongly in L2(Ω × (0, T)) and strongly in L2(Ω × (0, T)),
strongly in L2(Ω × (0,T))d and strongly in L2(Ω × (0,T))d and
strongly in L2(Ω × (0, T)) and strongly in L2(Ω × (0, T)).
For any solution (u, v) to the gradient scheme (3.1), writing and , we have for all ,
Step 2: We need to have estimates on the quantities , , , and .
Apply the inequality, , , to the first terms in the above equalities, sum on n = 0, …, m − 1, for some m = 0, …, N and apply the Cauchy–Schwarz inequality to obtain
From the Lipschitz continuous assumptions on F and G, one has, with letting L = max(LF, LG) and C0 = max(|F(0)|, |G(0)|),
Then, using the Young's inequality in the right-hand side of the inequalities (with ɛi > 0, i = 1, …, 9), we conclude
where , , and C1 depends on , and , which are bounded.
Thanks to the Gronwall inequality [[18], Lemma 5.1] and to the coercivity property, the above inequalities can be written as
Since the terms on the right hand side can be simplified with terms on the left hand side in the both relations, we can combine the above inequalities together and take the supremum on m = 0, …, N to obtain the desired estimates.
Step 3: We need to established estimates on and . Take generic test functions φ and ψ in (4.1). Use the Cauchy–Schwarz inequality to get, thanks to assumptions (2.1) and to the coercivity properties,
The desired estimates is then obtained by taking the supremum over with , multiplying by δt(n+1), summing over n = 0, …, N − 1, thanks to the estimates obtained in the previous step.
Step 4: Owing to these estimates and the strong convergence of , , and , the remaining of the proof is then similar to that of [[15], Theorem 3.2 ]. □
5. Numerical results
To measure the efficiency of the gradient scheme (3.1) for the continuous problem (2.1)–(2.6), we consider a particular choice of the gradient discrtisation method known as the HMM method, which is a kind of finite volume scheme and can be written in three different formats; the hybrid finite volume method [19], the (mixed-hybrid) mimetic finite differences methods [20] and the mixed finite volume methods [21]. For the sake of completeness we briefly recall the definition of this gradient discretisation. Let be the polytopal mesh of the spatial domain Ω used in the previous section and described in [[15], Definition 7.2]. The elements of the gradient discretisation are as follows:
The discrete spaces are
The non-conforming piecewise affine reconstruction is defined by
The reconstructed gradients is a piecewise constant on the cells (broken gradient), defined by
where a cell-wise constant gradient ∇Kφ and a stabilisation term RK(φ) are, respectively, defined by
in which xσ is centre of mass of σ, xK is the gravity centre of cell K, dK,σ is the orthogonal distance between xK and , nK,σ is the unit vector normal to σ outward to K and DK,σ is the convex hull of σ ∪ {xK}.
The interpolant is defined by
the interpolant is defined by
The HMM scheme for (3.1) is the gradient scheme (2.7) written with the gradient discretisation constructed above.
As a test, we consider the Brusselator reaction–diffusion model (2.1)–(2.6) with non-homogeneous Dirichlet boundary conditions over the domain Ω = [0,1]2. The reaction functions in the Brusselator system are defined as follows:
where a is positive constant, and b is a parameter that can be varied to result in a range of different patterns. With x = (x, y) ∈ Ω, the exact solution in such a case is given as [5].
with parameters chosen as a = 0, b = 1, μ1 = μ2 = 0.25.
The initial and the Dirichlet boundary conditions are extracted from the analytical solutions (5.2). The simulation is performed on a sequence of triangular meshes and is done up to T = 1. The chosen meshes are of size h = 0.125, h = 0.0625, h = 0.03125 and h = 0.015625, respectively, with time step is fixed as 0.0001. Table 1 shows the relative errors on and and the corresponding rates of convergence with respect to the mesh size h. The resultant errors on the solutions and are proportional to the mesh size h, indicating that the HMM scheme behaves very well.
The relative errors and convergence rates w.r.t. the mesh size h at time t = 1 for the Brusselator model with parameters chosen as a = 0, b = 1, μ1 = μ2 = 0.25
| h | Rate | Rate | ||
|---|---|---|---|---|
| 0.125 | 0.000720746 | – | 0.000561639 | – |
| 0.0625 | 0.000184132 | 1.968753 | 0.000140295 | 2.0011797 |
| 0.03125 | 0.0000501972 | 1.8750586 | 0.0000342813 | 2.03296997 |
| 0.015625 | 0.0000149187 | 1.750485 | 0.00000688301 | 2.31630842 |
| h | Rate | Rate | ||
|---|---|---|---|---|
| 0.125 | 0.000720746 | – | 0.000561639 | – |
| 0.0625 | 0.000184132 | 1.968753 | 0.000140295 | 2.0011797 |
| 0.03125 | 0.0000501972 | 1.8750586 | 0.0000342813 | 2.03296997 |
| 0.015625 | 0.0000149187 | 1.750485 | 0.00000688301 | 2.31630842 |
Moreover, the L2 relative errors on the gradients of the solutions with respect to the mesh size h are shown in log-log scale Figure 1a for and in Figure 1b for . A line of slope one is added in both figures as a reference. We observe that the relative errors on and scale linearly with h, giving a rate of convergence of one, which are compatible with behaviour expectations associated with the low-order methods such as the HMM method.
The relative errors on (a) and (b) w.r.t. the mesh size h at time t = 1 for the Brusselator model with parameters chosen as in Table 1
The relative errors on (a) and (b) w.r.t. the mesh size h at time t = 1 for the Brusselator model with parameters chosen as in Table 1
6. Conclusion
We developed the GDM for a system of reaction–diffusion equations, including non-homogeneous Dirichlet boundary conditions. Without non-physical assumptions on the model data, we proved the existence of the weak solution for the continuous model and established the strong convergence for the discrete solution. We showed through a numerical test the efficiency of mixed finite volume methods.
The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work [Grant Code: 19-SCI-1-01-0027].

