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)${\partial }_t\overline{u}\left(\mathit{x},t\right)-{\mu }_1\mathrm{d}\mathrm{i}\mathrm{v}(\nabla \overline{u}\left(\mathit{x},t\right))=F(\overline{u},\overline{v}),\left(\mathit{x},t\right)\in \mathrm{\Omega }\times \left(0,T\right),$

(1.1b)${\partial }_t\overline{v}\left(\mathit{x},t\right)-{\mu }_2\mathrm{d}\mathrm{i}\mathrm{v}(\nabla \overline{v}\left(\mathit{x},t\right))=G(\overline{u},\overline{v}),\left(\mathit{x},t\right)\in \mathrm{\Omega }\times \left(0,T\right),$

where μ₁ and μ₂ are the diffusion coefficients corresponding to the chemical concentrations $\overline{u}$ and $\overline{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 $\overline{v}$ and substrate $\overline{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)${\partial }_t\overline{u}\left(\mathit{x},t\right)-{\mu }_1\mathrm{d}\mathrm{i}\mathrm{v}(\nabla \overline{u}\left(\mathit{x},t\right))=F(\overline{u},\overline{v}),\left(\mathit{x},t\right)\in \mathrm{\Omega }\times \left(0,T\right),$

(2.2)${\partial }_t\overline{v}\left(\mathit{x},t\right)-{\mu }_2\mathrm{d}\mathrm{i}\mathrm{v}(\nabla \overline{v}(\mathit{x},t))=G(\overline{u},\overline{v}),\left(\mathit{x},t\right)\in \mathrm{\Omega }\times \left(0,T\right),$

(2.3)$\overline{u}\left(\mathit{x},t\right)=g,\left(\mathit{x},t\right)\in \partial \mathrm{\Omega }\times \left(0,T\right),$

(2.4)$\overline{u}\left(\mathit{x},t\right)=h,\left(\mathit{x},t\right)\in \partial \mathrm{\Omega }\times \left(0,T\right),$

(2.5)$\overline{u}\left(\mathit{x},0\right)=u_{\text{ini}}\left(\mathit{x}\right),\mathit{x}\in \mathrm{\Omega },$

(2.6)$\overline{v}\left(\mathit{x},0\right)=v_{\text{ini}}\left(\mathit{x}\right),\mathit{x}\in \mathrm{\Omega }.$

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 ${\mathbb{R}}^d\left(d>1\right)$, T > 0, ${\mu }_1,{\mu }_2\in {\mathbb{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:{\mathbb{R}}^2\to \mathbb{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

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

(2.7a)$\begin{array}{cc}\hfill -{\int }_0^T{\int }_{\mathrm{\Omega }}\overline{u}\left(\mathit{x},t\right){\partial }_t\phi \left(\mathit{x},t\right)\mathrm{d}\mathit{x}\mathrm{d}t& +{\mu }_1{\int }_0^T{\int }_{\mathrm{\Omega }}\nabla \overline{u}\left(\mathit{x},t\right)\cdot \nabla \phi \left(\mathit{x},t\right)\mathrm{d}\mathit{x}\mathrm{d}t,\hfill \\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

(2.7b)$\begin{array}{cc}\hfill -{\int }_0^T{\int }_{\mathrm{\Omega }}\overline{v}\left(\mathit{x},t\right){\partial }_t\psi \left(\mathit{x},t\right)\mathrm{d}\mathit{x}\mathrm{d}t& +{\mu }_2{\int }_0^T{\int }_{\mathrm{\Omega }}\nabla \overline{v}\left(\mathit{x},t\right)\cdot \nabla \psi \left(\mathit{x},t\right)\mathrm{d}\mathit{x}\mathrm{d}t\hfill \\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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 ${\mathbb{R}}^d$ (with d > 1) and T > 0. A gradient discretisation for time-dependent problems with non-homogeneous Dirichlet boundary conditions is $\mathcal{D}=\left(X_{\mathcal{D}},{\mathcal{I}}_{\mathcal{D},\partial },{\mathrm{\Pi }}_{\mathcal{D}},{\nabla }_{\mathcal{D}},J_{\mathcal{D}},{\left(t^{\left(n\right)}\right)}_{n=0,\dots ,N}\right)$), where

1. the set of discrete unknowns $X_{\mathcal{D}}=X_{\mathcal{D},0}\oplus X_{\mathcal{D},\partial \mathrm{\Omega }}$ is the direct sum of two finite dimensional spaces on $\mathbb{R}$, corresponding, respectively, to the interior unknowns and to the boundary unknowns,

2. the linear mapping ${\mathcal{I}}_{\mathcal{D},\partial }:H^{\frac{1}{2}}\left(\partial \mathrm{\Omega }\right)\to X_{\mathcal{D},\partial \mathrm{\Omega }}$ is an interpolation operator for the trace,

3. the function reconstruction ${\mathrm{\Pi }}_{\mathcal{D}}:X_{\mathcal{D}}\to L^2\left(\mathrm{\Omega }\right)$ is a linear operator,

4. the gradient reconstruction ${\nabla }_{\mathcal{D}}:X_{\mathcal{D}}\to L^2{\left(\mathrm{\Omega }\right)}^d$ is a linear operator and must be defined so that $\Vert \cdot {\Vert }_{\mathcal{D}}≔\Vert {\nabla }_{\mathcal{D}}\cdot {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}$ defines a norm on $X_{\mathcal{D},0}$,

5. $J_{\mathcal{D}}:L^{\infty }\left(\mathrm{\Omega }\right)\to X_{\mathcal{D}}$ 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 ${\mathrm{\Pi }}_{\mathcal{D}}\phi :\mathrm{\Omega }\times \left[0,T\right]\to \mathbb{R}$ and ${\nabla }_{\mathcal{D}}\phi :\mathrm{\Omega }\times \left[0,T\right]\to {\mathbb{R}}^d$ and the discrete time derivative ${\delta }_{\mathcal{D}}\phi :\left(0,T\right)\to L^2\left(\mathrm{\Omega }\right)$ for $\phi ={\left({\phi }^{\left(n\right)}\right)}_{n=0,\dots ,N}\in X_{\mathcal{D}}^N$.

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

$\begin{array}{c}{\mathrm{\Pi }}_{\mathcal{D}}\phi \left(\mathit{x},0\right)={\mathrm{\Pi }}_{\mathcal{D}}{\phi }^{\left(0\right)}\left(\mathit{x}\right),{\mathrm{\Pi }}_{\mathcal{D}}\phi \left(\mathit{x},t\right)={\mathrm{\Pi }}_{\mathcal{D}}{\phi }^{\left(n+1\right)}\left(\mathit{x}\right),\hfill \\ {\nabla }_{\mathcal{D}}\phi \left(\mathit{x},t\right)={\nabla }_{\mathcal{D}}{\phi }^{\left(n+1\right)}\left(\mathit{x}\right).\hfill \end{array}$

Set $\delta t^{\left(n+\frac{1}{2}\right)}=t^{\left(n+1\right)}-t^{\left(n\right)}$ and $\delta t_{\mathcal{D}}={\max }_{n=0,\dots ,N-1}\delta t^{\left(n+\frac{1}{2}\right)}$ to define

${\delta }_{\mathcal{D}}\phi \left(t\right)={\delta }_{\mathcal{D}}^{\left(n+\frac{1}{2}\right)}\phi ≔\frac{{\mathrm{\Pi }}_{\mathcal{D}}\left({\phi }^{\left(n+1\right)}-{\phi }^{\left(n\right)}\right)}{\delta t^{\left(n+\frac{1}{2}\right)}}.$

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 ${\left(u^{\left(n\right)},v^{\left(n\right)}\right)}_{n=0,\dots ,N}\in \left({\mathcal{I}}_{\mathcal{D},\partial }g+X_{\mathcal{D},0}^{N+1}\right)\times \left({\mathcal{I}}_{\mathcal{D},\partial }h+X_{\mathcal{D},0}^{N+1}\right)$ such that $\left(u^{\left(0\right)},v^{\left(0\right)}\right)=\left(J_{\mathcal{D}}{u}_{\text{ini}},J_{\mathcal{D}}{v}_{\text{ini}}\right)$, and for all n = 0, …, N − 1, u⁽ⁿ⁺¹⁾ and v⁽ⁿ⁺¹⁾ satisfy the following equalities:

(3.1a) $\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

(3.1b) $\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

In order to establish the stability and convergence of the above gradient scheme, sequences of gradient discretisations ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ described in Definition 3 are required to satisfy four properties: coercivity, consistency, limit–conformity and compactness.

Definition 3.3.

(Coercivity). Let $\mathcal{D}$ be a gradient discretisation and let $C_{\mathcal{D}}$ be defined by

(3.2) $C_{\mathcal{D}}=\underset{\phi \in X_{\mathcal{D},0}-\left\{0\right\}}{\max }\frac{\Vert {\mathrm{\Pi }}_{\mathcal{D}}\phi {\Vert }_{L^2\left(\mathrm{\Omega }\right)}}{\Vert {\nabla }_{\mathcal{D}}\phi {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}}.$

A sequence ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ of a gradient discretisation is coercive if $C_{\mathcal{D}}$ is bounded.

Definition 3.4.

(Consistency). If $\mathcal{D}$ is a gradient discretisation, define the function $S_{\mathcal{D}}:H^1\left(\mathrm{\Omega }\right)\to \left[0,+\infty \right)$ by, for φ ∈ H¹(Ω),

(3.3) $\begin{array}{c}\hfill S_{\mathcal{D}}\left(\phi \right)=\min \left\{\Vert {\mathrm{\Pi }}_{\mathcal{D}}w-\phi {\Vert }_{L^2\left(\mathrm{\Omega }\right)}+\Vert {\nabla }_{\mathcal{D}}w-\nabla \phi {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}:w\in X_{\mathcal{D}},\right.\\ \hfill \left.\text{\hspace{0.17em}such\hspace{0.17em}that\hspace{0.17em}}w-{\mathcal{I}}_{\mathcal{D},\partial }\gamma \phi \in X_{\mathcal{D},0}\right\},\end{array}$

A sequence ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ of gradient discretisations is consistent if, as m → ∞

1. for all φ ∈ H¹(Ω), $S_{{\mathcal{D}}_m}\left(\phi \right)\to 0$,

2. for all w ∈ L²(Ω), ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{J}_{{\mathcal{D}}_m}w\to w$ in L²(Ω) and

3. $\delta t_{{\mathcal{D}}_m}\to 0$.

Definition 3.5.

(Limit-conformity). If $\mathcal{D}$ is a gradient discretisation and the space H_div = {ψ ∈ L²(Ω)^d: divψ ∈ L²(Ω)}, define the function $W_{\mathcal{D}}:H_{\text{div}}\to \left[0,+\infty \right)$ by, for all ψ ∈ H_div,

(3.4) $W_{\mathcal{D}}\left(\psi \right)=\underset{w\in X_{\mathcal{D},0}\backslash \left\{0\right\}}{\sup }\frac{\left|{\int }_{\mathrm{\Omega }}\left({\nabla }_{\mathcal{D}}w\cdot \psi +{\mathrm{\Pi }}_{\mathcal{D}}w\mathrm{d}\mathrm{i}\mathrm{v}\left(\psi \right)\right)\mathrm{d}x\right|}{\Vert {\nabla }_{\mathcal{D}}w{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}}.$

A sequence ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ of gradient discretisations is limit-conforming if for all ψ ∈ H_div, $W_{{\mathcal{D}}_m}\left(\psi \right)\to 0$, as m → ∞.

Definition 3.6.

(Compactness). A sequence of gradient discretisation ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ in the sense of Definition 3.1 is compact if for any sequence ${\left({\phi }_m\right)}_{m\in \mathbb{N}}\in X_{{\mathcal{D}}_m,0}$, such that ${\left(\Vert {\nabla }_{{\mathcal{D}}_m}{\phi }_m{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}\right)}_{m\in \mathbb{N}}$ is bounded, the sequence ${\left({\mathrm{\Pi }}_{{\mathcal{D}}_m}{\phi }_m\right)}_{m\in \mathbb{N}}$ is relatively compact in L²(Ω).

Definition 3.7.

(Dual norm on ${\mathrm{\Pi }}_{\mathcal{D}}\left(X_{\mathcal{D}}\right)$). If $\mathcal{D}$ be a gradient discretisation, the dual norm $\Vert \cdot {\Vert }_{\star ,\mathcal{D}}$ on ${\mathrm{\Pi }}_{\mathcal{D}}\left(X_{\mathcal{D},0}\right)$ is defined by, for all $w\in {\mathrm{\Pi }}_{\mathcal{D}}\left(X_{\mathcal{D},0}\right)$,

(3.5) $\Vert w{\Vert }_{\star ,\mathcal{D}}=\sup \left\{{\int }_{\mathrm{\Omega }}w\left(\mathit{x}\right){\mathrm{\Pi }}_{\mathcal{D}}\phi \left(\mathit{x}\right)\mathrm{d}\mathit{x}:\phi \in X_{\mathcal{D},0},\Vert {\nabla }_{\mathcal{D}}\phi {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}=1\right\}.$

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 ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$ be a sequence of gradient discretisations, that is coercive, consistent, limit-conforming and compact. For $m\in \mathbb{N}$, let $\left(u_m,v_m\right)\in \left({\mathcal{I}}_{{\mathcal{D}}_m,\partial }g+X_{{\mathcal{D}}_m,0}^{N+1}\right)\times \left({\mathcal{I}}_{{\mathcal{D}}_m,\partial }h+X_{{\mathcal{D}}_m,0}^{N+1}\right)$ be a solution to the gradient scheme (3.1) with $\mathcal{D}={\mathcal{D}}_m$. Then there exists a weak solution $\left(\overline{u},\overline{v}\right)$ of (2.7)and a subsequence of gradient discretisations, still denoted by ${\left({\mathcal{D}}_m\right)}_{m\in \mathbb{N}}$, such that, as m → ∞,

1. ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{u}_m$ converges strongly to $\overline{u}$ in L^∞(0, T; L²(Ω)),

2. ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{v}_m$ converges strongly to $\overline{v}$ in L^∞(0, T; L²(Ω)),

3. ${\nabla }_{{\mathcal{D}}_m}{u}_m$ converges strongly to $\nabla \overline{u}$ in L²(Ω × (0,T))^d and

4. ${\nabla }_{{\mathcal{D}}_m}{v}_m$ converges strongly to $\nabla \overline{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 $\overline{g}\in L^2\left(0,T;H^1\left(\mathrm{\Omega }\right)\right)$ of g and $\overline{h}\in L^2\left(0,T;H^1\left(\mathrm{\Omega }\right)\right)$ of h such that $\gamma \overline{g}=g$ and $\gamma \overline{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 $\overline{g}$ and $\overline{h}$ in the following way: let $\ell \in \mathbb{N}$, (ϕ)_{i = 1,…,ℓ}, (ξ)_{i = 1,…,ℓ} ⊂ C^∞([0, T]) and ${\left(p_i\right)}_{i=1,\dots ,\ell },{\left(q_i\right)}_{i=1,\dots ,\ell }\subset H^1\left(\mathrm{\Omega }\right)$ such that

$\overline{g}\left(x,t\right)=\sum _{i=1}^{\ell }{\varphi }_i\left(t\right)p_i\left(\mathit{x}\right)\text{\hspace{0.17em}and\hspace{0.17em}a.e.\hspace{0.17em}}\mathit{x}\in \mathrm{\Omega }\text{\hspace{0.17em}and\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}}t\in \left(0,T\right),$

$\overline{h}\left(x,t\right)=\sum _{i=1}^{\ell }{\xi }_i\left(t\right)q_i\left(\mathit{x}\right)\text{\hspace{0.17em}and\hspace{0.17em}a.e.\hspace{0.17em}}\mathit{x}\in \mathrm{\Omega }\text{\hspace{0.17em}and\hspace{0.17em}for\hspace{0.17em}all\hspace{0.17em}}t\in \left(0,T\right).$

Let $g_{\mathcal{D}},h_{\mathcal{D}}\in X_{\mathcal{D}}^{N_m+1}$ be defined by $g_{\mathcal{D}}^{\left(n\right)}=\varphi \left(t^{\left(n\right)}\right)I_{\mathcal{D}}p$ and $h_{\mathcal{D}}^{\left(n\right)}=\xi \left(t^{\left(n\right)}\right)I_{\mathcal{D}}q$ for 0, …, N_m, where

$I_{\mathcal{D}}p=\underset{\phi \in {\mathcal{I}}_{\mathcal{D},\partial g}+X_{\mathcal{D},0}}{\mathrm{arg\; min}}\left(\Vert {\mathrm{\Pi }}_{\mathcal{D}}\phi -p{\Vert }_{L^2\left(\mathrm{\Omega }\right)}+\Vert {\nabla }_{\mathcal{D}}\phi -\nabla p{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}\right),$

$I_{\mathcal{D}}q=\underset{\phi \in {\mathcal{I}}_{\mathcal{D},\partial h}+X_{\mathcal{D},0}}{\mathrm{arg\; min}}\left(\Vert {\mathrm{\Pi }}_{\mathcal{D}}\phi -q{\Vert }_{L^2\left(\mathrm{\Omega }\right)}+\Vert {\nabla }_{\mathcal{D}}\phi -\nabla q{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}\right).$

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

1. ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{g}_{{\mathcal{D}}_m}\to \overline{g}$ strongly in L²(Ω × (0, T)) and ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{h}_{{\mathcal{D}}_m}\to \overline{h}$ strongly in L²(Ω × (0, T)),

2. ${\nabla }_{{\mathcal{D}}_m}{g}_{{\mathcal{D}}_m}\to \nabla \overline{g}$ strongly in L²(Ω × (0,T))^d and ${\nabla }_{{\mathcal{D}}_m}{h}_{{\mathcal{D}}_m}\to \nabla \overline{h}$ strongly in L²(Ω × (0,T))^d and

3. ${\delta }_{{\mathcal{D}}_m}^{\left(n+\frac{1}{2}\right)}{g}_{{\mathcal{D}}_m}\to {\partial }_t\overline{g}$ strongly in L²(Ω × (0, T)) and ${\delta }_{{\mathcal{D}}_m}^{\left(n+\frac{1}{2}\right)}{h}_{{\mathcal{D}}_m}\to {\partial }_t\overline{h}$ strongly in L²(Ω × (0, T)).

For any solution (u, v) to the gradient scheme (3.1), writing $w_1=u-g_{\mathcal{D}}\in X_{\mathcal{D},0}$ and $w_2=v-h_{\mathcal{D}}\in X_{\mathcal{D},0}$, we have for all $\phi ,\psi \in X_{\mathcal{D},0}$,

(4.1a)$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

(4.1b)$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

• Step 2: We need to have estimates on the quantities $\Vert {\mathrm{\Pi }}_{\mathcal{D}}{w}_1\left(t\right){\Vert }_{L^{\infty }\left(0,T;L^2\left(\mathrm{\Omega }\right)\right)}$, $\Vert {\mathrm{\Pi }}_{\mathcal{D}}{w}_2\left(t\right){\Vert }_{L^{\infty }\left(0,T;L^2\left(\mathrm{\Omega }\right)\right)}$, $\Vert {\nabla }_{\mathcal{D}}{w}_1{\Vert }_{L^2{\left(\mathrm{\Omega }\times \left(0,T\right)\right)}^d}$, and ${\Vert {\nabla }_{\mathcal{D}}{w}_2\Vert }_{L^2{\left(\mathrm{\Omega }\times \left(0,T\right)\right)}^d}$.

Let n ∈ {0, …, N − 1} and put $\phi ≔\delta t^{\left(n+\frac{1}{2}\right)}{w}_1^{\left(n+1\right)}$ in (4.1a) and $\psi ≔\delta t^{\left(n+\frac{1}{2}\right)}{w}_2^{\left(n+1\right)}$ in (4.1b). We have

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

Apply the inequality, $a,b\in \mathbb{R}$, $\left(a-b\right)a\ge \frac{1}{2}\left(|a{|}^2-|b{|}^2\right)$, 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

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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)|),

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

where $M_1≔L+{\sum }_{i=1}^5\frac{{\epsilon }_i}{2}$, $M_2≔L+{\sum }_{i=6}^{10}\frac{{\epsilon }_i}{2}$, and C₁ depends on $C_{\mathcal{D}}$, $\Vert {\nabla }_{\mathcal{D}}{g}_{\mathcal{D}}{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}$ and $\Vert {\nabla }_{\mathcal{D}}{h}_{\mathcal{D}}{\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^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

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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 $\Vert w_1{\Vert }_{\star ,\mathcal{D}}$ and $\Vert w_2{\Vert }_{\star ,\mathcal{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,

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill & +\left.\Vert {\mathrm{\Pi }}_{\mathcal{D}}{g}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}+\Vert {\mathrm{\Pi }}_{\mathcal{D}}{h}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}\right)+C_0\hfill \\ \hfill & +\left.\Vert {\delta }_{\mathcal{D}}^{\left(n+\frac{1}{2}\right)}{g}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}+{\mu }_1\Vert {\nabla }_{\mathcal{D}}{w}_1^{\left(n+1\right)}{\Vert }_{L^2{\left(\mathrm{\Omega }\times \left(0,T\right)\right)}^d}\right],\hfill \end{array}$

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\\ \hfill & \left.+\Vert {\mathrm{\Pi }}_{\mathcal{D}}{g}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}+\Vert {\mathrm{\Pi }}_{\mathcal{D}}{h}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}\right)+C_0\hfill \\ \hfill & +\left.\Vert {\delta }_{\mathcal{D}}^{\left(n+\frac{1}{2}\right)}{h}_{\mathcal{D}}^{\left(n+1\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\times \left(0,T\right)\right)}+{\mu }_1\Vert {\nabla }_{\mathcal{D}}{w}_2^{\left(n+1\right)}{\Vert }_{L^2{\left(\mathrm{\Omega }\times \left(0,T\right)\right)}^d}\right].\hfill \end{array}$

The desired estimates is then obtained by taking the supremum over $\phi ,\psi \in X_{\mathcal{D},0}$ with $\Vert {\nabla }_{\mathcal{D}}\phi {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^d}=\Vert {\nabla }_{\mathcal{D}}\psi {\Vert }_{L^2{\left(\mathrm{\Omega }\right)}^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 ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{g}_{{\mathcal{D}}_m}$, ${\mathrm{\Pi }}_{{\mathcal{D}}_m}{h}_{{\mathcal{D}}_m}$, ${\nabla }_{{\mathcal{D}}_m}{g}_{{\mathcal{D}}_m}$ and ${\nabla }_{{\mathcal{D}}_m}{h}_{{\mathcal{D}}_m}$, 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 $\mathcal{T}=\left(\mathcal{M},\mathcal{F},\mathcal{P},\mathcal{V}\right)$ 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

$X_{\mathcal{D},0}=\left\{v=\left({\left({\phi }_K\right)}_{K\in \mathcal{M}},{\left({\phi }_{\sigma }\right)}_{\sigma \in \mathcal{F}}\right):{\phi }_K,{\phi }_{\sigma }\in \mathbb{R},{\phi }_{\sigma }=0,\forall \sigma \in \mathcal{F}\cap \partial \mathrm{\Omega }\right\},$

$\begin{array}{cc}\hfill X_{\mathcal{D},\partial \mathrm{\Omega }}=& \left\{v=\left({\left(v_K\right)}_{K\in \mathcal{M}},{\left(v_{\sigma }\right)}_{\sigma \in \mathcal{F}}\right):v_K\in \mathbb{R},v_{\sigma }\in \mathbb{R},\right.\hfill \\ \hfill \text{mtd}\end{array}$

1. The non-conforming piecewise affine reconstruction ${\mathrm{\Pi }}_{\mathcal{D}}$ is defined by

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

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

${\nabla }_K\phi =\frac{1}{\left|K\right|}\sum _{\sigma \in {\mathcal{F}}_K}\left|\sigma \right|{\phi }_{\sigma }{\mathbf{n}}_{K,\sigma }\text{\hspace{0.17em}and\hspace{0.17em}}{R}_K\left(\phi \right)={\left({\phi }_{\sigma }-{\phi }_K-{\nabla }_K\phi \cdot \left(x_{\sigma }-x_K\right)\right)}_{\sigma \in {\mathcal{F}}_K},$

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 $\sigma \in {\mathcal{F}}_K$, n_K,σ is the unit vector normal to σ outward to K and D_K,σ is the convex hull of σ ∪ {x_K}.

1. The interpolant $J_{\mathcal{D}}:L^{\infty }\left(\mathrm{\Omega }\right)\to X_{\mathcal{D}}$ is defined by

$\begin{array}{c}\hfill \text{mtd}\\ \hfill \text{mtd}\end{array}$

1. the interpolant ${\mathcal{I}}_{\mathcal{D},\partial }:H^{\frac{1}{2}}\left(\partial \mathrm{\Omega }\right)\to X_{\mathcal{D},\partial \mathrm{\Omega }}$ is defined by

$\begin{array}{cc}\hfill \forall g\in H^{\frac{1}{2}}\left(\partial \mathrm{\Omega }\right):& {\left({\mathcal{I}}_{\mathcal{D},\partial }g\right)}_{\sigma }=\frac{1}{\left|\sigma \right|}{\int }_{\sigma }g\left(x\right)\mathrm{d}s\left(x\right),\hfill \\ \hfill \text{mtd}\end{array}$

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:

$\begin{array}{ll}\hfill F(\overline{u},\overline{v})& =a-\left(b+1\right)\overline{u}+{\overline{u}}^2\overline{v},\hfill \\ \hfill G(\overline{u},\overline{v})& =b\overline{u}-{\overline{u}}^2\overline{v},\hfill \end{array}$

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)$\overline{u}\left(\mathit{x},t\right)=\exp \left(-x-y-0.5t\right),$

(5.2b)$\overline{v}\left(\mathit{x},t\right)=\exp \left(x+y+0.5t\right),$

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 $\overline{u}$ and $\overline{v}$ and the corresponding rates of convergence with respect to the mesh size h. The resultant errors on the solutions $\overline{u}$ and $\overline{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 $\nabla \overline{u}$ and in Figure 1b for $\nabla \overline{v}$. A line of slope one is added in both figures as a reference. We observe that the relative errors on $\nabla \overline{u}$ and $\nabla \overline{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) $\nabla \overline{u}$ and (b) $\nabla \overline{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	$\frac{\Vert \overline{u}\left(\cdot ,t^{\left(n\right)}\right)-{\mathrm{\Pi }}_{\mathcal{D}}{u}^{\left(n\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\right)}}{\Vert \overline{u}\left(\cdot ,t^{\left(n\right)}\right){\Vert }_{L^2\left(\mathrm{\Omega }\right)}}$	Rate	$\frac{\Vert \overline{v}\left(\cdot ,t^{\left(n\right)}\right)-{\mathrm{\Pi }}_{\mathcal{D}}{v}^{\left(n\right)}{\Vert }_{L^2\left(\mathrm{\Omega }\right)}}{\Vert \overline{v}\left(\cdot ,t^{\left(n\right)}\right){\Vert }_{L^2\left(\mathrm{\Omega }\right)}}$	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