labelkeyrgb0.6,0,1 violetrgb0.580,0.,0.827

# A unified analysis for reaction diffusion models with application to the
spiral waves dynamics of the Barkley model

Yahya Alnashri  Department of Mathematics, Al-Qunfudah University College, Umm
Al-Qura University, Saudi Arabia<EMAIL_ADDRESS>and  Hasan Alzubaidi
Department of Mathematics, Al-Qunfudah University College, Umm Al-Qura
University, Saudi Arabia<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Applying the gradient discretisation method (GDM), the paper develops a
comprehensive numerical analysis for the reaction diffusion model. Using only
three properties, this analysis provides convergence results for several
conforming and non-conforming numerical schemes that align with the GDM. As an
application of this analysis, the hybrid mimetic mixed (HMM) method for the
reaction diffusion model is designed and its convergence established.
Numerical experiments using the HMM method are presented to facilitate study
of the creation of spiral waves in the Barkley model and the ways in which the
waves behave when interacting with the boundaries of their generating medium.

###### Key words and phrases:

Reaction diffusion models, Barkley model, spiral waves, homogeneous Neumann
boundary condition, gradient schemes, gradient discretisation, convergence
analysis, existence, hybrid mimetic mixed methods.

###### 2010 Mathematics Subject Classification:

35K57,65N12,65M08

##  1\. Introduction

Some arrhythmias can have fatal consequences if they are permitted to remain
in cardiac tissues without undergoing treatment [33, 13]. One method by which
it may be possible to develop a better understanding of how arrhythmias
develop and the most effective means by which they can be treated mathematical
modelling. A number of scholars have examined the mathematical modelling
approaches that are available and, following the application of these methods,
they have concluded that cardiac arrhythmias can be traced back to the free
rotation of spiral waves. These waves, reactivate an area of tissue at a
higher frequency than that associated with the normal sinoatrial node and,
thereby, generate a higher than average heartbeat [33, 43, 12, 29, 32]. In
some cases, the spiral waves break up into smaller spiral waves that stimulate
uncoordinated heart contractions that are referred to as fibrillation. When
fibrillation develops in the ventricles of the heart, it causes the heart to
tremor and its strength diminishes. Lacking the ability to pump blood at a
consistent rate, the heart ultimately suffers from a cardiac arrest [26].

It is possible to study the development and behaviour of spiral waves as a
form of spatio-temporal solution to reaction-diffusion equations of the
general form [34, 10],

$\displaystyle\partial_{t}\bar{u}(\boldsymbol{x},t)-\mu\mathop{\rm
div}(\nabla\bar{u}(\boldsymbol{x},t))$
$\displaystyle=f(\bar{u},\bar{v}),\quad(\boldsymbol{x},t)\in\Omega\times(0,T),$
(1.1) $\displaystyle\partial_{t}\bar{v}(\boldsymbol{x},t)$
$\displaystyle=g(\bar{u},\bar{v})\quad(\boldsymbol{x},t)\in\Omega\times(0,T),$
(1.2)

for diffusion coefficient $\mu$, with initial conditions and a pure Neumann
boundary condition are defined on the domain
$\Omega\subset{\mathbb{R}}^{d},d\geq 1$.

Although the model has been widely studied theoretically since it enters the
general theory of reaction diffusion equations [41], it is typically very rare
for researchers to be able to solve models that pertain to waves of this
nature. As such, there is an inherent need for numerical solutions that are
accurate and reliable. The existing literature describes the dynamics of
spiral waves using the reaction-diffusion models in depth. A variety of
mathematical approaches have been employed including finite difference methods
(FDM) [6, 5, 35, 16], pseudo-spectral methods [38, 31, 37], finite element
methods (FEM) [39, 9, 40], and finite volume methods (FVM) [30, 42, 14, 8, 7,
11, 15]. In this paper, we will focus on the literature related to FVM. In
this regard, the 1997 work of Harrild and Henriquez [30] is a particular note.
Harrild and Henriquez employed the elements of a FVM-based formulation to
examine the phenomenon of conduction in cardiac tissue. Trew et al. [42]
described the creation of a novel finite volume method that could effectively
be employed to model bidomain electrical activation in discontinuous cardiac
tissue. Coudire and Pierre[14] examined a 3D FVM that could be effectively
used to calculate the electrical activity that could be observed in the
myocardium on unstructured meshes and identified stability conditions for two
time-stepping methods in distinctive settings. They also generated error
estimates by which effective solutions of the monodomain reaction-diffusion
systems could be estimated. Bendahmane and Karlsen [8] successfully merged a
finite volume scheme with Dirichlet boundary conditions for the bidomain
model, representing a degenerate reaction-diffusion system that models the
electrophysiological waves that can be observed in cardiac tissue. Bendahmane
et al. [7] employed a comparable FVM with Neumann boundary conditions to
demonstrate the existence and uniqueness of the approximate solution for the
monodomain and bidomain models for the myocardial tissue electrical activity.
Burger et al. [11] presented some of fully space-time adaptive multiresolution
approaches that were based on a combination of FVM and Barkley’s approach for
modeling the complex dynamics of waves in excitable media. More recently in
2017, Coudire and Turpault [15], originated and evaluated a high-order FVM
approach in space combined with a high-order strong stability preserving (SSP)
Runge–Kutta technique to produce 2D simulations of spiral waves and simple
planar waves in cardiac tissue. Their research also yielded error estimates
for the given reaction-diffusion systems.

The chief elements of the spiral waves dynamics such as meandering and drift
under external perturbations are of great interest. External perturbations can
include light-induced drift, dual spiral interaction, and interaction of
spirals with a boundary [36]. When we scrutinize spiral-boundary interactions,
we can see two varieties of interaction, with the spiral waves being either
reflected or annihilated [36]. Examining spiral-boundary interactions and what
occurs in order to annihilate the spiral can assist us in creating techniques
to prevent arrhythmia generated by the spiral behaviour. For more details on
spiral drift due to the boundary effects, see [36, 28, 44] and the references
given there.

This paper’s main purpose is to use the gradient discretisation method (GDM)
to offer a complete and unified convergence analysis of the reaction-diffusion
model. To our knowledge, this analysis provides the first results that are
applicable to several conforming and non conforming methods. The GDM is an
abstract setting to study the numerical analysis for linear or non linear,
steady or time–dependent diffusion partial differential equations (PDEs).
Based on a limit number of properties, the GDM can establish convergence of
numerical schemes for various models with different boundary conditions.
Various studies have established that the GDM covers several families of
numerical schemes: conforming, non conforming and mixed finite elements
methods (including the non-conforming, Crouzeix Raviart method and the Raviart
Thomas method), the discontinuous Galerkin scheme, the vertex approximate
gradient (VAG) scheme, hybrid mimetic mixed methods (which contain hybrid
mimetic finite differences, hybrid finite volumes/SUSHI scheme and mixed
finite volumes), nodal mimetic finite differences, and finite volumes methods
(such as some multi-points flux approximation and discrete duality finite
volume methods). For more details, see [21, 24, 22, 25, 23, 18, 3, 2] and the
monograph [19] for a complete presentation.

This paper is organised as follows. Section 2 is devoted to the continuous
model. Section 3 provides the discrete elements of the GDM and two properties
required to analyse the studied model. It also states the discrete problem
(the gradient scheme) followed by the main results: existence and uniqueness
of the solution to the scheme, and its convergence to the weak solution of the
studied model. Section 4 explains the role that the GDM plays in applying the
the hybrid mixed mimetic (HMM) methods, and establishing its convergence for
the considered model. Section 5 is devoted to proving the main novelty of this
paper, Theorem 3.6, which is obtained via the compactness technique. In
Section 6 we study numerically using the HMM, the propagation of the spiral
waves for the Barkley model and their behavior when they interact with the
boundary of the medium where they spread.

##  2\. Continuous model

We consider the following system of partial differential equations,
representing varieties models such as the Barkley model [6], the
FitzHugh–Nagumo model [27], the Aliev and Panfilov model [1] and the
Belousov–Zhabotinsky reaction [45].

$\displaystyle\partial_{t}\bar{u}(\boldsymbol{x},t)-\mu\mathop{\rm
div}(\nabla\bar{u}(\boldsymbol{x},t))$
$\displaystyle=f(\bar{u},\bar{v}),\quad(\boldsymbol{x},t)\in\Omega\times(0,T),$
(2.1) $\displaystyle\partial_{t}\bar{v}(\boldsymbol{x},t)$
$\displaystyle=g(\bar{u},\bar{v}),\quad(\boldsymbol{x},t)\in\Omega\times(0,T),$
(2.2) $\displaystyle\nabla\bar{u}(\boldsymbol{x},t)\cdot\mathbf{n}$
$\displaystyle=0,\quad(\boldsymbol{x},t)\in\partial\Omega\times(0,T),$ (2.3)
$\displaystyle\bar{u}(\boldsymbol{x},0)$ $\displaystyle=u_{\rm
ini}(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega,$ (2.4)
$\displaystyle\bar{v}(\boldsymbol{x},0)$ $\displaystyle=v_{\rm
ini}(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega,$ (2.5)

where $\mathbf{n}$ is the outer normal to $\partial\Omega$. In this excitable
media system, the unknowns $\bar{u}(\boldsymbol{x},t)$ and
$\bar{v}(\boldsymbol{x},t)$ denote the excitation and recovery terms,
respectively. A dynamic of particular case of this model is explained in
Section (6).

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.5) are the following:

  1. (1)

$\Omega$ is an open bounded connected subset of ${\mathbb{R}}^{d}\;(d\geq 1)$,
with a Lipschitz boundary, $T>0$, $\mu\in{\mathbb{R}}$,

  2. (2)

$(u_{\rm ini},v_{\rm ini})$ are in $H^{1}(\Omega)\times L^{\infty}(\Omega)$,

  3. (3)

the functions $f,\;g:{\mathbb{R}}^{2}\to{\mathbb{R}}$ are in
$L^{2}(\Omega\times(0,T))$ and both polynomial functions satisfying the
following conditions:

     * •

they are affine with respect to $\xi$;

$f(s,\xi)=f_{1}(s)+f_{2}(s)\xi,\quad g(s,\xi)=g_{1}(s)+\alpha\xi,$ (2.6)

where $f_{1},f_{2},g_{1}$ are continuous functions defined on ${\mathbb{R}}$
and $\alpha$ is a constant.

     * •

there exists constants $c_{i}\geq 0(i=1...5)$ such that for any
$s\in{\mathbb{R}}$,

$|f_{1}(s)|\leq c_{1}+c_{2}|s|,$ (2.7a) $|f_{2}(s)|\leq c_{3},$ (2.7b)
$|g_{1}(s)|\leq c_{4}+c_{5}|s|.$ (2.7c)

######  Definition 2.2 (Weak solution).

Under Assumptions 2.1, a pair $(\bar{u},\bar{v})$ is said to be a weak
solution of Problem (2.1)–(2.5) if the following properties and equalities
hold:

$\displaystyle\bar{u}\in L^{2}(0,T;H^{1}(\Omega))\cap
C([0,T];L^{2}(\Omega)),\;\bar{v}\in L^{2}(0,T;L^{2}(\Omega)),$
$\displaystyle\partial_{t}\bar{u}\in
L^{2}(\Omega\times(0,T)),\;\bar{u}(\cdot,0)=u_{\rm ini},$
$\displaystyle\mbox{and for all $\bar{\varphi}\in L^{2}(0,T;H^{1}(\Omega))$
and for all $\bar{\psi}\in L^{2}(0,T;L^{2}(\Omega))$},$
$\displaystyle\mbox{such that $\partial_{t}\bar{\psi}\in
L^{2}(0,T;L^{2}(\Omega))$ and $\bar{\psi}(\cdot,T)=0$},$
$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}\bar{u}(\boldsymbol{x},t)\bar{\varphi}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t+\mu\int_{0}^{T}\int_{\Omega}\nabla\bar{u}(\boldsymbol{x},t)\cdot\nabla\bar{\varphi}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$
(2.8a)
$\displaystyle\quad\quad\quad\quad\quad=\int_{0}^{T}\int_{\Omega}f(\bar{u},\bar{v})\bar{\varphi}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$
$\displaystyle-\int_{0}^{T}\int_{\Omega}\bar{v}(\boldsymbol{x},t)\partial_{t}\bar{\psi}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t-\int_{\Omega}v_{\rm
ini}(\boldsymbol{x})\bar{\psi}(\boldsymbol{x},0)\,\mathrm{d}x$ (2.8b)
$\displaystyle\qquad\quad\quad\quad\quad=\int_{0}^{T}\int_{\Omega}g(\bar{u},\bar{v})\bar{\psi}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

######  Remark 2.3.

The existence of at least one solution $(\bar{u},\bar{v})$ to (2.8) will be a
consequence of the convergence analysis of the gradient discretisation method.

##  3\. Discrete Problem

As stated in the introduction section, the analysis of numerical schemes for
the approximation of solutions to the reaction diffusion model is performed
using the gradient discretisation method. This method first requires the
reconstruction of a set of discrete spaces and operators, which is called
gradient discretisation (GD).

######  Definition 3.1 (gradient discretisation for reaction diffusion model).

Let $\Omega$ be an open subset of ${\mathbb{R}}^{d}$ (with $d\geq 1$) and
$T>0$. A gradient discretisation for the reaction–diffusion model is
${\mathcal{D}}=(X_{{\mathcal{D}}},Y_{\mathcal{D}},\Pi_{\mathcal{D}},\Pi_{{\mathcal{D}}^{{}^{\prime}}},\nabla_{\mathcal{D}},J_{\mathcal{D}},J_{{\mathcal{D}}^{{}^{\prime}}},(t^{(n)})_{n=0,...,N}))$,
where

  * •

the two set of discrete unknowns $X_{{\mathcal{D}}}$ and $Y_{\mathcal{D}}$ are
finite dimensional vector spaces on ${\mathbb{R}}$,

  * •

the function reconstruction $\Pi_{\mathcal{D}}:X_{{\mathcal{D}}}\to
L^{2}(\Omega)$ is a linear,

  * •

the function reconstruction
$\Pi_{{\mathcal{D}}^{{}^{\prime}}}:Y_{{\mathcal{D}}}\to L^{2}(\Omega)$ is a
linear, and must be defined so that
$||\Pi_{{\mathcal{D}}^{{}^{\prime}}}\cdot||_{L^{2}(\Omega)}$ is a norm on
$Y_{\mathcal{D}}$,

  * •

the gradient reconstruction $\nabla_{\mathcal{D}}:X_{{\mathcal{D}}}\to
L^{2}(\Omega)^{d}$ is a linear and must be defined so that

$\|\varphi||_{\mathcal{D}}=||\Pi_{\mathcal{D}}\varphi||_{L^{2}(\Omega)}+||\nabla_{\mathcal{D}}\varphi||_{L^{2}(\Omega)^{d}}$
(3.1)

is a norm on $X_{{\mathcal{D}}}$.

  * •

$J_{\mathcal{D}}:H^{1}(\Omega)\to X_{\mathcal{D}}$ and
$J_{{\mathcal{D}}^{{}^{\prime}}}:L^{\infty}(\Omega)\to Y_{\mathcal{D}}$ are a
linear and continuous interpolation operator for the initial conditions,

  * •

$t^{(0)}=0<t^{(1)}<....<t^{(N)}=T$ are discrete times.

Let us introduce some notations to define the space–time reconstructions
$\Pi_{\mathcal{D}}\varphi:\Omega\times[0,T]\to{\mathbb{R}}$, and
$\nabla_{\mathcal{D}}\varphi:\Omega\times[0,T]\to{\mathbb{R}}^{d}$, and the
discrete time derivative $\delta_{\mathcal{D}}\varphi:(0,T)\to L^{2}(\Omega)$,
for $\varphi=(\varphi^{(n)})_{n=0,...,N}\in X_{\mathcal{D}}^{N}$ and
$\psi=(\psi^{(n)})_{n=0,...,N}\in Y_{\mathcal{D}}^{N}$

For a.e $\boldsymbol{x}\in\Omega$, for all $n\in\\{0,...,N-1\\}$ and for all
$t\in(t^{(n)},t^{(n+1)}]$, let

$\begin{split}&\Pi_{\mathcal{D}}\varphi(\boldsymbol{x},0)=\Pi_{\mathcal{D}}\varphi^{(0)}(\boldsymbol{x}),\quad\Pi_{\mathcal{D}}\varphi(\boldsymbol{x},t)=\Pi_{\mathcal{D}}\varphi^{(n+1)}(\boldsymbol{x}),\\\
&\nabla_{\mathcal{D}}\varphi(\boldsymbol{x},t)=\nabla_{\mathcal{D}}\varphi^{(n+1)}(\boldsymbol{x}),\\\
&\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi(\boldsymbol{x},0)=\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi^{(0)}(\boldsymbol{x}),\mbox{
and
}\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi(\boldsymbol{x},t)=\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi^{(n+1)}(\boldsymbol{x}).\end{split}$

Set $\delta t^{(n+\frac{1}{2})}=t^{(n+1)}-t^{(n)}$ and $\delta
t_{\mathcal{D}}=\max_{n=0,...,N-1}\delta t^{(n+\frac{1}{2})}$, to define

$\delta_{\mathcal{D}}\varphi(t)=\delta_{\mathcal{D}}^{(n+\frac{1}{2})}\varphi:=\frac{\Pi_{\mathcal{D}}(\varphi^{(n+1)}-\varphi^{(n)})}{\delta
t^{(n+\frac{1}{2})}}\mbox{ and
}\delta_{{\mathcal{D}}^{{}^{\prime}}}\psi(t)=\delta_{{\mathcal{D}}^{{}^{\prime}}}^{(n+\frac{1}{2})}\psi:=\frac{\Pi_{{\mathcal{D}}^{{}^{\prime}}}(\psi^{(n+1)}-\psi^{(n)})}{\delta
t^{(n+\frac{1}{2})}}.$

In order to build converging schemes for models including homogenous Neumann
boundary conditions, the elements of a GD must enjoy as much as possible the
properties of the continuous space and operators, coercivity, consistency and
limit-conformity. In our model, there is no need to the coercivity property,
which gives the discrete Poincare inequality, since the time derivative in our
model plays the same role in establishing $L^{2}$ estimate on the discrete
solution.

A consistency of the method (refers to an interplant error), which ensures the
accuracy of approximating smooth functions and their gradients by elements
defined from the discrete space, and the convergence of the time steps to
zero.

######  Definition 3.2 (Consistency).

If ${\mathcal{D}}$ is a gradient discretisation in the sense of Definition
3.1, define $S_{\mathcal{D}}:H^{1}(\Omega)\to[0,+\infty)$ and
$S_{{\mathcal{D}}^{{}^{\prime}}}:L^{2}(\Omega)\to[0,+\infty)$ by

$\forall\varphi\in H^{1}(\Omega),\;S_{\mathcal{D}}(\varphi)=\min_{w\in
X_{\mathcal{D}}}\left(\|\Pi_{\mathcal{D}}w-\varphi\|_{L^{2}(\Omega)}+\|\nabla_{\mathcal{D}}w-\nabla\varphi\|_{L^{2}(\Omega)^{d}}\right),$
(3.2)

and

$\forall\psi\in
L^{2}(\Omega),\;S_{{\mathcal{D}}^{{}^{\prime}}}(\psi)=\min_{w\in
Y_{\mathcal{D}}}\|\Pi_{{\mathcal{D}}^{{}^{\prime}}}w-\psi\|_{L^{2}(\Omega)}.$
(3.3)

A sequence $(\mathcal{D}_{m})_{m\in\mathbb{N}}$ of gradient discretisations is
_consistent_ if, as $m\to\infty$

  * •

for all $\varphi\in H^{1}(\Omega)$, $S_{{\mathcal{D}}_{m}}(\varphi)\to 0$, and
for all $\psi\in L^{2}(\Omega)$, $S_{{\mathcal{D}}_{m}^{{}^{\prime}}}(\psi)\to
0$,

  * •

for all $w\in L^{2}(\Omega)$,
$\Pi_{{\mathcal{D}}_{m}}J_{{\mathcal{D}}_{m}}w\to w$ in $L^{2}(\Omega)$,

  * •

for all $w\in L^{\infty}(\Omega)$,
$(\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}J_{{\mathcal{D}}_{m}^{{}^{\prime}}}w)_{m\in{\mathbb{N}}}$
is bounded in $L^{\infty}(\Omega)$ and converges to $w$ in $L^{2}(\Omega)$,

  * •

$\delta t_{{\mathcal{D}}_{m}}\to 0$.

The quantity $W_{\mathcal{D}}$ defined below measures how well the discrete
Stokes formula is satisfied, it is only exactly ensured in conforming methods.

######  Definition 3.3 (Limit-conformity).

If ${\mathcal{D}}$ is a gradient discretisation in the sense of Definition
3.1, and $H_{\rm div}=\\{{\boldsymbol{\psi}}\in L^{2}(\Omega)^{d}\;:\;{\rm
div}{\boldsymbol{\psi}}\in
L^{2}(\Omega),\;{\boldsymbol{\psi}}\cdot\mathbf{n}=0\mbox{ on
}\partial\Omega\\}$, define $W_{\mathcal{D}}:H_{\rm div}\to[0,+\infty)$ by

$\forall{\boldsymbol{\psi}}\in H_{\rm div},\quad
W_{\mathcal{D}}({\boldsymbol{\psi}})=\sup_{w\in
X_{\mathcal{D}}\setminus\\{0\\}}\frac{\Big{|}\displaystyle\int_{\Omega}(\nabla_{\mathcal{D}}w\cdot{\boldsymbol{\psi}}+\Pi_{\mathcal{D}}w\mathop{\rm
div}({\boldsymbol{\psi}}))\,\mathrm{d}x\Big{|}}{||w||_{\mathcal{D}}}.$ (3.4)

A sequence $({\mathcal{D}}_{m})_{m\in{\mathbb{N}}}$ of gradient
discretisations is _limit-conforming_ if for all ${\boldsymbol{\psi}}\in
H_{\rm div}$, $W_{{\mathcal{D}}_{m}}({\boldsymbol{\psi}})\to 0$, as
$m\to\infty$.

Finally, dealing with non linearity requires the operators $\Pi_{\mathcal{D}}$
and $\nabla_{\mathcal{D}}$ to provide the compactness properties, defined
below. Since $f$ and $g$ are assumed to be affine functions, there is no need
for the operator $\Pi_{{\mathcal{D}}^{{}^{\prime}}}$ to afford this property.

######  Definition 3.4 (Compactness).

A sequence of gradient discretisation ${\mathcal{D}}_{m}$ in the sense of
Definition 3.1 is _compact_ if for any sequence
$(\varphi_{m})_{m\in{\mathbb{N}}}\in X_{{\mathcal{D}}_{m}}$, such that
$(||\varphi_{m}||_{{\mathcal{D}}_{m}})_{m\in{\mathbb{N}}}$ is bounded, the
sequence $(\Pi_{{\mathcal{D}}_{m}}\varphi_{m})_{m\in{\mathbb{N}}}$ is
relatively compact in $L^{2}(\Omega)$.

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 (GS).

######  Definition 3.5 (Gradient scheme for (2.8)).

Find families $(u^{(n)})_{n=0,...,N}\in X_{\mathcal{D}}^{N+1}$ and
$(v^{(n)})_{n=0,...,N}\in Y_{\mathcal{D}}^{N+1}$, such that
$u^{(0)}=J_{\mathcal{D}}u_{\rm ini}$ and
$v^{(0)}=J_{{\mathcal{D}}^{{}^{\prime}}}v_{\rm ini}$, and for all
$n=0,...,N-1$, $u^{(n+1)}$ and $v^{(n+1)}$ satisfy

$\displaystyle\int_{\Omega}\delta_{\mathcal{D}}^{(n+\frac{1}{2})}u(\boldsymbol{x})\Pi_{\mathcal{D}}\varphi(\boldsymbol{x})+\mu\int_{\Omega}\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\cdot\nabla_{\mathcal{D}}\varphi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
(3.5a)
$\displaystyle\qquad=\int_{\Omega}f_{\mathcal{D}}^{(n+1)}\Pi_{\mathcal{D}}\varphi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x},\quad\forall\varphi\in
X_{\mathcal{D}},$
$\displaystyle\int_{\Omega}\delta_{{\mathcal{D}}^{{}^{\prime}}}^{(n+\frac{1}{2})}v(\boldsymbol{x})\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi(\boldsymbol{x})=\displaystyle\int_{\Omega}g_{\mathcal{D}}^{(n+1)}\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x},\quad\forall\psi\in
Y_{\mathcal{D}},$ (3.5b)

where the discrete reaction terms are defined by

$f_{\mathcal{D}}^{(n+1)}=f(\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x}),\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x}))\mbox{
and
}g_{\mathcal{D}}^{(n+1)}=g(\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x}),\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x})).$

The convergence results of this gradient scheme is stated in the following
theorem, whose proof is detailed in Section 5.

######  Theorem 3.6 (Convergence of the GS for the reaction diffusion model).

Assume (2.1) and Let $({\mathcal{D}}_{m})_{m\in{\mathbb{N}}}$ be a sequence of
gradient discretisations in the sense of Definition 3.1, that is consistent,
limit-conforming and compact in the sense of Definitions 3.2, 3.3 and 3.4. For
$m\in{\mathbb{N}}$, let $(u_{m},v_{m})$ be a solution to the gradient scheme
(3.5) with ${\mathcal{D}}={\mathcal{D}}_{m}$ and
$||\nabla_{{\mathcal{D}}_{m}}J_{{\mathcal{D}}_{m}}u_{\rm
ini}||_{L^{2}(\Omega)^{d}}$ is bounded. Then there exists a weak solution
$(\bar{u},\bar{v})$ of (2.8) and a subsequence of gradient discretisations,
still denoted by $({\mathcal{D}}_{m})_{m\in{\mathbb{N}}}$, such that, as
$m\to\infty$,

  1. (1)

$\Pi_{{\mathcal{D}}_{m}}u_{m}$ converges strongly to $\bar{u}$ in
$L^{\infty}(0,T;L^{2}(\Omega))$,

  2. (2)

$\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}$ converges weakly to $\bar{v}$ in
$L^{\infty}(0,T;L^{2}(\Omega))$,

  3. (3)

$\nabla_{{\mathcal{D}}_{m}}u_{m}$ converges strongly to $\nabla\bar{u}$ in
$L^{2}(\Omega\times(0,T))^{d}$.

##  4\. The hybrid mixed mimetic (HMM) methods

We show here that hybrid mixed mimetic (HMM) methods can be expressed as
gradient schemes formulation when applied to the reaction diffusion model. It
is shown in [20] that the HMM methods is a generic setting gathering three
different families of methods, namely the hybrid finite volume method, the
(mixed-hybrid) mimetic finite differences methods, and the mixed finite volume
methods. In order to construct this mesh–based method, we recall here the
notion of polytopal mesh [20].

######  Definition 4.1 (Polytopal mesh).

Let $\Omega$ be a bounded polytopal open subset of ${\mathbb{R}}^{d}$ ($d\geq
1$). A polytopal mesh of $\Omega$ is given by
${\mathcal{T}}=({\mathcal{M}},{\mathcal{E}},\mathcal{P})$, where:

  1. (1)

${\mathcal{M}}$ is a finite family of non empty connected polytopal open
disjoint subsets of $\Omega$ (the cells) such that
$\overline{\Omega}=\displaystyle{\cup_{K\in{\mathcal{M}}}\overline{K}}$. For
any $K\in{\mathcal{M}}$, $|K|>0$ is the measure of $K$ and $h_{K}$ denotes the
diameter of $K$.

  2. (2)

${\mathcal{E}}$ is a finite family of disjoint subsets of $\overline{\Omega}$
(the edges of the mesh in 2D, the faces in 3D), such that any
$\sigma\in{\mathcal{E}}$ is a non empty open subset of a hyperplane of
${\mathbb{R}}^{d}$ and $\sigma\subset\overline{\Omega}$. We assume that for
all $K\in{\mathcal{M}}$ there exists a subset ${{{\mathcal{E}}}_{K}}$ of
${\mathcal{E}}$ such that $\partial
K=\displaystyle{\cup_{\sigma\in{{{\mathcal{E}}}_{K}}}}\overline{\sigma}$. We
then set
${\mathcal{M}}_{\sigma}=\\{K\in{\mathcal{M}}\,:\,\sigma\in{{{\mathcal{E}}}_{K}}\\}$
and assume that, for all $\sigma\in{\mathcal{E}}$, ${\mathcal{M}}_{\sigma}$
has exactly one element and $\sigma\subset\partial\Omega$, or
${\mathcal{M}}_{\sigma}$ has two elements and $\sigma\subset\Omega$.
${{{\mathcal{E}}}_{\rm int}}$ is the set of all interior faces, i.e.
$\sigma\in{\mathcal{E}}$ such that $\sigma\subset\Omega$, and
${{{\mathcal{E}}}_{\rm ext}}$ the set of boundary faces, i.e.
$\sigma\in{\mathcal{E}}$ such that $\sigma\subset\partial\Omega$. For
$\sigma\in{\mathcal{E}}$, the $(d-1)$-dimensional measure of $\sigma$ is
$|\sigma|$, the centre of mass of $\sigma$ is $\overline{x}_{\sigma}$, and the
diameter of $\sigma$ is $h_{\sigma}$.

  3. (3)

$\mathcal{P}=(x_{K})_{K\in{\mathcal{M}}}$ is a family of points of $\Omega$
indexed by ${\mathcal{M}}$ and such that, for all $K\in{\mathcal{M}}$,
$x_{K}\in K$ ($x_{K}$ is sometimes called the “centre” of $K$). We then assume
that all cells $K\in{\mathcal{M}}$ are strictly $x_{K}$-star-shaped, meaning
that if $x\in\overline{K}$ then the line segment $[x_{K},x)$ is included in
$K$.

For a given $K\in{\mathcal{M}}$, let $\mathbf{n}_{K,\sigma}$ be the unit
vector normal to $\sigma$ outward to $K$ and denote by $d_{K,\sigma}$ the
orthogonal distance between $x_{K}$ and $\sigma\in\mathcal{E}_{K}$. The size
of the discretisation is
$h_{\mathcal{M}}=\sup\\{h_{K}\,:\;K\in{\mathcal{M}}\\}$.

Let ${\mathcal{T}}$ be a polytopal mesh. A gradient discretisation
$(X_{\mathcal{D}},Y_{\mathcal{D}},\Pi_{\mathcal{D}},\Pi_{{\mathcal{D}}^{{}^{\prime}}},\nabla_{\mathcal{D}},J_{\mathcal{D}},J_{{\mathcal{D}}^{{}^{\prime}}})$
in the setting of HMM method formats are then constructed by setting

$\displaystyle
X_{{\mathcal{D}}}=\\{\varphi=((\varphi_{K})_{K\in\mathcal{M}},(\varphi_{\sigma})_{\sigma\in\mathcal{E}})\;:\;\varphi_{K}\in{\mathbb{R}},\,\varphi_{\sigma}\in{\mathbb{R}}\\},$
$\displaystyle
Y_{{\mathcal{D}}}=\\{\varphi=(\varphi_{K})_{K\in\mathcal{M}}\;:\;\varphi_{K}\in{\mathbb{R}}\\},$
$\displaystyle\forall
K\in\mathcal{M}\,:\,\Pi_{\mathcal{D}}\varphi=\Pi_{{\mathcal{D}}^{{}^{\prime}}}\varphi=\varphi_{K}\mbox{
on $K$},$ $\displaystyle\forall w\in
L^{2}(\Omega)\;:\;J_{\mathcal{D}}w=((w_{K})_{K\in{\mathcal{M}}},(w_{\sigma})_{\sigma\in{\mathcal{E}}_{K}})\in
X_{\mathcal{D}},$ $\displaystyle\forall w\in
L^{\infty}(\Omega),\;J_{{\mathcal{D}}^{{}^{\prime}}}w=(w_{K})_{K\in{\mathcal{M}}}\in
Y_{\mathcal{D}},$ $\displaystyle\mbox{where
}w_{k}=\frac{1}{|K|}\int_{K}w(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\mbox{
and }w_{\sigma}=0,$ $\displaystyle\forall\varphi\in X_{\mathcal{D}},\;\forall
K\in\mathcal{M},\,\forall\sigma\in\mathcal{E}_{K},$
$\displaystyle\nabla_{\mathcal{D}}\varphi=\nabla_{K}\varphi+\frac{\sqrt{d}}{d_{K,\sigma}}R_{K}(\varphi)\mathbf{n}_{K,\sigma}\mbox{
on }D_{K,\sigma},$

where where a cell–wise constant gradient $\nabla_{K}(\varphi)$ and a
stabilisation term $R_{K}(\varphi)$ are respectively defined by:

$\nabla_{K}\varphi=\frac{1}{|K|}\sum_{\sigma\in{{{\mathcal{E}}}_{K}}}|\sigma|\varphi_{\sigma}\mathbf{n}_{K,\sigma}\mbox{
and
}R_{K}(\varphi)=(\varphi_{\sigma}-\varphi_{K}-\nabla_{K}\varphi\cdot(\overline{x}_{\sigma}-x_{K}))_{\sigma\in{\mathcal{E}}_{K}}.\\\
$

The gradient scheme (3.5) coming from such a GD can be given by: find
$(u^{(n)})_{n=0,...,N}\in X_{\mathcal{D}}^{N+1}$ and $(v^{(n)})_{n=0,...,N}\in
Y_{\mathcal{D}}^{N+1}$, such that $u^{(0)}=J_{\mathcal{D}}u_{\rm ini}$ and
$v^{(0)}=J_{{\mathcal{D}}^{{}^{\prime}}}v_{\rm ini}$, and for all
$n=0,...,N-1$, $u^{(n+1)}$ and $v^{(n+1)}$ satisfy, for all $\varphi\in
X_{\mathcal{D}}$ and for all $\psi\in Y_{\mathcal{D}}$,

$\left.\begin{aligned} &\displaystyle\sum_{K\in{\mathcal{M}}}\frac{|K|}{\delta
t^{(n+\frac{1}{2})}}\Big{(}u_{K}^{(n+1)}-u_{K}^{(n)}\Big{)}\varphi_{K}+\displaystyle\sum_{K\in{\mathcal{M}}}|K|\Lambda_{K}\nabla_{K}u^{(n+1)}\cdot\nabla_{K}\varphi\\\
&\quad+\displaystyle\sum_{K\in{\mathcal{M}}}(R_{K}\varphi)^{T}\mathbb{B}_{K}R_{K}(u^{(n+1)})=\displaystyle\sum_{K\in{\mathcal{M}}}\varphi_{K}\displaystyle\int_{K}f(u_{K}^{(n+1)},v_{K}^{(n+1)})\,\mathrm{d}\boldsymbol{x},\\\
&\displaystyle\sum_{K\in{\mathcal{M}}}\frac{|K|}{\delta
t^{(n+\frac{1}{2})}}\Big{(}v_{K}^{(n+1)}-v_{K}^{(n)}\Big{)}\psi_{K}=\displaystyle\sum_{K\in{\mathcal{M}}}\psi_{K}\displaystyle\int_{K}g(u_{K}^{(n+1)},v_{K}^{(n+1)})\,\mathrm{d}\boldsymbol{x},\end{aligned}\right.$
(4.1)

where $\mathbb{B}_{K}$ is a symmetric positive definite matrix of size
$\mbox{Card}({\mathcal{E}}_{K})$.

The above HMM scheme can be formulated as a finite volume scheme (a format
used in the computational processes). To do so, introduce the linear fluxes
$u\mapsto F_{K,\sigma}(u)$ (for $K\in{\mathcal{M}}$ and
$\sigma\in{\mathcal{E}}_{K}$) defined in [2]: for all $K\in{\mathcal{M}}$ and
all $u,w\in X_{\mathcal{D}}$,

$\displaystyle\sum_{\sigma\in\mathcal{E}_{K}}|\sigma|F_{K,\sigma}(u)(w_{K}-w_{\sigma})={}$
$\displaystyle\mu\int_{K}\nabla_{\mathcal{D}}u\cdot\nabla_{\mathcal{D}}w\,\mathrm{d}\boldsymbol{x}.$

Then Problem (4.1) can be recast as, for all $n=0,...,N-1$,

$\displaystyle\frac{|K|}{\delta
t^{(n+\frac{1}{2})}}\Big{(}u^{(n+1)}-u^{(n)}\Big{)}$
$\displaystyle+\sum_{\sigma\in\mathcal{E}_{K}}|\sigma|F_{K,\sigma}(u^{(n+1)})=|K|f_{K}^{(n+1)},\quad\forall
K\in{\mathcal{M}}$ $\displaystyle\frac{|K|}{\delta
t^{(n+\frac{1}{2})}}\Big{(}v^{(n+1)}-v^{(n)}\Big{)}$
$\displaystyle=|K|g_{K}^{(n+1)},\quad\forall K\in{\mathcal{M}}$ $\displaystyle
F_{K,\sigma}(u^{(n+1)})+F_{L,\sigma}(u^{(n+1)})$
$\displaystyle=0,\quad\forall\sigma\in\mathcal{E}_{\rm int}\mbox{ with
}{\mathcal{M}}_{\sigma}=\\{K,L\\},$ $\displaystyle F_{K,\sigma}(u^{(n+1)})$
$\displaystyle=0,\quad\forall
K\in{\mathcal{M}}\,,\forall\sigma\in\mathcal{E}_{K}\mbox{ such that
}\sigma\subset\partial\Omega,$

where
$f_{K}^{(n+1)}=\frac{1}{|K|}\int_{K}f(u_{K}^{(n+1)},v_{K}^{(n+1)})\,\mathrm{d}\boldsymbol{x}$
and
$g_{K}^{(n+1)}=\frac{1}{|K|}\int_{K}g(u_{K}^{(n+1)},v_{K}^{(n+1)})\,\mathrm{d}\boldsymbol{x},\;\forall
K\in{\mathcal{M}}$.

Now, let us discuss the convergence of the scheme (4.1). Assume the existence
of $\theta>0$ such that, for any $m\in{\mathbb{N}}$,

$\displaystyle{\theta_{\mathcal{M}}:=}\max_{K\in{\mathcal{M}}_{m}}$
$\displaystyle\left(\max_{\sigma\in{\mathcal{E}}_{K}}\frac{h_{K}}{d_{K,\sigma}}+{\rm
Card}({\mathcal{E}}_{K})\right)+\max_{\stackrel{{\scriptstyle\mbox{\scriptsize$\sigma\in{\mathcal{E}}_{m,\rm
int}$}}}{{{\mathcal{M}}_{\sigma}=\\{K,L\\}}}}\left(\frac{d_{K,\sigma}}{d_{L,\sigma}}+\frac{d_{L,\sigma}}{d_{K,\sigma}}\right)\leq\theta$
(4.2)

and, for all $K\in{\mathcal{M}}_{m}$ and
$\gamma\in{\mathbb{R}}^{{\mathcal{E}}_{K}}$,

$\displaystyle\frac{1}{\theta}\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|\left|\frac{R_{K,\sigma}(\gamma)}{d_{K,\sigma}}\right|^{2}\leq{}$
$\displaystyle\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|\left|\frac{(A_{K}R_{K}(\gamma))_{\sigma}}{d_{K,\sigma}}\right|^{2}$
(4.3) $\displaystyle\leq{}$
$\displaystyle\theta\sum_{\sigma\in{\mathcal{E}}_{K}}|D_{K,\sigma}|\left|\frac{R_{K,\sigma}(\gamma)}{d_{K,\sigma}}\right|^{2}.$

Under the above boundedness assumptions on the mesh regularity parameter
$\theta_{\mathcal{M}}$, the limit-conformity and the compactness properties
defined in Section 3 directly follow from the results in [19, Theorem 13.18].
It remains to prove the consistency. Proving that, for all $\bar{\varphi}\in
H^{1}(\Omega)$, $\lim_{m\to\infty}S_{{\mathcal{D}}_{m}}(\bar{\varphi})=0$ is
as in the case of the HMM method for PDEs, see [19, Theorem 4.14]. Also, with
the same manner, we can show that, for all $\bar{\psi}\in L^{2}(\Omega)$,
$\widehat{S}_{{\mathcal{D}}_{m}}(\bar{\psi})\to 0$, as $m\to\infty$. Let
$\varphi_{m}=((\varphi_{K})_{K\in{\mathcal{M}}},(\varphi_{\sigma})_{\sigma\in{\mathcal{E}}})\in
X_{{\mathcal{D}}_{m}}$ and $\psi_{m}=(\psi_{K})_{K\in{\mathcal{M}}}\in
Y_{{\mathcal{D}}_{m}}$ be the two interplants such that
$\varphi_{m}=J_{{\mathcal{D}}_{m}}u_{\rm{ini}}$ and
$\psi_{m}=J_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{\rm{ini}}$. Using [19,
Estimate (B.11), in Lemma B.6] with $p=2$, we can shows that $||u_{\rm
ini}-\Pi_{{\mathcal{D}}_{m}}J_{{\mathcal{D}}_{m}}u_{\rm
ini}||_{L^{2}(\Omega)}\to 0$ and $||v_{\rm
ini}-\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}J_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{\rm
ini}||_{L^{2}(\Omega)}\to 0$, as $m\to\infty$, which completes the consistency
property.

It is proved in [19, Theorem 13.14] that, for $w\in W^{1,p}(\Omega)$, there
exists an interpolant
$\varphi_{m}=((\varphi_{K})_{K\in{\mathcal{M}}},(\varphi_{\sigma})_{\sigma\in{\mathcal{E}}})\in
X_{{\mathcal{D}}_{m}}$, such that there is a constant $C_{1}>0$ not depending
on $m$, and

$||\nabla_{{\mathcal{D}}_{m}}\varphi_{m}||_{L^{p}(\Omega)^{d}}\leq
C_{1}||\nabla w||_{L^{p}(\Omega)^{d}}.$

The boundedness of $||\nabla_{{\mathcal{D}}_{m}}J_{{\mathcal{D}}_{m}}u_{\rm
ini}||_{L^{p}(\Omega)^{d}}$ is a consequence of applying the above estimate to
$w=u_{\rm ini}$ (with $p=2$). The convergence of the HMM scheme for the
reaction diffusion model is a consequence of Theorem 3.6 if the discrete time
steps $\delta t_{{\mathcal{D}}_{m}}^{(n+\frac{1}{2})}$ tends to 0, as
$m\to\infty$.

##  5\. Proof of the convergence results

In order to prove the convergence results, we first establish some
preliminaries estimates on the solution and its gradient of the scheme.

######  Lemma 5.1 (Estimates).

Under Assumptions (2.1), let ${\mathcal{D}}$ be a gradient discretisation in
the sense of Definition 3.1 and let $(u,v)\in X_{\mathcal{D}}\times
Y_{\mathcal{D}}$ be a solution of the gradient scheme (3.5). Then there exists
a constant $C_{2}\geq 0$ depending only on $C_{\rm
ini}>\max(||\Pi_{\mathcal{D}}u^{(0)}||_{L^{2}(\Omega)},||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(0)}||_{L^{2}(\Omega)})$,
such that

$||\Pi_{\mathcal{D}}u||_{L^{\infty}(0,T;L^{2}(\Omega))}+||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v||_{L^{\infty}(0,T;L^{2}(\Omega))}+||\nabla_{\mathcal{D}}u||_{L^{2}(\Omega\times(0,T))^{d}}\leq
C_{2}.$ (5.1)

###### Proof.

Setting, as test functions in Scheme (3.5), the functions $\varphi:=\delta
t^{(n+\frac{1}{2})}u^{(n+1)}$ and $\psi:=\delta t^{(n+\frac{1}{2})}v^{(n+1)}$
leads to

$\displaystyle\int_{\Omega}\Big{(}\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})$
$\displaystyle-\Pi_{\mathcal{D}}u^{(n)}(\boldsymbol{x})\Big{)}\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}+\mu\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad=\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}f_{\mathcal{D}}^{(n+1)}\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$

and

$\displaystyle\int_{\Omega}\Big{(}\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x})-\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n)}(\boldsymbol{x})\Big{)}\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}=\displaystyle\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}g_{\mathcal{D}}^{(n+1)}\Pi_{\mathcal{D}}v^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Apply the inequality $(y-z)y\geq\frac{1}{2}(|y|^{2}-|z|^{2})$ to
$y=\Pi_{\mathcal{D}}u^{(n+1)}$ and $z=\Pi_{\mathcal{D}}u^{(n)}$ (resp.
$y=\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}$ and
$z=\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n)}$) to obtain

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})|^{2}-|\Pi_{\mathcal{D}}u^{(n)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle+\mu\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}f_{\mathcal{D}}^{(n+1)}\Pi_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$

and

$\frac{1}{2}\displaystyle\int_{\Omega}\Big{[}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x})|^{2}-|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}\leq\displaystyle\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}g_{\mathcal{D}}^{(n+1)}\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Sum on $n=0,...,m-1$, for some $m=0,...,N$:

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}-|\Pi_{\mathcal{D}}u^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle+\mu\int_{0}^{t^{(m)}}\int_{\Omega}|\nabla_{\mathcal{D}}u(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\int_{0}^{t^{(m)}}\int_{\Omega}f_{\mathcal{D}}^{(m)}\Pi_{\mathcal{D}}u^{(m)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$

and

$\frac{1}{2}\displaystyle\int_{\Omega}\Big{[}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}(\boldsymbol{x})|^{2}-|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}\leq\displaystyle\int_{0}^{t^{(m)}}\int_{\Omega}g_{\mathcal{D}}^{(m)}\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Applying the Cauchy–Schwarz inequality to the right–hand side of both
inequalities leads to

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}$
$\displaystyle-|\Pi_{\mathcal{D}}u^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}+\mu\int_{0}^{t^{(m)}}\int_{\Omega}|\nabla_{\mathcal{D}}u(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle\quad\leq||f_{\mathcal{D}}^{(m)}||_{L^{2}(\Omega\times(0,T))}\;||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))},$

and

$\frac{1}{2}\displaystyle\int_{\Omega}\Big{[}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}(\boldsymbol{x})|^{2}-|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}\leq||g_{\mathcal{D}}^{(m)}||_{L^{2}(\Omega\times(0,T))}\;||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}.$

Thanks to the assumptions on $f$ and $g$ given in Assumptions 2.1, one writes

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}$
$\displaystyle-|\Pi_{\mathcal{D}}u^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}+\mu\int_{0}^{t^{(m)}}\int_{\Omega}|\nabla_{\mathcal{D}}u(\boldsymbol{x},t)|^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle\leq
c_{1}||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}+c_{2}||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2}$
$\displaystyle+c_{3}|\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))},$

and

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}(\boldsymbol{x})|^{2}$
$\displaystyle-|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\leq
c_{4}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}+c_{6}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2}$
$\displaystyle+c_{5}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}.$

Apply the Young’s inequality, with $\varepsilon$ satisfying
$M_{1}+\frac{1}{2\varepsilon}-\frac{1}{2}>0$ and
$M_{2}+\frac{1}{2\varepsilon}-\frac{1}{2}>0$, to the right–hand side of both
inequalities to obtain

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}$
$\displaystyle-|\Pi_{\mathcal{D}}u^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}+\mu\int_{0}^{t^{(m)}}\int_{\Omega}|\nabla_{\mathcal{D}}u(\boldsymbol{x},t)|^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle\quad\leq\frac{c_{1}^{2}}{2\varepsilon}+M_{1}||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2}+\frac{1}{2\varepsilon}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2}$

and

$\displaystyle\frac{1}{2}\int_{\Omega}\Big{[}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}(\boldsymbol{x})|^{2}$
$\displaystyle-|\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(0)}(\boldsymbol{x})|^{2}\Big{]}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\frac{c_{4}^{2}}{2\varepsilon}+M_{2}||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2}+\frac{1}{2\varepsilon}||\Pi_{\mathcal{D}}u^{(m)}||_{L^{2}(\Omega\times(0,T))}^{2},$

where $M_{1}:=\frac{\varepsilon}{2}+c_{2}+\frac{c_{3}^{2}}{2}$ and
$M_{2}:=\frac{\varepsilon}{2}+c_{6}+\frac{1}{2\varepsilon}$.

Combine the above inequalities together and take the supremum on $m=0,...,N$
establishes Estimate (5.1), since
$\sup_{m=0,...,N}\int_{\Omega}|\Pi_{\mathcal{D}}U^{(m)}|^{2}\,\mathrm{d}\boldsymbol{x}=||\Pi_{\mathcal{D}}U^{(m)}||_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}$
(resp.
$\sup_{m=0,...,N}\int_{\Omega}|\Pi_{{\mathcal{D}}^{{}^{\prime}}}U^{(m)}|^{2}\,\mathrm{d}\boldsymbol{x}=||\Pi_{{\mathcal{D}}^{{}^{\prime}}}U^{(m)}||_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}$).
∎

######  Corollary 5.2.

Assume (2.1) and let ${\mathcal{D}}$ be a gradient discretisation. Then there
exists at least one solution $(u,v)$ to the gradient scheme (3.5).

###### Proof.

At each time step $n+1$, (3.5) provides square non linear equations on
$u^{(n+1)}$ and $v^{(n+1)}$. For a given $w=(w_{1},w_{2})\in
X_{\mathcal{D}}\times Y_{\mathcal{D}}$, $(u,v)\in X_{\mathcal{D}}\times
Y_{\mathcal{D}}$ is is the solution to

$\displaystyle\int_{\Omega}\Pi_{{\mathcal{D}}}\frac{u^{(n+1)}-u^{(n)}}{\delta
t^{(n+\frac{1}{2})}}(x)\Pi_{\mathcal{D}}\varphi(\boldsymbol{x})+\mu\int_{\Omega}\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\cdot\nabla_{\mathcal{D}}\varphi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
(5.2)
$\displaystyle\qquad=\int_{\Omega}f(\Pi_{\mathcal{D}}w_{1},\Pi_{{\mathcal{D}}^{{}^{\prime}}}w_{2})\Pi_{\mathcal{D}}\varphi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x},\quad\forall\varphi\in
X_{\mathcal{D}},$
$\displaystyle\int_{\Omega}\Pi_{{\mathcal{D}}^{{}^{\prime}}}\frac{v^{(n+1)}-v^{(n)}}{\delta
t^{(n+\frac{1}{2})}}\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi(\boldsymbol{x})=\int_{\Omega}g(\Pi_{\mathcal{D}}w_{1},\Pi_{{\mathcal{D}}^{{}^{\prime}}}w_{2})\Pi_{\mathcal{D}}\psi(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x},\quad\forall\psi\in
Y_{\mathcal{D}}.$

This problem describes a linear square system, whose a right hand-side is
built from the terms
$\int_{\Omega}f(\Pi_{\mathcal{D}}w_{1},\Pi_{{\mathcal{D}}^{{}^{\prime}}}w_{2})\Pi_{\mathcal{D}}\varphi\,\mathrm{d}\boldsymbol{x}$,
$\int_{\Omega}\Pi_{\mathcal{D}}u^{(n)}\Pi_{\mathcal{D}}\varphi\,\mathrm{d}\boldsymbol{x}$,
$\int_{\Omega}g(\Pi_{\mathcal{D}}w_{1},\Pi_{{\mathcal{D}}^{{}^{\prime}}}w_{2})\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi\,\mathrm{d}\boldsymbol{x}$,
and
$\int_{\Omega}\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n)}\Pi_{{\mathcal{D}}^{{}^{\prime}}}\psi\,\mathrm{d}\boldsymbol{x}$.

Using arguments similar to the proof of Lemma (5.3), we can obtain

$\displaystyle||\Pi_{\mathcal{D}}u^{(n+1)}||_{L^{2}(\Omega)}+||\nabla_{\mathcal{D}}u^{(n+1)}||_{L^{2}(\Omega)^{d}}\leq
C_{4}||f||_{L^{2}(\Omega)}+||\Pi_{{\mathcal{D}}}u^{(n)}||_{L^{2}(\Omega)},$
$\displaystyle||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n+1)}||_{L^{2}(\Omega)}\leq
C_{5}||g||_{L^{2}(\Omega)}+||\Pi_{{\mathcal{D}}^{{}^{\prime}}}v^{(n)}||_{L^{2}(\Omega)}.$

where $C_{4}$ and $C_{5}$ not depending on $u^{(n+1)}$ or $v^{(n+1)}$. Hence,
the kernel of the matrix associated with the above linear square system is
reduced to $\\{0\\}$, and the matrix is invertible. We then can define the
mapping $T:X_{\mathcal{D}}\times Y_{\mathcal{D}}\to X_{\mathcal{D}}\times
Y_{\mathcal{D}}$ by $T(w)=(u,v)$ with $(u,v)$ is the solution to (5.2). Since
$T$ is continuous, Brouwer’s fixed point establishes the existence of a
solution $u^{(n+1)}$ to the system at time step $n+1$. ∎

######  Lemma 5.3 (Estimates on the discrete gradient and time derivative).

Under Assumptions (2.1), let ${\mathcal{D}}$ be a gradient discretisation in
the sense of Definition 3.1 and let $(u,v)$ be a solution of the gradient
scheme (3.5). If $||\nabla_{\mathcal{D}}u^{(0)}||_{L^{2}(\Omega)}$ is bounded
by $C_{6}$, then there exists a constant $C_{7}$ only depending on $C_{2}$ and
$C_{6}$, such that

$||\delta_{\mathcal{D}}u||_{L^{2}(\Omega\times(0,T))}+||\nabla_{\mathcal{D}}u||_{L^{\infty}(0,T;L^{2}(\Omega)^{d})}\leq
C_{7}$ (5.3)

###### Proof.

Setting, as a test function in Scheme (3.5), the function
$\varphi=u^{(n+1)}-u^{(n)}$ leads to

$\displaystyle\delta
t^{(n+\frac{1}{2})}\int_{\Omega}|\delta_{\mathcal{D}}^{(n+\frac{1}{2})}u(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle+\mu\int_{\Omega}\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})\cdot\nabla_{\mathcal{D}}(u^{(n+1)}(\boldsymbol{x})-u^{(n)}(\boldsymbol{x}))\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad=\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}f_{\mathcal{D}}^{(n+1)}\delta_{\mathcal{D}}^{(n+\frac{1}{2})}u(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Applying the relation $(r-q)\cdot r\geq\frac{1}{2}|r|^{2}-\frac{1}{2}|q|^{2}$
to the second term on the left hand side yields

$\displaystyle\delta
t^{(n+\frac{1}{2})}\int_{\Omega}|\delta_{\mathcal{D}}^{(n+\frac{1}{2})}u(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle+\frac{\mu}{2}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(n+1)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}-\frac{\mu}{2}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(n)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\int_{t^{(n)}}^{t^{(n+1)}}\int_{\Omega}f_{\mathcal{D}}^{(n+1)}\delta_{\mathcal{D}}^{(n+\frac{1}{2})}u(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Sum this inequality on $n=0,...,m-1$, for some $m=0,...,N$ to get

$\displaystyle\int_{0}^{t^{(m)}}\int_{\Omega}|\delta_{\mathcal{D}}u(t)|^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle+\frac{\mu}{2}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}-\frac{\mu}{2}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(0)}|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\int_{0}^{t^{(m)}}\int_{\Omega}f_{\mathcal{D}}^{(m)}\delta_{\mathcal{D}}u(t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Thanks to Assumptions 2.1, Cauchy–Schwars inequality and Young’s inequalities,
we have

$\displaystyle\int_{0}^{t^{(m)}}\int_{\Omega}|\delta_{\mathcal{D}}u(t)|^{2}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle+\frac{\mu}{2}\int_{\Omega}|\nabla_{\mathcal{D}}u^{(m)}(\boldsymbol{x})|^{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad\leq\frac{1}{2}\Big{(}c_{1}+c_{2}||\Pi_{\mathcal{D}}u||_{L^{2}(\Omega\times(0,t^{(m)}))}+c_{3}||\Pi_{\mathcal{D}}v||_{L^{2}(\Omega\times(0,t^{(m)}))}\Big{)}^{2}$
$\displaystyle\qquad+\frac{1}{2}||\delta_{\mathcal{D}}u||_{L^{2}(\Omega\times(0,t^{(m)}))}^{2}+\frac{\mu}{2}C_{6}.$

The proof is concluded by taking the supremum on $m=0,...,N$, using the
inequality $\sup_{n}(a_{n}+b_{n})\leq\sup_{n}(a_{n})+\sup{{}_{n}}(b_{n})$ and
invoking Lemma (5.1) to estimate the second and third terms on the right hand
side. ∎

### Proof of Theorem 3.6

The proof follows the compactness technique detailed in [17] and it is split
into three steps:

Step 1: Compactness results.
Due to Estimate (5.1) and to the consistency and the limit–conformity of GD,
[19, Lemma 4.8] provides $\bar{u}\in L^{2}(0,T;H^{1}(\Omega))$ and $\bar{v}\in
L^{2}(0,T;L^{2}(\Omega))$, such that, up to a subsequence,
$\Pi_{{\mathcal{D}}_{m}}u_{m}\to\bar{u}$ weakly in $L^{2}(0,T,L^{2}(\Omega))$,
$\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}\to\bar{v}$ weakly in
$L^{2}(0,T,L^{2}(\Omega))$, and
$\nabla_{{\mathcal{D}}_{m}}u_{m}\to\nabla\bar{u}$ weakly in
$L^{2}(0,T,L^{2}(\Omega)^{d})$, as $m\to\infty$. Since
$||\nabla_{\mathcal{D}}u^{(0)}||_{L^{2}(\Omega)}$ is assumed to bounded,
Estimate (5.3) holds. This estimate with the three properties ( consistency,
limit–conformity and compactness) show that the assumptions of [19, Theorem
4.18] are satisfied. This theorem proves that $\Pi_{{\mathcal{D}}_{m}}u_{m}$
converges strongly to $\bar{u}$ in $L^{\infty}(0,T;L^{2}(\Omega))$, as
$m\to\infty$. In fact, $\partial_{t}\bar{u}\in L^{2}(0,T;L^{2}(\Omega))$ and
$\delta_{{\mathcal{D}}_{m}}u_{m}$ converges weakly to $\partial_{t}\bar{u}$ in
$L^{2}(0,T;L^{2}(\Omega))$, as $m\to\infty$. Estimates (5.1) yields the weak
convergence of $\Pi_{{\mathcal{D}}_{m}}v_{m}$ to $\bar{v}$ in
$L^{\infty}(0,T;L^{2}(\Omega))$, as $m\to\infty$.

Step 2: $(\bar{u},\bar{v})$ is a solution to the continuous problem.
The a.e. convergence of $\Pi_{{\mathcal{D}}_{m}}u_{m}$, the assumptions on
$f_{1}$ and $f_{2}$ and the dominated convergence theorem show that
$f_{1}(\Pi_{{\mathcal{D}}_{m}}u_{m})\to f_{1}(\bar{u})$ and
$f_{2}(\Pi_{{\mathcal{D}}_{m}}u_{m})\to f_{2}(\bar{u})$ in
$L^{2}(\Omega\times(0,T))$. Consider $\bar{\varphi}\in
L^{2}(0,T;H^{1}(\Omega))$. [19, Lemma 4.10] shows the existence of
$w_{m}=(w_{m}^{(n)})_{n=0,...,N_{m}}\in X_{{\mathcal{D}}_{m}}^{N_{m}+1}$, such
that $\Pi_{{\mathcal{D}}_{m}}w_{m}\to\bar{\varphi}$ in
$L^{2}(0,T;L^{2}(\Omega))$ and
$\nabla_{{\mathcal{D}}_{m}}w_{m}\to\nabla\bar{\varphi}$ in
$L^{2}(0,T;L^{2}(\Omega)^{d})$. Take $\varphi=\delta
t_{m}^{(n+\frac{1}{2})}w_{m}^{(n)}$ as a test function in (3.5a), sum on
$n=0,...,N_{m}-1$, and pass to the limit $m\to\infty$ to see that
$(\bar{u},\bar{v})$ satisfies (2.8a), thanks to the strong–weak convergence of
$\Pi_{{\mathcal{D}}_{m}}u_{m}$, $\Pi_{{\mathcal{D}}_{m}^{,}}v_{m}$ and
$\delta_{{\mathcal{D}}_{m}}u_{m}$.

Let us now verify (2.8b). Take $\bar{\psi}\in L^{2}(0,T;L^{2}(\Omega))$ such
that $\partial_{t}\bar{\psi}\in L^{2}(\Omega\times(0,T))$ and
$\bar{\psi}(T,\cdot)=0$. By [19, Lemma 4.10] (with a slight change), we can
find $w_{m}=(w_{m}^{(n)})_{n=0,...,N_{m}}\in Y_{{\mathcal{D}}_{m}}^{N_{m}+1}$,
such that $\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}\to\bar{\psi}$ in
$L^{2}(0,T;L^{2}(\Omega))$ and
$\delta_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}\to\partial_{t}\bar{\psi}$
strongly in $L^{2}(\Omega\times(0,T))$. Set $\psi=\delta
t_{m}^{(n+\frac{1}{2})}w_{m}^{(n)}$ as a test function in (3.5b) and sum on
$n=0,...,N_{m}-1$ to obtain

$\displaystyle\sum_{n=0}^{N_{m}-1}\int_{\Omega}\left(\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}^{(n+1)}(\boldsymbol{x})-\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v^{(n)}(\boldsymbol{x})\right)$
$\displaystyle\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}^{(n)}(\boldsymbol{x})\,\mathrm{d}x$
(5.4)
$\displaystyle=\int_{0}^{T}\int_{\Omega}g_{{\mathcal{D}}_{m}}^{(n+1)}\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}^{(n)}(\boldsymbol{x})\,\mathrm{d}x\,\mathrm{d}t.$

For two families of real numbers $(a_{n})_{n=0,...,{\mathbb{N}}}$ and
$(b_{n})_{n=0,...,N}$, the form of discrete integration by part [19, Eq.
(D.15)] is

$\displaystyle\sum_{n=0}^{N-1}\Big{(}b_{n+1}-b_{n}\Big{)}a_{n}=b_{N}a_{N}-b_{0}a_{0}-\displaystyle\sum_{n=0}^{N-1}b_{n+1}\Big{(}a_{n+1}-a_{n}\Big{)}.$

Apply this relation to the right hand side in (5.4) and use the fact
$w^{(N)}=0$ to deduce

$\displaystyle-\int_{0}^{T}\int_{\Omega}\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}(\boldsymbol{x},t)\delta_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle-\int_{\Omega}\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}^{(0)}(\boldsymbol{x})\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}^{(0)}(\boldsymbol{x})\,\mathrm{d}\boldsymbol{x}$
$\displaystyle\quad=\int_{0}^{T}\int_{\Omega}g_{{\mathcal{D}}_{m}}^{(n+1)}\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}w_{m}^{(n)}(\boldsymbol{x})\,\mathrm{d}x\,\mathrm{d}t.$

By the space–time consistency,
$\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{m}^{(0)}=\Pi_{{\mathcal{D}}_{m}^{{}^{\prime}}}J_{{\mathcal{D}}_{m}^{{}^{\prime}}}v_{\rm
ini}\to v_{\rm ini}$ in $L^{2}(\Omega)$. The strong convergence of
$\Pi_{{\mathcal{D}}_{m}}u_{m}$, the assumptions on $g_{1}$ and the dominated
convergence theorem lead to $g_{1}(\Pi_{{\mathcal{D}}_{m}}u_{m})\to
g_{1}(\bar{u})$ in $L^{2}(\Omega\times(0,T))$. Hence, passing to the limit
$m\to\infty$ in the above equation implies that $\bar{v}$ satisfies (3.5b).

Step 3: $\nabla_{{\mathcal{D}}_{m}}u_{m}$ converges strongly.
Taking $\bar{\varphi}:=u_{m}$ in (3.5a) and passing to the superior limit
gives (choose $\varphi=\bar{u}$ in (2.8a))

$\displaystyle\lim\sup_{m\to\infty}\int_{0}^{T}\int_{\Omega}\nabla_{{\mathcal{D}}_{m}}$
$\displaystyle
u_{m}(\boldsymbol{x},t)\cdot\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)\,\mathrm{d}x\,\mathrm{d}t$
(5.5)
$\displaystyle=\int_{0}^{T}\int_{\Omega}f(\bar{u},\bar{v})\bar{u}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t-\int_{0}^{T}\int_{\Omega}\partial_{t}\bar{u}(\boldsymbol{x},t)\bar{u}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle=\int_{0}^{T}\int_{\Omega}\nabla\bar{u}(\boldsymbol{x},t)\cdot\nabla\bar{u}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

One can write

$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)-\nabla\bar{u}(\boldsymbol{x},t)\right)\cdot(\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)-\nabla\bar{u}(\boldsymbol{x},t))\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle\quad=\int_{0}^{T}\int_{\Omega}\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)\cdot\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t-\int_{0}^{T}\int_{\Omega}\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)\cdot\nabla\bar{u}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle\qquad-\int_{0}^{T}\int_{\Omega}\nabla\bar{u}(\boldsymbol{x},t)\cdot(\nabla_{{\mathcal{D}}_{m}}u_{m}(\boldsymbol{x},t)-\nabla\bar{u}(\boldsymbol{x},t))\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$

Thanks to the weak convergence of $\nabla_{{\mathcal{D}}_{m}}u_{m}$ and (5.5),
passing to the limit $m\to\infty$ in each of the terms above completes the
proof.

##  6\. Numerical results

We examine here the validity of the HMM method in solving the reaction-
diffusion models by exploring the dynamics of the spiral waves of the Barkley
model. This model is considered as an example of the reaction-diffusion model
(2.1)–(2.5) with reaction terms taking the form [6]:

$\displaystyle f(\bar{u},\bar{v})$
$\displaystyle=\frac{1}{\rho}\bar{u}(1-\bar{u})(\bar{u}-\frac{\bar{v}+b}{a}),$
(6.1) $\displaystyle g(\bar{u},\bar{v})$ $\displaystyle=\bar{u}-\bar{v},$

for $a,b>0$. The small parameter $0<\rho\ll 1$ represents the time scale
separation of fast variable $\bar{u}$ and slow variable $\bar{v}$. Before
proceeding to do that, we first review the dynamics of the underlying model to
understand the mechanism of nucleation and propagation of spiral waves in
excitable media.

###  6.1. Dynamics of the Barkley model

A pair of essential features are shared by all reaction-diffusion models:
spatially localized excitation that undergoes diffusion in space and quickly
reverts to a recovery state. Precisely, if the stimulus is strong enough, an
excitable system will switch from a state of rest to a state of excitement,
then quickly revert to a refractory state and back to rest. A new excitement
cannot be created until a certain period (refractory time) has gone by. So an
excitable system has the ability to support the wave propagation caused by
strong disturbances when in a rest state driven by local nonlinearity coupled
with diffusion [10, 6]. Figure 1 [4] shows the dynamics of the Barkley model
or reaction kinetics when diffusion is not present, with nullclines pictures
of $\bar{u}$ and $\bar{v}$. The $\bar{u}$-nullclines ($f(\bar{u},\bar{v})=0$)
are represented by three straight lines:
$\bar{u}=0,\bar{u}=1,\bar{u}=\frac{\bar{v}+b}{a}$, whereas
$\bar{v}$-nullclines ($g(\bar{u},\bar{v})=0$) is the line $\bar{u}=\bar{v}$,
where $a,b>0$. Excitable media systems are characterized by dynamic states of
excitation and recovery. Thus, by setting a small boundary, say $\delta$,
bordering the line $\bar{u}=0$, a given point $(\bar{u},\bar{v})$ is said to
be excited if $\bar{u}>\delta$ and recovering otherwise.

Figure 6.1. Illustration of local dynamics of Barkley model in absence of
diffusion ($\mu=0$) taken from [4].

The detail of local dynamics is given by the physical parameters $\rho$, $a$
and $b$. $\rho$ is selected very small, so that the dynamics of the activator
$\bar{u}$ is much faster than the inhibitor $\bar{v}$ inside the region of
excitement. However, $\bar{u}\approx 0$ within recovery region and therefore
the exponentially decay of the inhibitor $\bar{v}$ affects only the local
dynamics [6].

Increasing $a$ would cause a lengthening of the excitation period and
increasing $\frac{b}{a}$ would raise the excitation threshold [11].
Intersection of all nullclines yields the fixed points $(0,0)$ and $(1,1)$.
The origin $(0,0)$ is stable and excitable fixed point of the model with
excitation threshold $u_{th}=\frac{\bar{v}+b}{a}$.

To be precise, the initial conditions found on the left of the threshold
$u_{th}$ near the $(0,0)$ decay straight to the fixed point of origin.
Conversely, the initial conditions found on the right of $u_{th}$ take a large
excursion prior to shrinking down to $(0,0)$, see [6] for more details. When
spatial diffusion is added to the reaction kinetics we gain spiral wave
dynamics [11]. This allows us to examine the ways in which spiral waves are
created, their movement, the way they collapse and disappear, and this offers
us insights into the appearance, modification, and ways of stopping
arrhythmias.

###  6.2. Spiral waves

The wave fronts and backs of spiral waves meet to form a tip which pivots
around the region named the core. If the core of the spiral wave moves then it
is regarded as being drifting. In order to examine the ways in which spiral
waves behave, employing the basis of the HMM method, we consider a pair of
examples for the Barkley model at the edge of the simulation domains that have
no-flux boundary conditions. In both cases, the HMM scheme is tasted on a
square domain $\Omega=[-L,L]^{2}$, which dividing into a uniform mesh of
$3584$ triangular elements.

In all test cases, we consider the Barkley model defined in the introduction
of this section (6.1). Letting $\boldsymbol{x}=(x_{1},x_{2})$, the initial
conditions, required to initiate the spiral waves, are chosen in general form
by

$\displaystyle\bar{u}_{0}(x_{1},x_{2})$
$\displaystyle=(1+\exp(4(|x_{1}|-\alpha_{1})))^{-2}-(1+\exp(4(|x_{1}|-\alpha_{2})))^{-2},$
$\displaystyle\bar{v}_{0}(x_{1},x_{2})$ $\displaystyle=\alpha_{3},$ (6.2)

where $\alpha_{1}$,$\alpha_{2}$ and $\alpha_{3}$ are real constants. For the
time–discretisation, the forward Euler method is considered with time step
$\delta t>0$.

### Example 1

For the first example, the parameters of the model are chosen according to
[37] as $\rho=0.005$, $a=0.3$, $b=0.01$, $\mu=1$ and $L=30$. The initial
conditions are given as

$u_{\rm ini}(x_{1},x_{2})=\begin{cases}\hfill 0&,\quad\mbox{for }x_{1}<0\mbox{
or }x_{2}>5,\\\ \bar{u}_{0}&,\quad\mbox{otherwise},\end{cases}$

and

$v_{\rm ini}(x_{1},x_{2})=\begin{cases}\hfill\bar{v}_{0}&,\quad\mbox{for
}x_{1}<1\mbox{ and }x_{2}<10,\\\ 0&,\quad\mbox{otherwise},\end{cases}$

where $\bar{u}_{0}$ and $\bar{v}_{0}$ are defined in Equation (6.2), by
setting $\alpha_{1}=5$, $\alpha_{2}=1$ and $\alpha_{3}=0.1$.

|  |
---|---|---
Time: $t=5$ |  Time: $t=10$ |  Time: $t=20$
|  |
Time: $t=30$ |  Time: $t=40$ |  Time: $t=60$
|  |
Time: $t=80$ |  Time: $t=100$ |
Figure 6.2. Example 1, propagation of spiral waves of the Barkley model at
different time levels with parameters detailed in text.

Figure 6.2 shows the way that spiral waves were propagated in this instance,
plotted at different times, with the simulation is done up to $t=100$. It can
be seen that the development of the spiral wave involves a core at the
medium’s center, and this develops into a periodically rotating single spiral
wave. The trajectory of this spiral wave runs symmetrically from origin to
boundary, moving along the length to turn into stationary spirals. The
observed behavior is called the reflection of the spiral at the boundary.

### Example 2

However, in certain instances, the spiral tip may collide with the boundary
and then the spiral waves are annihilated instead of being reflected,
especially if the initial displacement of the tip is close to the boundary as
can be seen in the second example. In this case, we choose the parameters of
the model as [37]: $\rho=0.0208$, $a=0.52$, $b=0.05$ and $\mu=1$ on the square
domain, which is of the length $L=7.5$. The initial conditions are given as

$u_{\rm ini}(x_{1},x_{2})=\begin{cases}\hfill 0&,\quad\mbox{for }x_{1}<0\mbox{
or }x_{2}>5,\\\ \bar{u}_{0}&,\quad\mbox{otherwise},\end{cases}$

and

$v_{\rm ini}(x_{1},x_{2})=\begin{cases}\hfill\bar{v}_{0}&,\quad\mbox{for
}x_{1}<-1\mbox{ and }x_{2}<3,\\\ 0&,\quad\mbox{otherwise},\end{cases}$

where $\bar{u}_{0}$ and $\bar{v}_{0}$ are as in Equation (6.2), by putting
$\alpha_{1}=3$, $\alpha_{2}=1$ and $\alpha_{3}=0.25$. Spiral solution for the
Barkley model under the considered parameters is shown in Figure 6.3 at
different time levels, from $t=0.2$ to $t=3$. In this instance, with a small
domain, the spiral commences taking shape at $t=1$, although the rotation is
not yet complete. Following this, each part of the spiral will drift over to
the nearest boundary and be absorbed.

A dormant state will be achieved in both instances as when there are
interactions between spiral waves and boundaries within excitable media, the
trajectory of the wave is either reflected or destroyed by the interaction
with the boundary, according to the parameters that have been chosen. We
conclude that our numerical observations for the Barkley model obtained using
HMM method are consistent with those observations in the previous studies [37,
9] in terms of spiral waves dynamics.

|  |
---|---|---
Time: $t=5$ |  Time: $t=10$ |  Time: $t=20$
|  |
Time: $t=30$ |  Time: $t=40$ |  Time: $t=60$
|  |
Time: $t=80$ |  Time: $t=100$ |
Figure 6.3. Example 2, propagation of spiral waves of the Barkley model at
different time levels with parameters detailed in text.

## References

  * [1] R. R. Aliev and A. V. Panfilov, A simple two-variable model of cardiac excitation, Chaos Solitons Fractals, 7 (1996), pp. 293–301.
  * [2] Y. Alnashri and J. Droniou, Gradient schemes for the Signorini and the obstacle problems, and application to hybrid mimetic mixed methods, Computers and Mathematics with Applications, 72 (2016), pp. 2788–2807.
  * [3] , A gradient discretization method to analyze numerical schemes for nonlinear variational inequalities, application to the seepage problem, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2375–2405.
  * [4] H. Alzubaidi and T. Shardlow, Interaction of waves in a one dimensional stochastic pde model of excitable media, Discrete And Continuous Dynamical Systems Series B, 18 (2013), pp. 1735–1754.
  * [5] F. Amdjadi and J. Gomatam, Spiral waves on static and moving spherical domains, Journal of Computational and Applied Mathematics, 182 (2005), pp. 472–486.
  * [6] D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D, 49 (1991), pp. 61–70.
  * [7] M. Bendahmane, R. Bùrger, and R. Ruiz-Baier, A finite volume scheme for cardiac propagation in media with isotropic conductivities, Mathematics and Computers in Simulation, 80 (2010), pp. 1821–1840.
  * [8] M. Bendahmane and K. H. Karlsen, Convergence of a finite volume scheme for the bidomain model of cardiac tissue, Applied Numerical Mathematics, 59 (2009), pp. 2266–2284.
  * [9] Y. Bourgault, M. Ethier, and V. LeBlanc, Simulation of electrophysiological waves with an unstructured finite element method, Mathematical Modelling and Numerical Analysis, 37 (2003), pp. 649–661.
  * [10] N. Britton, Reaction-Diffusion Equations and Their Applications to Biology, Academic Press, New York, London, 1986.
  * [11] R. Bürger, R. Ruiz-Baier, and K. Schneider, Adaptive multiresolution methods for the simulation of waves in excitable media, Journal of Scientific Computing, 43 (2010), pp. 261–290.
  * [12] B. H. C. Diks and J. Degoede, Spiral wave dynamics, Chaos Solitom and Fractals, 5 (1995), pp. 645–660.
  * [13] D. J. Christini and L. Glass, Introduction: Mapping and control of complex cardiac arrhythmias, Chaos, 12 (2002), pp. 732–739.
  * [14] Y. Coudière and C. Pierre, Stability and convergence of a finite volume method for two systems of reaction-diffusion equations in electro-cardiology, Nonlinear Analysis: Real World Applications, 7 (2006), pp. 916 – 935.
  * [15] Y. Coudière and R. Turpault, Very high order finite volume methods for cardiac electrophysiology, Computers and Mathematics with Applications, 74 (2017), pp. 684–700.
  * [16] Devanand and J. C. Kalita, Hoc simulation of barkley model in excitable media, AIP Conference Proceedings, 1975 (2018), p. 030011.
  * [17] J. Droniou, Introduction to discrete functional analysis techniques for the numerical study of diffusion equations with irregular data, in Proceedings of the 17th Biennial Computational Techniques and Applications Conference, CTAC-2014, J. Sharples and J. Bunder, eds., vol. 56 of ANZIAM J., 2015, pp. C101–C127.
  * [18] J. Droniou and R. Eymard, Uniform-in-time convergence of numerical methods for non-linear degenerate parabolic equations, Numer. Math., 132 (2016), pp. 721–766.
  * [19] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin, The gradient discretisation method, Mathematics & Applications, Springer, Heidelberg, 2018.
  * [20] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin, A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methodss, Mathematical Models and Methods in Applied Sciences, 20 (2010), pp. 265–295.
  * [21] J. Droniou, R. Eymard, and R. Herbin, Gradient schemes: a generic framework for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic problems, Mathematical Models and Methods in Applied Sciences, 23 (2013), pp. 2395–2432.
  * [22] J. Droniou, R. Eymard, and R. Herbin, Gradient schemes: generic tools for the numerical analysis of diffusion equations, M2AN Math. Model. Numer. Anal., 50 (2016), pp. 749–781.
  * [23] R. Eymard, P. Feron, T. Gallouët, R. Herbin, and C. Guichard, Gradient schemes for the Stefan problem, International Journal On Finite Volumes, 10s (2013).
  * [24] R. Eymard, C. Guichard, and R. Herbin, Small-stencil 3D schemes for diffusive flows in porous media, ESAIM. Mathematical Modelling and Numerical Analysis, 46 (2012), pp. 265–290.
  * [25] R. Eymard, C. Guichard, R. Herbin, and R. Masson, Gradient schemes for two-phase flow in heterogeneous porous media and Richards equation, ZAMM Z. Angew. Math. Mech., 94 (2014), pp. 560–585.
  * [26] F. H. Fentona, E. M. Cherry, H. M. Hastings, and S. J. Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity, Chaos, 12 (2002), pp. 852–892.
  * [27] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophysical Journal, 1 (1961), pp. 445–466.
  * [28] M. Gomez-Gesteira, A. M. nuzuri, V. P.-M. nuzuri, and V. Perez-Villar, Boundary imposed spiral drift, Physical Review E, 53 (1996), pp. 5480–5483.
  * [29] A. Goryachev and R. Kapral, Spiral waves in chaotic systems, Physical Review Letters, 76 (1996), pp. 1619–1622.
  * [30] D. Harrild and C. Henriquez, A finite volume model of cardiac propagation, Annals of Biomedical Engineering, 25 (1997), pp. 315–334.
  * [31] A. Karma, Meandering transition in two-dimensional excitable media, Physical Review Letters, 65 (1990), pp. 2824–2828.
  * [32] J. Keener, A geometrical theory for spiral waves in excitable media, SIAM Journal on Applied Mathematics, 46 (1986), pp. 1039–1059.
  * [33] , Arrhythmias by dimension, Proceedings of Symposia in Applied Mathematics, 59 (2002), pp. 57–81.
  * [34] J. P. Keener and J. Sneyd, Mathematical Physiology, Springer-Verlag, 1998.
  * [35] J. Li and J. Li, High-order compact difference methods for simulating wave propogation in excitable media, International Journal Of Numerical Analysis And Modeling, Series B, 5 (2014), pp. 339–346.
  * [36] D. Olmos and B. Shizgal, Annihilation and reflection of spiral waves at a boundary for the beeler-reuter model, Physical Review E, 77 (2008), pp. 33–50.
  * [37] , Pseudospectral method of solution of the fitzhugh nagumo equation, Mathematics and Computers in Simulation, 79 (2009), pp. 2258–2278.
  * [38] J. Ramos, Spiral wave break-up and planar front formation in two-dimensional reactive diffusive media with straining, Applied Mathematics and Computation, 154 (2004), pp. 697–711.
  * [39] J. Rogers, M. Courtemanche, and A. McCulloch, Finite element methods for modelling impulse propagation in the heart., in Computational biology of the heart. Workshop on Whole Heart Modeling, 11th–13th February 1994, Utrecht, the Netherlands, Chichester: Wiley, 1997, pp. 217–233.
  * [40] J. M. Rogers and A. D. McCulloch, A collocation-galerkin finite element model of cardiac action potential propagation, IEEE Transactions on Biomedical Engineering, 41 (1994), pp. 743–757.
  * [41] J. Smoller, Shock Waves and Reaction–Diffusion Equations, Springer, 1983.
  * [42] M. Trew, I. L. Grice, B. Smaill, and A. Pullan, A finite volume method for modeling discontinuous electrical activation in cardiac tissue, Annals of Biomedical Engineering, 33 (2005), pp. 590–602.
  * [43] A. Winfree, Varieties of spiral wave behavior: an experimentalist’s approach to the theory of excitable media, Chaos, 1 (1991), pp. 303–334.
  * [44] Y. A. Yermakova and A. M. Pertsov, Interaction of rotating spiral waves with a boundary, Biophys, 31 (1986), pp. 932–940.
  * [45] A. M. Zhabotinskya, A history of chemical oscillations and waves, CHAOS, 1 (1991), pp. 379–386.