The gradient discretisation method for the chemical reactions of biochemical systems

Abstract

Purpose

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.

Design/methodology/approach

In this paper, the authors use the GDM to discretise a system of reaction diffusion equations with non-homogeneous Dirichlet boundary conditions.

Findings

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.

Originality/value

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.

Citation

Alnashri, Y. and Alzubaidi, H. (2024), "The gradient discretisation method for the chemical reactions of biochemical systems", Arab Journal of Mathematical Sciences, Vol. 30 No. 1, pp. 67-80. https://doi.org/10.1108/AJMS-01-2022-0021

1. Introduction

In this paper, we study a system of reaction diffusion equations:

(1.1a)∂tu¯(x,t)−μ1div(∇u¯(x,t))=F(u¯,v¯),(x,t)∈Ω×(0,T),

(1.1b)∂tv¯(x,t)−μ2div(∇v¯(x,t))=G(u¯,v¯),(x,t)∈Ω×(0,T),

where μ₁ and μ₂ are the diffusion coefficients corresponding to the chemical concentrations u¯ and v¯, 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 v¯ and substrate u¯, 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:

(2.1)∂tu¯(x,t)−μ1div(∇u¯(x,t))=F(u¯,v¯),(x,t)∈Ω×(0,T),

(2.2)∂tv¯(x,t)−μ2div(∇v¯(x,t))=G(u¯,v¯),(x,t)∈Ω×(0,T),

(2.3)u¯(x,t)=g,(x,t)∈∂Ω×(0,T),

(2.4)u¯(x,t)=h,(x,t)∈∂Ω×(0,T),

(2.5)u¯(x,0)=uini(x),x∈Ω,

(2.6)v¯(x,0)=vini(x),x∈Ω.

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.

Assumptions 2.1.

The assumptions on the data in Problem (2.1)–(2.6) are the following:

1. Ω is an open bounded connected subset of Rd(d>1), T > 0, μ1,μ2∈R+,

2. (u_ini, v_ini) are in L^∞(Ω) × L^∞(Ω),

3. g and h are traces of functions in L²(0, T; H¹(Ω)) whose time derivatives are in L²(0, T; H⁻¹(Ω)) and

4. the functions F,G:R2→R are Lipschitz continuous with Lipschitz constants L_F and L_G, respectively.

Under Assumptions 2.1, the weak solution of (2.1)–(2.6) is seeking

u¯∈{w∈L2(0,T;H1(Ω)):γ(u¯)=g},v¯∈{w∈L2(0,T;H1(Ω)):γ(v¯)=h} such that u¯∈L2(0,T;H1(Ω))∩C([0,T];L2(Ω)),v¯∈L2(0,T;H1(Ω))∩C([0,T];L2(Ω)),∀φ∈L2(0,T;H01(Ω)),∂tφ∈L2(0,T;L2(Ω)),φ(⋅,T)=0,∀ψ∈L2(0,T;H01(Ω)),∂tψ∈L2(0,T;L2(Ω)),ψ(⋅,T)=0,

(2.7a)−∫0T∫Ωu¯(x,t)∂tφ(x,t)dxdt+μ1∫0T∫Ω∇u¯(x,t)⋅∇φ(x,t)dxdt,−∫Ωuini(x)φ(x,0)dx=∫0T∫ΩF(u¯,v¯)φ(x,t)dxdt,

(2.7b)−∫0T∫Ωv¯(x,t)∂tψ(x,t)dxdt+μ2∫0T∫Ω∇v¯(x,t)⋅∇ψ(x,t)dxdt−∫Ωvini(x)ψ(x,0)dx=∫0T∫ΩG(u¯,v¯)ψ(x,t)dxdt.

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.

Definition 3.1.

(Gradient discretisation for time-dependent problems with non-homogeneous Dirichlet boundary conditions). Let Ω be an open subset of Rd (with d > 1) and T > 0. A gradient discretisation for time-dependent problems with non-homogeneous Dirichlet boundary conditions is D=(XD,ID,∂,ΠD,∇D,JD,(t(n))n=0,…,N)), where

1. the set of discrete unknowns XD=XD,0⊕XD,∂Ω is the direct sum of two finite dimensional spaces on R, corresponding, respectively, to the interior unknowns and to the boundary unknowns,

2. the linear mapping ID,∂:H12(∂Ω)→XD,∂Ω is an interpolation operator for the trace,

3. the function reconstruction ΠD:XD→L2(Ω) is a linear operator,

4. the gradient reconstruction ∇D:XD→L2(Ω)d is a linear operator and must be defined so that ‖⋅‖D≔‖∇D⋅‖L2(Ω)d defines a norm on XD,0,

5. JD:L∞(Ω)→XD is a linear and continuous interpolation operator for the initial conditions and

6. t⁽⁰⁾ = 0 < t⁽¹⁾ < …. < t^(N) = T are discrete times.

Let us introduce some notations to define the space–time reconstructions ΠDφ:Ω×[0,T]→R and ∇Dφ:Ω×[0,T]→Rd and the discrete time derivative δDφ:(0,T)→L2(Ω) for φ=(φ(n))n=0,…,N∈XDN.

For a.e x ∈ Ω, for all n ∈ {0, …, N − 1} and for all t ∈ (t⁽ⁿ⁾, t⁽ⁿ⁺¹⁾], let

ΠDφ(x,0)=ΠDφ(0)(x),ΠDφ(x,t)=ΠDφ(n+1)(x),∇Dφ(x,t)=∇Dφ(n+1)(x).

Set δt(n+12)=t(n+1)−t(n) and δtD=maxn=0,…,N−1δt(n+12) to define

δDφ(t)=δDn+12φ≔ΠD(φ(n+1)−φ(n))δtn+12.

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.

Definition 3.2.

(Gradient scheme). The gradient scheme for the continuous problem (2.7) is to find families of pair (u(n),v(n))n=0,…,N∈(ID,∂g+XD,0N+1)×(ID,∂h+XD,0N+1) such that (u(0),v(0))=(JDuini,JDvini), and for all n = 0, …, N − 1, u⁽ⁿ⁺¹⁾ and v⁽ⁿ⁺¹⁾ satisfy the following equalities:

(3.1a) ∫ΩδDn+12u(x)ΠDφ(x)+μ1∫Ω∇Du(n+1)(x)⋅∇Dφ(x)dx=∫ΩF(ΠDu(n+1)(x),ΠDv(n+1)(x))ΠDφ(x)dx,∀φ∈XD,0, and 

(3.1b) ∫ΩδDn+12v(x)ΠDψ(x)+μ2∫Ω∇Dv(n+1)(x)⋅∇Dψ(x)dx=∫ΩG(ΠDu(n+1)(x),ΠDv(n+1)(x))ΠDψ(x)dx,∀ψ∈XD,0.

In order to establish the stability and convergence of the above gradient scheme, sequences of gradient discretisations (Dm)m∈N described in Definition 3 are required to satisfy four properties: coercivity, consistency, limit–conformity and compactness.

Definition 3.3.

(Coercivity). Let D be a gradient discretisation and let CD be defined by

(3.2) CD=maxφ∈XD,0−{0}‖ΠDφ‖L2(Ω)‖∇Dφ‖L2(Ω)d.

A sequence (Dm)m∈N of a gradient discretisation is coercive if CD is bounded.

Definition 3.4.

(Consistency). If D is a gradient discretisation, define the function SD:H1(Ω)→[0,+∞) by, for φ ∈ H¹(Ω),

(3.3) SD(φ)=min‖ΠDw−φ‖L2(Ω)+‖∇Dw−∇φ‖L2(Ω)d:w∈XD, such that w−ID,∂γφ∈XD,0,

A sequence (Dm)m∈N of gradient discretisations is consistent if, as m → ∞

1. for all φ ∈ H¹(Ω), SDm(φ)→0,

2. for all w ∈ L²(Ω), ΠDmJDmw→w in L²(Ω) and

3. δtDm→0.

Definition 3.5.

(Limit-conformity). If D is a gradient discretisation and the space H_div = {ψ ∈ L²(Ω)^d: divψ ∈ L²(Ω)}, define the function WD:Hdiv→[0,+∞) by, for all ψ ∈ H_div,

(3.4) WD(ψ)=supw∈XD,0\{0}∫Ω(∇Dw⋅ψ+ΠDwdiv(ψ))dx‖∇Dw‖L2(Ω)d.

A sequence (Dm)m∈N of gradient discretisations is limit-conforming if for all ψ ∈ H_div, WDm(ψ)→0, as m → ∞.

Definition 3.6.

(Compactness). A sequence of gradient discretisation (Dm)m∈N in the sense of Definition 3.1 is compact if for any sequence (φm)m∈N∈XDm,0, such that (‖∇Dmφm‖L2(Ω)d)m∈N is bounded, the sequence (ΠDmφm)m∈N is relatively compact in L²(Ω).

Definition 3.7.

(Dual norm on ΠD(XD)). If D be a gradient discretisation, the dual norm ‖⋅‖⋆,D on ΠD(XD,0) is defined by, for all w∈ΠD(XD,0),

(3.5) ‖w‖⋆,D=sup∫Ωw(x)ΠDφ(x)dx:φ∈XD,0,‖∇Dφ‖L2(Ω)d=1.

4. Convergence results

Our convergence results are stated in the following theorem.

Theorem 4.1.

(Convergence of the gradient scheme). Assume (2.1) and let (Dm)m∈N be a sequence of gradient discretisations, that is coercive, consistent, limit-conforming and compact. For m∈N, let (um,vm)∈(IDm,∂g+XDm,0N+1)×(IDm,∂h+XDm,0N+1) be a solution to the gradient scheme (3.1) with D=Dm. Then there exists a weak solution (u¯,v¯) of (2.7) and a subsequence of gradient discretisations, still denoted by (Dm)m∈N, such that, as m → ∞,

1. ΠDmum converges strongly to u¯ in L^∞(0, T; L²(Ω)),

2. ΠDmvm converges strongly to v¯ in L^∞(0, T; L²(Ω)),

3. ∇Dmum converges strongly to ∇u¯ in L²(Ω × (0,T))^d and

4. ∇Dmvm converges strongly to ∇v¯ in L²(Ω × (0,T))^d.

Proof.

The proof relies on the compactness arguments as in Ref. [15] and is divided into four stages.

• Step 1: Take liftings g¯∈L2(0,T;H1(Ω)) of g and h¯∈L2(0,T;H1(Ω)) of h such that γg¯=g and γh¯=h. Thanks to the density of space–time tensorial functions in the space L²(0, T; H¹(Ω)) established in [[17], Corollary 1.3.1], we can express the liftings g¯ and h¯ in the following way: let ℓ∈N, (ϕ)_{i = 1,…,ℓ}, (ξ)_{i = 1,…,ℓ} ⊂ C^∞([0, T]) and (pi)i=1,…,ℓ,(qi)i=1,…,ℓ⊂H1(Ω) such that

g¯(x,t)=∑i=1ℓϕi(t)pi(x) and a.e. x∈Ω and for all t∈(0,T),

h¯(x,t)=∑i=1ℓξi(t)qi(x) and a.e. x∈Ω and for all t∈(0,T).

Let gD,hD∈XDNm+1 be defined by gD(n)=ϕ(t(n))IDp and hD(n)=ξ(t(n))IDq for 0, …, N_m, where

IDp=arg minφ∈ID,∂g+XD,0‖ΠDφ−p‖L2(Ω)+‖∇Dφ−∇p‖L2(Ω)d,

IDq=arg minφ∈ID,∂h+XD,0‖ΠDφ−q‖L2(Ω)+‖∇Dφ−∇q‖L2(Ω)d.

From the consistency property, as m → ∞, we have.

1. ΠDmgDm→g¯ strongly in L²(Ω × (0, T)) and ΠDmhDm→h¯ strongly in L²(Ω × (0, T)),

2. ∇DmgDm→∇g¯ strongly in L²(Ω × (0,T))^d and ∇DmhDm→∇h¯ strongly in L²(Ω × (0,T))^d and

3. δDm(n+12)gDm→∂tg¯ strongly in L²(Ω × (0, T)) and δDm(n+12)hDm→∂th¯ strongly in L²(Ω × (0, T)).

For any solution (u, v) to the gradient scheme (3.1), writing w1=u−gD∈XD,0 and w2=v−hD∈XD,0, we have for all φ,ψ∈XD,0,

(4.1a)∫ΩδDn+12w1(x)ΠDφ(x)+μ1∫Ω∇Dw1(n+1)(x)⋅∇Dφ(x)dx=∫ΩF(ΠD(w1+gD),ΠD(w2+hD))ΠDφ(x)dx−∫ΩδDn+12gD(x)ΠDφ(x)−μ1∫Ω∇DgD(n+1)(x)⋅∇Dφ(x)dx,

(4.1b)∫ΩδDn+12w2(x)ΠDψ(x)+μ2∫Ω∇Dw2(n+1)(x)⋅∇Dψ(x)dx=∫ΩG(ΠD(w1+gD),ΠD(w2+hD))ΠDψ(x)dx−∫ΩδDn+12hD(x)ΠDψ(x)−μ2∫Ω∇DhD(n+1)(x)⋅∇Dψ(x)dx.

• Step 2: We need to have estimates on the quantities ‖ΠDw1(t)‖L∞(0,T;L2(Ω)), ‖ΠDw2(t)‖L∞(0,T;L2(Ω)), ‖∇Dw1‖L2(Ω×(0,T))d, and ‖∇Dw2‖L2(Ω×(0,T))d.

Let n ∈ {0, …, N − 1} and put φ≔δt(n+12)w1(n+1) in (4.1a) and ψ≔δt(n+12)w2(n+1) in (4.1b). We have

∫ΩΠDw1(n+1)(x)−ΠDw1(n)(x)ΠDw1(n+1)(x)dx+μ1δtn+12∫Ω∇Dw1(n+1)(x)⋅∇Dw1(n+1)(x)dx=δtn+12∫ΩF(ΠD(w1(n+1)(x)+gD),ΠD(w2(n+1)(x)+hD))ΠDw1(n+1)(x)dx−δtn+12∫ΩδDn+12gD(x)ΠDw1(n+1)(x)dx−μ1δtn+12∫Ω∇DgD(n+1)(x)⋅∇Dw1(n+1)(x)dx,and

∫ΩΠDw2(n+1)(x)−ΠDw2(n)(x)ΠDw2(n+1)(x)dx+μ2δtn+12∫Ω∇Dw2(n+1)(x)⋅∇Dw2(n+1)(x)dxdt=δtn+12∫ΩG(ΠD(w1(n+1)(x)+gD),ΠD(w2(n+1)(x)+hD))ΠDw2(n+1)(x)dx−δtn+12∫ΩδDn+12hD(x)ΠDw2(n+1)(x)dx−μ2δtn+12∫Ω∇DhD(n+1)(x)⋅∇Dw2(n+1)(x)dx.

Apply the inequality, a,b∈R, (a−b)a≥12(|a|2−|b|2), 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

12‖ΠDw1(m)‖L2(Ω)2−12‖ΠDw1(0)‖2+μ1∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d2≤∑n=0m−1δtn+12‖F(ΠD(w1(n+1)+gD(n+1)),ΠD(w2(n+1)+hD(n+1)))‖L2(Ω)‖ΠDw1(n+1)‖L2(Ω)+∑n=0m−1δtn+12‖ΠDw1(n+1)‖L2(Ω)‖δDn+12gD‖L2(Ω)+μ1∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d‖∇DgD(n)‖L2(Ω)d, and 

12‖ΠDw2(m)‖L2(Ω)2−12‖ΠDw2(0)‖2+μ2∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d2≤∑n=0m−1δtn+12‖G(ΠD(w1(n+1)+gD(n+1)),ΠD(w2(n+1)+hD(n+1)))‖L2(Ω)‖ΠDw2(n+1)‖L2(Ω)+∑n=0m−1‖ΠDw2(n+1)‖L2(Ω)‖δDn+12hD(n+1)‖L2(Ω)+μ2∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d‖∇DhD(n)‖L2(Ω)d.

From the Lipschitz continuous assumptions on F and G, one has, with letting L = max(L_F, L_G) and C₀ = max(|F(0)|, |G(0)|),

12‖ΠDw1(m)‖L2(Ω)2−12‖ΠDw1(0)‖2+μ1∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d2≤∑n=0m−1δtn+12L‖ΠDw1(n)‖L2(Ω)2+‖ΠDw1(n)‖L2(Ω)‖ΠDgD(n)‖L2(Ω)+‖ΠDw1(n)‖L2(Ω)‖ΠDw2(n)‖L2(Ω)+‖ΠDw1(n)‖L2(Ω)‖ΠDhD(n)‖L2(Ω)+C0‖ΠDw1(n)‖L2(Ω)+∑n=0m−1δtn+12‖ΠDw1(n+1)‖L2(Ω)‖δDn+12gD‖L2(Ω)+μ1∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d‖∇DgD(n)‖L2(Ω)d, and 

12‖ΠDw2(m)‖L2(Ω)2−12‖ΠDw2(0)‖2+μ2∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d2≤∑n=0m−1δtn+12L‖ΠDw2(n)‖L2(Ω)2+‖ΠDw2(n)‖L2(Ω)‖ΠDgD(n)‖L2(Ω)+‖ΠDw2(n)‖L2(Ω)‖ΠDw1(n)‖L2(Ω)+‖ΠDw2(n)‖L2(Ω)‖ΠDhD(n)‖L2(Ω)+C0‖ΠDw2(n)‖L2(Ω)+∑n=0m−1δtn+12‖ΠDw2(n+1)‖L2(Ω)‖δDn+12hD‖L2(Ω)+μ2∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d‖∇DhD(n)‖L2(Ω)d.

Then, using the Young's inequality in the right-hand side of the inequalities (with ɛ_i > 0, i = 1, …, 9), we conclude

12‖ΠDw1(m)‖L2(Ω)2−12‖ΠDw1(0)‖2+μ12∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d2≤M1∑n=0m−1δtn+12‖ΠDw1(n)‖L2(Ω)2+12ε2∑n=0m−1δtn+12‖ΠDw2(n)‖L2(Ω)2+T2ε4C02+TM1C1,and

12‖ΠDw2(m)‖L2(Ω)2−12‖ΠDw2(0)‖2+μ22∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d2≤M2∑n=0m−1δtn+12‖ΠDw2(n)‖L2(Ω)2+12ε7∑n=0m−1δtn+12‖ΠDw1(n)‖L2(Ω)2+T2ε9C02+TM2C1,

where M1≔L+∑i=15εi2, M2≔L+∑i=610εi2, and C₁ depends on CD, ‖∇DgD‖L2(Ω)d and ‖∇DhD‖L2(Ω)d, which are bounded.

Thanks to the Gronwall inequality [[18], Lemma 5.1] and to the coercivity property, the above inequalities can be written as

12‖ΠDw1(m)‖L2(Ω)2+μ12∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)d2≤expM1CD2ε2∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)2+T2ε4C02+TM1C1+12‖ΠDw1(0)‖2,and

12‖ΠDw2(m)‖L2(Ω)2+μ22∑n=0m−1δtn+12‖∇Dw2(n)‖L2(Ω)d2≤expM2CD2ε7∑n=0m−1δtn+12‖∇Dw1(n)‖L2(Ω)2+T2ε9C02+TM2C1+12‖ΠDw2(0)‖2.

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 ‖w1‖⋆,D and ‖w2‖⋆,D. 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,

∫ΩδDn+12w1(x)ΠDφ(x)dx≤CD‖∇Dφ‖L2(Ω)dL‖ΠDw1(n+1)‖L2(Ω×(0,T))+‖ΠDw2(n+1)‖L2(Ω×(0,T))+‖ΠDgD(n+1)‖L2(Ω×(0,T))+‖ΠDhD(n+1)‖L2(Ω×(0,T))+C0+‖δDn+12gD(n+1)‖L2(Ω×(0,T))+μ1‖∇Dw1(n+1)‖L2(Ω×(0,T))d,

∫ΩδDn+12w2(x)ΠDψ(x)dx≤CD‖∇Dψ‖L2(Ω)dL‖ΠDw1(n+1)‖L2(Ω×(0,T))+‖ΠDw2(n+1)‖L2(Ω×(0,T))+‖ΠDgD(n+1)‖L2(Ω×(0,T))+‖ΠDhD(n+1)‖L2(Ω×(0,T))+C0+‖δDn+12hD(n+1)‖L2(Ω×(0,T))+μ1‖∇Dw2(n+1)‖L2(Ω×(0,T))d.

The desired estimates is then obtained by taking the supremum over φ,ψ∈XD,0 with ‖∇Dφ‖L2(Ω)d=‖∇Dψ‖L2(Ω)d=1, multiplying by δt⁽ⁿ⁺¹⁾, 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 ΠDmgDm, ΠDmhDm, ∇DmgDm and ∇DmhDm, 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 T=(M,F,P,V) 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:

1. The discrete spaces are

XD,0={v=((φK)K∈M,(φσ)σ∈F):φK,φσ∈R,φσ=0,∀σ∈F∩∂Ω},

XD,∂Ω=v=((vK)K∈M,(vσ)σ∈F):vK∈R,vσ∈R,vK=0 for all K∈M,vσ=0 for all σ∈Fint.

1. The non-conforming piecewise affine reconstruction ΠD is defined by

∀φ∈XD,∀K∈M, for a.e. x∈K,ΠDφ=φK on K.

1. The reconstructed gradients is a piecewise constant on the cells (broken gradient), defined by

∀φ∈XD,∀K∈M,∀σ∈FK,∇Dφ=∇Kφ+ddK,σRK(φ)nK,σ on DK,σ,

where a cell-wise constant gradient ∇_Kφ and a stabilisation term R_K(φ) are, respectively, defined by

∇Kφ=1|K|∑σ∈FK|σ|φσnK,σ and RK(φ)=(φσ−φK−∇Kφ⋅(xσ−xK))σ∈FK,

in which x_σ is centre of mass of σ, x_K is the gravity centre of cell K, d_K,σ is the orthogonal distance between x_K and σ∈FK, n_K,σ is the unit vector normal to σ outward to K and D_K,σ is the convex hull of σ ∪ {x_K}.

1. The interpolant JD:L∞(Ω)→XD is defined by

∀w∈L∞(Ω):JDw=((wK)K∈M,(wσ)σ∈F),∀K∈M,wK=1|K|∫Kw(x)dx and ∀σ∈F,wσ=0.

1. the interpolant ID,∂:H12(∂Ω)→XD,∂Ω is defined by

∀g∈H12(∂Ω):(ID,∂g)σ=1|σ|∫σg(x)ds(x), for all σ∈Fext such that σ⊂∂Ω.

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]². The reaction functions in the Brusselator system are defined as follows:

F(u¯,v¯)=a−(b+1)u¯+u¯2v¯,G(u¯,v¯)=bu¯−u¯2v¯,

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].

(5.2a)u¯(x,t)=exp(−x−y−0.5t),

(5.2b)v¯(x,t)=exp(x+y+0.5t),

with parameters chosen as a = 0, b = 1, μ₁ = μ₂ = 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 u¯ and v¯ and the corresponding rates of convergence with respect to the mesh size h. The resultant errors on the solutions u¯ and v¯ are proportional to the mesh size h, indicating that the HMM scheme behaves very well.

Moreover, the L₂ relative errors on the gradients of the solutions with respect to the mesh size h are shown in log-log scale Figure 1a for ∇u¯ and in Figure 1b for ∇v¯. A line of slope one is added in both figures as a reference. We observe that the relative errors on ∇u¯ and ∇v¯ 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.

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.

Figures

Figure 1

The relative errors on (a) ∇u¯ and (b) ∇v¯ w.r.t. the mesh size h at time t = 1 for the Brusselator model with parameters chosen as in Table 1

Table 1

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, μ₁ = μ₂ = 0.25

h	‖u¯(⋅,t(n))−ΠDu(n)‖L2(Ω)‖u¯(⋅,t(n))‖L2(Ω)	Rate	‖v¯(⋅,t(n))−ΠDv(n)‖L2(Ω)‖v¯(⋅,t(n))‖L2(Ω)	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