A posteriori error estimates for the monodomain model in cardiac electrophysiology

TL;DR

This paper develops and validates a posteriori error estimators for the monodomain model in cardiac electrophysiology, enabling accurate assessment of numerical solution errors from discretization and linearization.

Contribution

It introduces reliable and efficient a posteriori error estimators for the monodomain model, accounting for multiple error sources and supported by numerical validation.

Findings

Error estimators are reliable and efficient under certain conditions.

Numerical experiments confirm the theoretical error bounds.

The approach improves error control in cardiac electrophysiology simulations.

Abstract

We consider the monodomain model, a system of a parabolic semilinear reaction-diffusion equation coupled with a nonlinear ordinary differential equation, arising from the (simplified) mathematical description of the electrical activity of the heart. We derive a posteriori error estimators accounting for different sources of error (space/time discretization and linearization). We prove reliability and efficiency (this latter under a suitable assumption) of the error indicators. Finally, numerical experiments assess the validity of the theoretical results.

Tables

Table 1. Table 1: Values of the main parameters of the model

$M$	$A$	$ϵ$	$a$
$1$	$8$	$0.2$	$0.15$

Equations

$$\{\begin{array}{rlr}{\partial }_{t}u-\nabla \cdot (M\nabla u)+f(u,w)& =0& \text{in }\Omega \times (0,T),\\ M{\partial }_{\nu }u& =0& \text{on }\partial \Omega \times (0,T),\\ u{|}_{t=0}& =u_0& \text{in }\Omega ,\\ {\partial }_{t}w+g(u,w)& =0& \text{in }\Omega \times (0,T),\\ w{|}_{t=0}& =w_0& \text{in }\Omega ,\end{array}$$

$$f(u,w)=Au(u-a)(u-1)+uw,g(u,w)=ϵ(Au(u-1-a)+w),$$

$$0\le u(x,t)\le 1,0\le w(x,t)\le \frac{A(1+a{)}^2}{4}\forall (x,t)\in \Omega \times (0,T).$$

$$\{\begin{array}{rl}{\int }_{\Omega }{\partial }_{t}u\phi dx+{\int }_{\Omega }M\nabla u\cdot \nabla \phi dx+{\int }_{\Omega }f(u,w)\phi dx& =0\forall \phi \in H^1(\Omega ),\\ {\int }_{\Omega }{\partial }_{t}w\psi dx+{\int }_{\Omega }g(u,w)\psi dx& =0\forall \psi \in L^2(\Omega ).\end{array}$$

$$X(t_a,t_b)$$

$$Y(t_a,t_b)$$

$${\Vert u\Vert }_{X(t_a,t_b)}$$

$${\Vert w\Vert }_{Y(t_a,t_b)}$$

$$\{t_n{\}}_{n=0}^N\subset [0,T];t_0=0,t_N=T;t_n-t_{n-1}={\tau }_n>0,$$

$$u^0=u_0;w^0=w_0;$$

$${\int }_{\Omega }\frac{u^n-u^{n-1}}{{\tau }^n}\phi dx+{\int }_{\Omega }M\nabla u^n\cdot \nabla \phi dx+{\int }_{\Omega }f(u^n,w^n)\phi dx=0$$

$${\int }_{\Omega }\frac{w^n-w^{n-1}}{{\tau }^n}\psi dx+{\int }_{\Omega }g(u^n,w^n)\psi dx=0$$

$$⟨{\mathcal{F}}^1(u,w),\phi ⟩$$

$$⟨{\mathcal{F}}^2(u,w),\psi ⟩$$

$$\left[\begin{array}{cc}{\mathcal{F}}_u^1(u_{k-1}^n,w_{k-1}^n)& {\mathcal{F}}_w^1(u_{k-1}^n,w_{k-1}^n)\\ {\mathcal{F}}_u^2(u_{k-1}^n,w_{k-1}^n)& {\mathcal{F}}_w^2(u_{k-1}^n,w_{k-1}^n)\end{array}\right]\left[\begin{array}{c}\delta u\\ \delta w\end{array}\right]=\left[\begin{array}{c}-{\mathcal{F}}^1(u_{k-1}^n,w_{k-1}^n)\\ -{\mathcal{F}}^2(u_{k-1}^n,w_{k-1}^n)\end{array}\right]\text{ in }{H}^{\ast }\times L^2(\Omega );$$

$${\int }_{\Omega }\frac{u_k^n}{{\tau }^n}\phi dx+{\int }_{\Omega }M\nabla u_k^n\cdot \nabla \phi dx+{\int }_{\Omega }\left[f_u(u_{k-1}^n,w_{k-1}^n)u_k^n+f_w(u_{k-1}^n,w_{k-1}^n)w_k^n\right]\phi dx={\int }_{\Omega }\frac{u^{n-1}}{{\tau }^n}\phi dx$$

$$+{\int }_{\Omega }\left[f_u(u_{k-1}^n,w_{k-1}^n)u_{k-1}^n+f_w(u_{k-1}^n,w_{k-1}^n)w_{k-1}^n-f(u_{k-1}^n,w_{k-1}^n)\right]\phi dx\forall \phi \in H^1(\Omega )$$

$${\int }_{\Omega }\frac{w_k^n}{{\tau }^n}\psi dx+{\int }_{\Omega }\left[g_u(u_{k-1}^n,u_{k-1}^n)u_k^n+g_w(u_{k-1}^n,u_{k-1}^n)w_k^n\right]\psi dx={\int }_{\Omega }\frac{w^{n-1}}{{\tau }^n}\psi dx$$

$$+{\int }_{\Omega }\left[g_u(u_{k-1}^n,w_{k-1}^n)u_{k-1}^n+g_w(u_{k-1}^n,w_{k-1}^n)w_{k-1}^n-g(u_{k-1}^n,w_{k-1}^n)\right]\psi dx\forall \psi \in L^2(\Omega ).$$

$$\rho (K^{\prime },K)=\left\{\frac{h_{K^{\prime }}}{h_K},K^{\prime }\in {\mathcal{T}}_h^n,K\in {\tilde{\mathcal{T}}}_h^n:K\subset K^{\prime }\right\},$$

$$V_h^n=\{v_h\in C(\bar{\Omega }),v_h{|}_K\in {\mathbb{P}}_1(K)\forall K\in {\tilde{\mathcal{T}}}_h^n\}$$

$$\{\begin{array}{rl}{\int }_{\Omega }{\partial }_t{u}_h{\phi }_{h}dx+{\int }_{\Omega }M\nabla u_h\cdot \nabla {\phi }_{h}dx+{\int }_{\Omega }f(u_h,w_h){\phi }_{h}dx& =0\forall {\phi }_h\in V_h,\\ {\int }_{\Omega }{\partial }_t{w}_h{\psi }_{h}dx+{\int }_{\Omega }g(u_h,w_h){\psi }_{h}dx& =0\forall {\psi }_h\in V_h.\end{array}$$

$${\left\{{\Vert u-u_h\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2+{\Vert w-w_h\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2\right\}}^{\frac{1}{2}}\le c{h}^2.$$

$${\left\{{\Vert u-u_h\Vert }_{X(0,T)}^2+{\Vert w-w_h\Vert }_{Y(0,T)}^2\right\}}^{\frac{1}{2}}\le ch.$$

$$\frac{1}{2}\frac{d}{dt}$$

$$+{\int }_{\Omega }{\partial }_t(w-w_h)(w-{\psi }_h)dx+{\int }_{\Omega }\nabla (u-u_h)\cdot \nabla (u-{\phi }_h)dx$$

$$+{\int }_{\Omega }(f(u,w)-f(u_h,w_h)(u_h-{\phi }_h)dx+{\int }_{\Omega }(g(u,w)-g(u_h,w_h)(w_h-{\psi }_h)dx.$$

$${\int }_{\Omega }{\partial }_t(u-u_h)(u-{\Pi }_{h}u)dx$$

$$=\frac{1}{2}\frac{d}{dt}{\Vert u-{\Pi }_{h}u\Vert }_{L^2(\Omega )}^2.$$

$$\frac{1}{2}\frac{d}{dt}$$

$$\frac{1}{2}\frac{d}{dt}\left({\Vert u-{\Pi }_{h}u\Vert }_{L^2(\Omega )}^2+{\Vert w-{\Pi }_{h}w\Vert }_{L^2(\Omega )}^2\right)+\frac{1}{2}{\Vert \nabla (u-{\Pi }_{h}u)\Vert }_{L^2(\Omega )}^2+|E_{{\Pi }_h}|,$$

$$E_{{\Pi }_h}={\int }_{\Omega }(f(u,w)-f(u_h,w_h))(u_h-{\Pi }_{h}u)dx+{\int }_{\Omega }(g(u,w)-g(u_h,w_h))(w_h-{\Pi }_{h}w)dx.$$

$${\Vert u(t)-u_h(t)\Vert }_{L^2(\Omega )}^2+{\Vert w(t)-w_h(t)\Vert }_{L^2(\Omega )}^2+{\int }_0^t{\Vert \nabla (u(s)-u_h(s))\Vert }_{L^2(\Omega )}^{2}ds$$

$$\le {\Vert u_0-{\Pi }_h{u}_0\Vert }_{L^2(\Omega )}^2+{\Vert w_0-{\Pi }_h{w}_0\Vert }_{L^2(\Omega )}^2+{\Vert u(t)-{\Pi }_{h}u(t)\Vert }_{L^2(\Omega )}^2+{\Vert w(t)-{\Pi }_{h}w(t)\Vert }_{L^2(\Omega )}^2$$

$$+{\int }_0^t{\Vert \nabla (u(s)-{\Pi }_{h}u(s))\Vert }_{L^2(\Omega )}^{2}ds+2{\int }_0^t|E_{{\Pi }_h}(s)|ds.$$

$${\Vert u-u_h\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2+{\Vert w-w_h\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2\le {\Vert u_0-{\Pi }_h{u}_0\Vert }_{L^2(\Omega )}^2+{\Vert w_0-{\Pi }_h{w}_0\Vert }_{L^2(\Omega )}^2$$

$$+{\Vert u-{\Pi }_{h}u\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2+{\Vert w-{\Pi }_{h}w\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2+{\Vert u-{\Pi }_{h}u\Vert }_{L^2(0,T;H^1(\Omega ))}^2$$

Full text

1 Introduction

The main goal of this paper is the a posteriori numerical analysis of the monodomain model, a system of a parabolic semilinear reaction-diffusion equation coupled with a nonlinear ordinary differential equation, arising from the mathematical description of the electrical activity of the heart. The monodomain model represents a simplified version of the more realistic bidomain model which has been object in recent years of an intense research activity, see e.g. [9] and references therein. For the purpose of the paper, we first recall [12], where a careful a priori analysis of the Galerkin semidiscrete space approximation of the bidomain system is performed, investigating convergence properties and stability estimates for the semidiscrete solution. This result, coupled with the argument regarding the time-discretization analysis provided in [10], allows for an exhaustive a priori error analysis for the bidomain model. Moreover, in [8] the authors introduce a space-time adaptive algorithm for the solution of the bidomain model by resorting to a stepsize control for the temporal adaptivity, whereas spatial adaptivity is performed by virtue of a posteriori local error estimators. However, a complete a posteriori error analysis is missing.

With the aim of contributing to fill this gap, in this paper we focus on the simpler monodomain model and provide a detailed a posteriori analysis. In particular, we consider a Newton-Galerkin approximation of the monodomain system and look for a posteriori indicators of the error involving the $L^{2}(0,T;H^{1}(\Omega))$ norm. Inspired by the seminal work [17] and by the recent papers [11, 1], we derive a posteriori error bounds by providing a suitable splitting of the total residual into three operators, accounting for different sources of error entailed by the discretization process. Specifically, we introduce a linearization residual, a time discretization residual, and a space discretization residual, with the additional difficulty with respect, e.g., to [1] represented by the coupled structure of the system of differential equations.

The a posteriori analysis is complemented with an a priori analysis which relies on previous results obtained in [13], where error estimates with respect to the $L^{\infty}(0,T;L^{2}(\Omega))$ norm of the error are obtained. Here, we derive a priori estimates for the semidiscrete problem in a different norm involving the $L^{2}(0,T;H^{1}(\Omega))$ one.

The a posteriori error estimators obtained in this paper can be employed to derive fully space-time adaptive algorithms that can be of particular importance, for instance, in the solution of inverse problems like the identification of ischemic regions (i.e. areas in which the coefficient of the system are altered from the reference values) by means of boundary voltage. An iterative algorithm (as the one proposed in [3] for a simplified model) would greatly benefit from an adaptive approach that would drastically reduce the computational cost.

The paper is organized as follows: in Section 2 we introduce the Newton-Galerkin full discretization of the monodomain model, whereas Section 3 is devoted to the a priori estimates for the problem. In Section 4 we introduce the residual operators associated to the discrete solution and prove the equivalence between the error and the residual (in suitable norms). In Section 5 we define three a posteriori estimators and employ them to prove an upper bound for the approximation error. We also provide, under a suitable assumption, a lower estimate for the error in terms of the same indicators, thus assessing their efficiency. Finally, Section 6 reports some numerical experiments assessing the validity of the derived estimates and investigating convergence rates both of the error and of the estimators as the discretization parameters are reduced.

2 A Newton-Galerkin scheme for the approximation of the monodomain model

Let $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, be an open bounded domain. Consider the monodomain model (see [9, 14])

[TABLE]

being $u$ the trasmembrane electrical potential in the cardiac tissue and $M:\Omega\rightarrow\mathbb{R}^{d\times d}$ the conductivity tensor. In particular, according to the biological application, we assume that $M$ is constant in time, and in each point $x\in\Omega$ the tensor $M(x)$ is a symmetric positive definite matrix, with positive eigenvalues $\mu_{i}$, $i=1,\ldots,d$. Moreover, we suppose that $\mu_{i}(x)$ are uniform in space and denote by $\mu_{min}$ and $\mu_{max}$ the minimum and the maximum eigenvalue, respectively. The associated eigenvectors may instead vary in space, and we assume that the overall matrix function $M(x)$ is smooth. The nonlinear term $f(u,w)$ models the current induced by the motion of ions across the membrane, and is addressed as ionic current. According to a well established phenomenological approach (see, e.g., [14]), $f$ is a function of the potential $u$ and of a recovery variable $w$, whose dynamics is governed by a coupled nonlinear ordinary differential equation involving a nonlinear term $g$. We focus in particular on the Aliev-Panfilov model of the cardiac tissue, according to the version reported, e.g., in [4]; namely, the nonlinear terms $f$ and $g$ are as follows:

[TABLE]

with $A,\epsilon>0$, $0<a<1$. Such a problem is showed to be well-posed: in particular, we refer to [2], which extends the results contained in [13] to the model of interest, and guarantees the following existence, uniqueness and comparison result:

Proposition 2.1**.**

Let the initial data $u_{0}\in C^{2+\alpha}(\bar{\Omega})$, $w_{0}\in C^{\alpha}(\bar{\Omega})$ satisfy the bound $0\leq u_{0}\leq 1$ and $0\leq w_{0}\leq\frac{A(1+a)^{2}}{4}$, consider $M\in C^{2}(\Omega)$ and let the following compatibility conditions hold: $M\nabla u_{0}\cdot\nu=0$, being $\partial\Omega\in C^{2+\alpha}$. Then, there exists a unique classical solution $(u,w)$ of (2.1), $u\in C^{2+\alpha,1+\alpha/2}(\overline{\Omega}\times[0,T])$ and $w\in C^{\alpha,1+\alpha/2}(\overline{\Omega}\times[0,T])$. Moreover, it holds that

[TABLE]

Remark 2.1**.**

When considering $w_{0}\in C^{2+\alpha}(\bar{\Omega})$, one can easily conclude (see [2]) that also $w\in C^{2+\alpha,1+\alpha/2}(\overline{\Omega}\times[0,T])$. In particular, both $u(\cdot,t)$ and $w(\cdot,t)$ belong to the Sobolev’s space $H^{2}(\Omega)$ for each $t\in[0,T]$.

Remark 2.2**.**

If the conductivity tensor $M$ is only in $L^{\infty}(\Omega)$ (as in the case of an ischemic heart), we can nevertheless show ([2]) the existence and uniqueness of the weak solution $(u,w)$ s.t. $u\in L^{2}(0,T;H^{1}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}u\in L^{2}(0,T;H^{*})$, $w\in L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}w\in L^{2}((0,T)\times\Omega)$, where $H^{*}=(H^{1}(\Omega))^{*}$. Moreover, the same bounds on $u,w$ hold as above and it is possible to guarantee additional regularity on the solution, namely $u\in C^{\alpha,\alpha/2}(\overline{\Omega}\times[0,T])$, $w\in C^{\alpha,1+\alpha/2}(\overline{\Omega}\times[0,T])$.

The weak formulation of (2.1) reads

[TABLE]

For each time interval $(t_{a},t_{b})\subset(0,T)$, we introduce the following functional spaces:

[TABLE]

which are Banach spaces endowed with the norms:

[TABLE]

To ease the notation, we denote with $X$ and $Y$ the spaces $X(0,T)$ and $Y(0,T)$, respectively.

We now introduce a time semidiscretization of the problem by employing an implicit Euler scheme: consider a partition of the time interval

[TABLE]

and define the semidiscrete solution as the couple of collections $(\{u^{n}\}_{n=0}^{N},\{w^{n}\}_{n=0}^{N})$, being $u^{n}\in H^{1}(\Omega)$, $w^{n}\in L^{2}(\Omega)$ with $n=0,\ldots,N$ such that

[TABLE]

Consider the operators $\mathcal{F}^{1}:H^{1}(\Omega)\times L^{2}(\Omega)\rightarrow(H^{1}(\Omega))^{*}$, $\mathcal{F}^{2}:H^{1}(\Omega)\times L^{2}(\Omega)\rightarrow L^{2}(\Omega)$, which are defined interval-wise as follows: for $t\in(t_{n-1},t_{n}]$

[TABLE]

The functionals $\mathcal{F}^{1}$ and $\mathcal{F}^{2}$ are (Fréchet) differentiable with respect to the $H^{1}(\Omega)$ norm in the variable $u$ and with respect to the $L^{2}(\Omega)$ norm in the variable $w$, respectively. This allows to define a Newton scheme for the solution of the nonlinear system (2.5)-(2.6) as follows:

Computing the expression of the derivatives of $\mathcal{F}^{1}$ and $\mathcal{F}^{2}$, and substituting $\delta u=u_{k}^{n}-u_{k-1}^{n}$, $\delta w=w_{k}^{n}-w_{k-1}^{n}$, the system (2.7) can be rewritten as

[TABLE]

[TABLE]

Following [17], we introduce an affinely equivalent, admissible, and shape-regular tessellation $\mathcal{T}_{h}^{n}$ for each instant $t_{n}$. For each element $K$ of $\mathcal{T}_{h}^{n}$, we denote by $h_{K}$ its diameter, and require $h_{K}\leq h$. We moreover require the following conditions to hold:

1. i)

$\forall n\leq 1$, there exists a common refinement $\widetilde{\mathcal{T}}_{h}^{n}$ of both $\mathcal{T}_{h}^{n}$ and $\mathcal{T}_{h}^{n-1}$; 2. ii)

$\exists\rho_{*},\rho^{*}>0$ independent of $n$ and $h$ s.t., defined

[TABLE]

then $\rho_{*}\leq\rho(K^{\prime},K)\leq\rho^{*}$ $\forall K\in\mathcal{T}_{h}^{n}$, $\forall n=1,\cdots,N$;

Taking advantage of $\widetilde{\mathcal{T}}_{h}^{n}$, we introduce the Finite Element discrete space $V_{h}^{n}\subset H^{1}(\Omega)$

[TABLE]

and the $L^{2}$ orthogonal projection $\Pi_{h}^{n}:L^{2}(\Omega)\rightarrow V_{h}^{n}$.

The fully discrete solution of (2.1) consists in the pair of collections $(\{u_{h,k}^{n}\},\{u_{h,k}^{n}\})$, with $n=0,\ldots,N$ and $k=0,\ldots,K_{n}$, being $K_{n}$ the maximum number of iterations performed in each timestep (possibly varying with $n$). In particular, $\{u_{h,K_{n}}^{n}\}$ and $\{w_{h,K_{n}}^{n}\}$ are such that:

• •

$u_{h}^{0}=\Pi_{h}^{0}u_{0}$, $w_{h}^{0}=\Pi_{h}^{0}w_{0}$, the projections of the initial data on $\widetilde{\mathcal{T}}_{h}^{0}=\mathcal{T}_{h}^{0}$ ;

• •

for each $n=1,\cdots,N$, we initialize the Newton algorithm with $u_{h,0}^{n}=\Pi_{h}^{n}u_{h,K_{n-1}}^{n-1}$;

• •

for each $n=1,\cdots,N$, for each $k=1,\cdots,K_{n}$, the couple $(u_{h,k}^{n},w_{h,k}^{n})\in V_{h}^{n}\times V_{h}^{n}$ solves the system (2.8)-(2.9) for all $(\varphi_{h},\psi_{h})\in V_{h}^{n}\times V_{h}^{n}$.

3 A priori estimates for the space semidiscretization

In this section we consider a priori error estimates for the space semidiscretized problem under the assumption that the same tessellation $\mathcal{T}_{h}$ is considered in each instant, together with the discrete space $V_{h}$ of linear finite elements. We refer to the space semidiscrete solution as to the couple of functions $(u_{h},w_{h}):[0,T]\rightarrow V_{h}\times V_{h}$ satisfying $u_{h}(0)=u_{h}^{0}$, $w_{h}(0)=w_{h}^{0}$ and

[TABLE]

Taking advantage of standard inverse estimates and approximation results (see [5]), it is possible to prove the following result:

Theorem 3.1**.**

There exists a unique solution $(u_{h},w_{h})$ of problem (3.1) in $(C^{1}(0,T;V_{h}))^{2}$. Moreover, for any fixed $h_{0}$ there exists a positive $\delta_{0}$ such that $(u_{h},w_{h})\in S_{\delta_{0}}$ $\forall x,t\in\Omega\times[0,T]$, being $S_{\delta_{0}}=[-\delta_{0},1+\delta_{0}]\times[-\delta_{0},\frac{A(1+a)^{2}}{4}+\delta_{0}]$. Finally, there exists a constant $c$ depending on $u,w,u_{0},w_{0},f,g,\Omega,T$ and independent of $h$ such that

[TABLE]

The proof of this theorem relies on techniques introduced in [16]: with minor modifications, it is possible to adapt the proof of [13, Theorem 4.4] to the present context where the Aliev-Panfilov electrophysiological model is considered. We are moreover interested in establishing the convergence rate of the $X$ and $Y$ norms of the error. This is the object of the following result:

Theorem 3.2**.**

There exists a constant $c$ depending on $u,w,u_{0},w_{0},f,g,\Omega,T$ and independent of $h$ such that

[TABLE]

Proof.

Consider the equations of system (2.3), test them with the functions $u_{h}-\varphi_{h}$ and $w_{h}-\psi_{h}$ respectively, being $\varphi_{h}$, $\psi_{h}\in V_{h}$, and sum them. Repeating the same procedure on system (3.1) and subracting the two equations obtained, we get

[TABLE]

Consider now $\varphi_{h}=\Pi_{h}u$ and $\psi_{h}=\Pi_{h}w$, being $\Pi_{h}$ the $L^{2}$ orthogonal projection on $V_{h}$ operator, and observe that

[TABLE]

A similar result hold for the second term on the rihgt-hand side of (3.4). By Cauchy-Schwarz and Young inequalities, we conclude that

[TABLE]

being

[TABLE]

Integrating from [math] to $t$, and employing the fundamental theorem of calculus, together with the choice $u_{h}(0)=u_{h}^{0}=\Pi_{h}u_{0}$, $w_{h}(0)=w_{h}^{0}=\Pi_{h}w_{0}$, we get

[TABLE]

It immediately follows that

[TABLE]

Now, we observe that, since both $(u,w)$ and $(u_{h},w_{h})$ belong to $S_{\delta_{0}}$ for a suitable value of $\delta_{0}$ (see Theorem 3.2) and since the functions $f,g$ are Lipschitz continuous on $S_{\delta_{0}}$ with constants bounded by $c_{\delta_{0}}>0$, it holds

[TABLE]

In conclusion, we have

[TABLE]

Applying standard approximation properties of $V_{h}$, taking advantage of the fact that both $u(\cdot,t)$ and $w(\cdot,t)$ belong to $H^{2}(\Omega)$ for $t\in[0,T]$ (see Remark 2.2), we can conclude that the following suboptimal estimate holds:

[TABLE]

In view of this estimate, from (3.5) we infer that ${\left\|u-u_{h}\right\|}_{L^{2}(0,T;H^{1}(\Omega))}\leq ch$.

To conclude, we need to consider the terms involving the derivative in time. This requires the introduction of the elliptic projection operator associated to the bilinear form $\int_{\Omega}M\nabla u\cdot\nabla v+\mu_{min}\int_{\Omega}uv$, i.e., the map $R_{h}:H^{1}(\Omega)\rightarrow V_{h}$ such that

[TABLE]

According to the properties of $R_{h}$ (see, e.g., [15]), we know that $\forall u\in H^{1}(\Omega)$ it holds

[TABLE]

By employing the first equation in system (3.1), for each $\varphi\in H^{1}(\Omega)$ it holds

[TABLE]

According to (2.3), and in view of (3.6), we can conclude that $\forall\varphi\in H^{1}(\Omega)$

[TABLE]

Via Cauchy-Schwarz inequality we obtain

[TABLE]

Now, we show that ${\left\|\partial_{t}u_{h}\right\|}_{L^{2}((0,T)\times\Omega)}$ is bounded by a constant independent of $h$. Indeed, considering $\varphi_{h}=\partial_{t}u_{h}$ in the first equation of (3.1), we obtain

[TABLE]

Integrating from [math] to $T$, we get

[TABLE]

Thus, it holds that

[TABLE]

and by solving the second-order inequality, we conclude that

[TABLE]

In view of (3.8), employing (3.7), the above estimate for ${\left\|u-u_{h}\right\|}_{L^{2}(0,T;H^{1}(\Omega))}$, together with the estimates for ${\left\|u-u_{h}\right\|}_{L^{\infty}(0,T;L^{2}(\Omega))}$ and ${\left\|w-w_{h}\right\|}_{L^{\infty}(0,T;L^{2}(\Omega))}$ in Theorem 3.2, we get

[TABLE]

An analogous argument holds for ${\left\|\partial_{t}(w-w_{h})\right\|}_{L^{2}((0,T)\times\Omega)}$, and the thesis follows. ∎

Remark 3.1**.**

When stating the discrete problem (3.1), we have neglected any error introduced by the computation of the integral $\int_{\Omega}M\nabla u_{h}\cdot\nabla w_{h}dx$. When $M$ is a polynomial function, the integration can be performed exactly by choosing a suitable quadrature rule. In case $M$ is not a polynomial but still sufficiently smoot (e.g., $M\in C^{1+\alpha}(\Omega)$), the quadrature error do not affect the results contained in Theorems 3.2 and 3.3, as can be verified by an application of Strang’s lemma. When considering the case of a piecewse smooth coefficient $M$ (which occurs, e.g., when modeling an ischemic cardiac tissue), one should adopt a different strategy, as suggested, e.g., in [6].

4 Residual operators

We now move towards the introduction of a posteriori estimators. Consider the fully discrete solution $(\{u_{h,k}^{n}\},\{w_{h,k}^{n}\})$ as introduced in Section 1, being again $\{\mathcal{T}_{h}^{n}\}_{n}$ possibly different tessellations among the different discrete instants. Collecting all the final indices $K_{n}$ in a multi-index $\mathbf{k}=[K_{n}]_{n=1}^{N}$, the associated linear interpolated solution $(u_{h,\tau}^{(\mathbf{k})},w_{h,\tau}^{(\mathbf{k})})$ is a couple of continuous functions on $[0,T]$, defined timestep-wise as follows: for each $t\in(t_{n-1},t_{n}]$, $n=1,\ldots,N$,

[TABLE]

We now define for almost each instant $t$ the residual operator $R(t)$ in the product space $(H^{1}(\Omega)\times L^{2}(\Omega))^{*}=H^{*}\times L^{2}(\Omega)$, being $H^{*}$ the dual space of $H^{1}(\Omega)$:

[TABLE]

It is now possible to prove a result of equivalence between the $X,Y$ norms of the error and the dual norms of the residual operators. More precisely, it holds:

Theorem 4.1**.**

The functions ${\left\|R_{1}(t)\right\|}_{H^{*}}$ and ${\left\|R_{2}(t)\right\|}_{L^{2}(\Omega)}$ are square integrable on each interval $(t_{a},t_{b})\subset(0,T)$, and moreover

[TABLE]

where $c_{*}$ and $c^{*}$ depend on $\Omega,\mu_{max},\mu_{min},f,g$ and $T$.

Proof.

By employing equation (2.3) together with the expressions of $R_{1}(t)$ and $R_{2}(t)$ we have, $\forall\varphi\in H^{1}(\Omega)$, $\forall\psi\in L^{2}(\Omega)$, $a.e.$ $t\in(0,T)$

[TABLE]

Fixing $\psi=0$ and employing the Cauchy-Schwarz inequality and the fact that $f$ is Lipschitz continuous with constant $K_{f}$,

[TABLE]

Thus, computing the $L^{2}$ norm on $(t_{a},t_{b})$ we obtain

[TABLE]

Analogously, when taking $\varphi=0$, we get

[TABLE]

being $K_{g}$ the Lipschitz constant of $g$. Summing (4.5) and (4.6) we obtain (4.3a).

To prove (4.3b), consider (4.4) and take $\varphi=u-u_{h,\tau}^{(\mathbf{k})}$, $\psi=w-w_{h,\tau}^{(\mathbf{k})}$; by mean value theorem 222 Applied to the real valued function $\textit{h}:\zeta\in\mathbb{R}\rightarrow h(\zeta)=\int_{\Omega}f(u+\zeta(u_{h,\tau}^{(\mathbf{k})}-u),w+\zeta(w_{h,\tau}^{(\mathbf{k})}-w))(u_{h,\tau}^{(\mathbf{k})}-u)dx$, it guarantees that there exists $\zeta^{*}\in[0,1]$ $s.t.$ $\textit{h}(1)-\textit{h}(0)=h^{\prime}(\zeta*)$. it holds that $\exists(\xi_{1},\eta_{1})=(u,w)+\zeta_{1}(u_{h,\tau}^{(\mathbf{k})}-u,w_{h,\tau}^{(\mathbf{k})}-1),(\xi_{2},\eta_{2})=(u,w)+\zeta_{2}(u_{h,\tau}^{(\mathbf{k})}-u,w_{h,\tau}^{(\mathbf{k})}-1)$, such that

[TABLE]

Consider now the quadratic form $\mathcal{Q}:H^{1}(\Omega)\times L^{2}(\Omega)\rightarrow\mathbb{R}$,

[TABLE]

which allows to rewrite the previous equation as

[TABLE]

It clearly holds that $|\mathcal{Q}(m,n)|\leq\lambda_{max}({\left\|m\right\|}_{L^{2}(\Omega)}+{\left\|n\right\|}_{L^{2}(\Omega)})$, being $\lambda_{max}$ a continuous function of $f_{u}(\xi_{1},\eta_{1}),f_{w}(\xi_{1},\eta_{1}),g_{u}(\xi_{2},\eta_{2}),g_{w}(\xi_{2},\eta_{2})$. Hence, $\lambda_{max}$ depends both on $x$ and $t$, but thanks to a priori bounds on $(\xi_{1},\eta_{1})$ and $(\xi_{2},\eta_{2})$ (inherited from Proposition 2.1 and Theorem 3.2), we can ensure it is bounded from above on $\Omega\times(0,T)$ by a positive constant $\Lambda$. Via Cauchy-Schwarz and Young inequalities,

[TABLE]

hence

[TABLE]

Let us now take a fixed $t\in(0,T)$ and integrate from [math] to $t$, obtaining

[TABLE]

Via Gronwall’s inequality, we obtain

[TABLE]

whence the bound on ${\left\|u-u_{h,\tau}^{(\mathbf{k})}\right\|}_{L^{\infty}(0,t,L^{2}(\Omega))}$ and ${\left\|w-w_{h,\tau}^{(\mathbf{k})}\right\|}_{L^{\infty}(0,t,L^{2}(\Omega))}$. Moreover, from (4.7) we get

[TABLE]

Finally, taking $\psi=0$ in (4.4), by Cauchy-Schwarz inequality we get

[TABLE]

thus

[TABLE]

A similar strategy allows to conclude that an analogous bound holds for $\partial_{t}(w-w_{h,\tau}^{(\mathbf{k})})$, hence every part of the norms ${\left\|u-u_{h,\tau}^{(\mathbf{k})}\right\|}_{X(0,t)}$, ${\left\|w-w_{h,\tau}^{(\mathbf{k})}\right\|}_{Y(0,t)}$ is bounded as in the thesis. ∎

According to the strategy proposed in [1], it is now possible to perform a decomposition of the residual operators, by distinguishing the contribution from space discretization, time discretization and linearization as follows :

[TABLE]

[TABLE]

It is immediate to verify that $R_{1}(t)=R_{1}^{h}(t)+R_{1}^{\tau}(t)+R_{1}^{k}(t)$ in $H^{*}$ and $R_{2}(t)=R_{2}^{h}(t)+R_{2}^{\tau}(t)+R_{2}^{k}(t)$ in $L^{2}(\Omega)$; moreover, in view of the discrete problem (2.8)-(2.9), the following orthogonality property holds:

[TABLE]

5 A posteriori estimators

We denote by $\widetilde{\mathcal{E}}_{h}^{n}$ the set of all faces of $\mathcal{T}_{h}^{n}$ and distinguish between the set of boundary faces $\widetilde{\mathcal{E}}_{h,\partial\Omega}^{n}$ and the set $\widetilde{\mathcal{E}}_{h,int}^{n}$ of the interior ones. Each face $E\in\widetilde{\mathcal{E}}_{h,int}^{n}$ is shared by two distinct elements, which we denote as $K_{E,1}$ and $K_{E,2}$; we define the jump of the conormal derivative across $E$ as

[TABLE]

where $\nu_{E,1}$ and $\nu_{E,2}$ are outer the normals of $E$ with respect to $K_{E,1}$ and $K_{E,2}$, hence $\nu_{E,1}=-\nu_{E,2}$. For each face $E$ of $\widetilde{\mathcal{E}}_{h,\partial\Omega}^{n}$ (which belongs to a single element $K$ of the tessellation), we set

[TABLE]

We now introduce the following computable quantities which will appear in the a posteriori estimates:

Space indicators

[TABLE]

Time indicators

[TABLE]

Linearization indicators

[TABLE]

The first main result of this section is the following a posteriori upper bound:

Theorem 5.1**.**

For each discrete solution $(\{u_{h,k}^{n}\},\{w_{h,k}^{n}\}$ with $n=1,\ldots,N$, $k=1,\ldots,K_{n}$, collecting all $K_{n}$ in the multi-index $\mathbf{k}=[K_{n}]_{n=1}^{N}$ and definining $u_{h,\tau}^{(\mathbf{k})},w_{h,\tau}^{(\mathbf{k})}$ as in (4.1), it holds that for each $n=1,\ldots,N$:

[TABLE]

where the symbol $\lesssim$ denotes that an inequality holds up to a positive multiplicative constant independent of the space discretization parameter $h$.

In order to prove Theorem 5.1, we need a preliminary results dealing with the spatial residual operators only.

Lemma 5.1**.**

There exist two positive constants $c_{\dagger}$, $c^{{\dagger}}$ independent of $n$ s.t., for almost every $t\in(t_{n-1},t_{n})$ and for each $n=1,\ldots,N$, it holds:

[TABLE]

Proof.

We follow the strategy outlined in [17, Lemma 5.1] (see also [1]). In particular, since $R_{1}^{k}(t)$ and $R_{2}^{k}(t)$ are constant in time within each interval $(t_{n-1},t_{n})$, estimates (5.2) can be proved by similar arguments as the ones employed for elliptic problems. We now consider $t\in(t_{n-1},t_{n})$ and neglect the dependence of $R_{1}^{h}$, $R_{2}^{h}$ on $t$. Integrating by parts the expression of $R_{1}^{h}$ , we obtain that for each $\varphi\in H^{1}(\Omega)$

[TABLE]

We now introduce the Clément interpolation operator $I_{h}:H^{1}(\Omega)\rightarrow V_{h}^{n}$ (see [7], [5]); proceeding in a standard way (see, e.g., [18]) and employing the orthogonality properties in (4.10) we have

[TABLE]

where $\widetilde{\omega}_{K}$ (respectively, $\widetilde{\omega}_{E}$) is the union of all the elements of $\widetilde{\mathcal{T}}_{h}^{n}$ containing at least a vertex of $K$ (respectively, $E$). This entails that

[TABLE]

By an application of the Cauchy-Schwarz inequality it follows that ${\left\|R_{2}^{h}\right\|}_{L^{2}(\Omega)}\leq{\left\|R_{\Omega,2}\right\|}_{L^{2}(\Omega)}$, hence the estimate from above in (5.2) holds with $c^{\dagger}=\left(\max\{1,(C^{\dagger})^{2}\}\right)^{\frac{1}{2}}$.

In order to prove the lower bound, we introduce

[TABLE]

with $\alpha,\beta>0$, $\phi_{K}$, $\phi_{E}$ the baricentrical bubble functions respectively on $K$ and $\omega_{E}=K_{E,1}\cup K_{E,2}$. Analogously to [17, Lemma 5.1], we can show that

[TABLE]

and

[TABLE]

which entails that

[TABLE]

Regarding $R_{2}^{h}$, the following equality clearly holds

[TABLE]

and this, together with (5.3) allows to conclude the lower bound in (5.2) with $\frac{1}{c_{\dagger}}=\left(\min\left\{1,\frac{1}{C_{\dagger}^{2}}\right\}\right)^{\frac{1}{2}}$. ∎

It is now possible to prove the upper bound (5.1).

Proof of Theorem 5.1.

In view of (4.3b), we only need to prove that, for each $n=1,\ldots,N$, it holds

[TABLE]

According to Lemma 5.2,

[TABLE]

and since by definition both $R_{1}^{h}$ and $R_{2}^{h}$ are constant in each interval $(t_{n-1},t_{n})$, we conclude that

[TABLE]

Moreover, it is immediate to verify via Cauchy-Schwarz inequality that

[TABLE]

which, integrating on $(t_{n-1},t_{n})$ yields

[TABLE]

Eventually, again by the Cauchy-Schwarz inequality and employing (4.1), for each $t\in(t_{n-1},t_{n})$

[TABLE]

Since $\int_{t_{n-1}}^{t_{n}}\left(\frac{t_{n}-t}{\tau_{n}}\right)^{2}=\frac{\tau_{n}}{3}$, we get

[TABLE]

By means of the triangular inequality, (5.5), (5.6) and (5.7) we obtain (5.4), and hence (5.1). ∎

5.1 Efficiency of the estimators

The upper estimate provided in (5.1) holds for any choice of $\mathbf{k}$, i.e., the total number of Newton iterations $K_{n}$ performed in each interval $(t_{n-1},t_{n})$ can be selected arbitrarily. We now prove a result of efficiency for our a posteriori estimators, which holds true when a specific condition on the indices $K_{n}$ is satisfied. In particular, for each $n\geq 1$, we assume as in [11, equation (3.12)] that there exists $K_{n}$ such that

[TABLE]

being $\sigma<\frac{1}{c_{\dagger}}$, where $c_{\dagger}$ is the constant appearing in Lemma 5.2. Such an hypothesis can be understood as a stopping cryterion for the Newton algorithm associated to each timestep $n$. In particular, (5.8) prescribes that an iteration $K_{n}$ is considered acceptable if the correspondent computable indicator of the linearization error is sufficiently smaller than the one associated to the space error.

Moreover, we need to introduce the following assumption on the nonlinear terms $f$ and $g$: $\exists\lambda>0$ (without loss of generality, we assume $\lambda\leq\mu_{min}$) such that, $\forall u_{1},u_{2},w_{1},w_{2}\in\mathbb{R}$,

[TABLE]

This assumption is verified under small modifications of the original problem by a large class of models, including Aliev-Panfilov, see Remark 5.1.

Theorem 5.2**.**

Let $f,g$ satisfy (5.9) and let $(\{u_{h,k}^{n}\}\{w_{h,k}^{n}\})$, $n=0,\ldots,N$, $k=0,\ldots,K_{n}$ be the fully discrete solution of (2.1) obtained by the Newton scheme (2.8)-(2.9), satisfying assumption (5.8) on the choice of $K_{n}$. Then,

[TABLE]

being $u_{h,\tau}^{(\mathbf{k})}$, $w_{h,\tau}^{(\mathbf{k})}$ the interpolants defined in (4.1).

Proof.

First of all, we exploit the assumption (5.9) on $f,g$ to obtain a useful inequality. Consider the temporal residual operators $R_{1}^{\tau}$, $R_{2}^{\tau}$ with test functions $\varphi_{1}=u_{h,\tau}^{(\mathbf{k})}-u_{h,K_{n}}^{n}=-\frac{t_{n}-t}{\tau_{n}}(u_{h,K_{n}}^{n}-u_{h,K_{n-1}}^{n-1})$, $\psi_{1}=w_{h,\tau}^{(\mathbf{k})}-w_{h,K_{n}}^{n}=-\frac{t_{n}-t}{\tau_{n}}(w_{h,K_{n}}^{n}-w_{h,K_{n-1}}^{n-1})$:

[TABLE]

We recall that

[TABLE]

when integrating in time, we can bound the right-hand side by considering two terms at a time as follows:

[TABLE]

where we set ${\left\|err\right\|}_{XY}\mathrel{\mathop{:}}=\left({\left\|u-u_{h,\tau}^{(\mathbf{k})}\right\|}_{X(t_{n-1},t_{n})}^{2}+{\left\|w-w_{h,\tau}^{(\mathbf{k})}\right\|}_{Y(t_{n-1},t_{n})}^{2}\right)^{\frac{1}{2}}$ and we made use of (4.3a) and of the Jensen inequality $A+B\leq\sqrt{2}(A^{2}+B^{2})^{\frac{1}{2}}$. Moreover, via (5.2) we get

[TABLE]

Eventually, by the definition of $R_{1}^{k},R_{2}^{k}$ and $\gamma_{k}^{n}$,

[TABLE]

This allows to conclude that

[TABLE]

We focus now on the spatial estimator $\eta_{k}^{n}$. According to the proof of Lemma 5.2, for the particular choice of test functions $\varphi_{2}=W_{n}$, $\psi_{2}=R_{\Omega,2}$, it holds that

[TABLE]

whence

[TABLE]

By the decomposition of the residual, $R_{1}^{h}=R_{1}-R_{1}^{\tau}-R_{1}^{k}$ and $R_{2}^{h}=R_{2}-R_{2}^{\tau}-R_{2}^{k}$. Moreover,

[TABLE]

where $K_{f}$ and $K_{g}$ are the Lipschitz constants of $f$ and $g$ and $K_{fg}=max\{\mu_{max},K_{f},K_{g}\}$. Exploiting the Cauchy-Schwarz and the Jensen inequalities and the definition of $\gamma_{k}^{n}$,

[TABLE]

and since $u_{h,\tau}^{(\mathbf{k})}-u_{h,K_{n}}^{n}=\frac{t_{n}-t}{\tau_{n}}(u_{h,K_{n}}^{n}-u_{h,K_{n-1}}^{n-1})$, we have

[TABLE]

Now, we take advantage of the strategy used in the proof of the lower bound in [17], in particular, choosing a positive $\alpha$, we multiply the inequality (5.12) by $(\alpha+1)\left(\frac{t-t_{n-1}}{\tau_{n}}\right)^{\alpha}$ and integrate from $t_{n-1}$ to $t_{n}$. We observe that

[TABLE]

Thus, we obtain (applying (4.3a) and (5.11))

[TABLE]

Taking advantage of the assumtpion (5.8) and dividing by $\sqrt{\tau_{n}}$, we get

[TABLE]

Since by assumption (5.8) $\frac{1}{c_{\dagger}}-\sigma>0$, selecting

[TABLE]

we can ensure that

[TABLE]

Thus, we deduce

[TABLE]

from now on, we omit the explicit expression of the constants in front of each term in the inequality. As an immediate consequence, again by (5.8), we infer

[TABLE]

We now focus on $\theta_{k}^{n}$. By definition,

[TABLE]

Therefore, in view of (5.11)

[TABLE]

and eventually (using (5.14) and (5.15))

[TABLE]

Eventually, collecting the results (5.14), (5.15), (5.16) we conclude that

[TABLE]

∎

Remark 5.1**.**

Assumption (5.9) is in general not satisfied by $f$ and $g$ as in (2.2). In particular, inequality (5.9) holds with a possibly negative constant, $-\tilde{K}$. This can be deduced by mean value theorem, exploiting the fact that $f,g$ in (2.2) are continuously differentiable and take values on a bounded subset of $\mathbb{R}^{2}$ due to the uniform a priori bounds on the solutions prescribed in Proposition 2.1. However, we can introduce a change of variable in the original problem (2.1): for a positive $\lambda$, we set $\tilde{u}=e^{-(\tilde{K}+\lambda)t}u$ and $\tilde{w}=e^{-(\tilde{K}+\lambda)t}w$. It holds $\partial_{t}\tilde{u}=-(\tilde{K}+\lambda)\tilde{u}+e^{-(\tilde{K}+\lambda)t}\partial_{t}u$, and $(\tilde{u},\tilde{w})$ is the solution of

[TABLE]

where $\tilde{f}=e^{-(\tilde{K}+\lambda)t}f(e^{(\tilde{K}+\lambda)t}\tilde{u},e^{(\tilde{K}+\lambda)t}\tilde{w})+(\tilde{K}+\lambda)\tilde{u}$ and $\tilde{g}$ (analogously defined) satisfy (5.9).

Remark 5.2**.**

In the particular case where the source of error coming from the linearization process is disregarded, the simplified counterpart of Theorem 5.1 holds with the only estimators $\eta^{n}$, $\theta^{n}$ defined as

[TABLE]

being $u_{h,\tau}=\frac{t_{n}-t}{\tau_{n}}u_{h}^{n-1}+\frac{t-t_{n-1}}{\tau_{n}}u_{h}^{n}$ and $w_{h,\tau}=\frac{t_{n}-t}{\tau_{n}}w_{h}^{n-1}+\frac{t-t_{n-1}}{\tau_{n}}w_{h}^{n}$. An efficiency result analogous to Theorem 5.2 holds with the same estimators, clearly without requiring (5.8).

6 Numerical experiments

We now numerically assess the validity of the derived a posteriori estimates. We consider the following two-dimensional setup: the domain $\Omega$ is the square $(0,1)^{2}$, whereas the time interval is set equal to $(0,16)$. All the experiments are performed in an isotropic tissue, whence $M$ is a scalar coefficient. We consider the initial data

[TABLE]

whereas the value of the constants of the problem are reported in Table 1.

We report in Figure 1 several snapshots of the evolution of the electrical potential $u$ throughout time. The results are obtained via the Newton-Galerkin scheme in (2.8)-(2.9), making use of the same computational mesh $\mathcal{T}_{h}$ for each instant, with maximum diameter $h=0.0125$ and a fixed timestep $\tau=0.025$. As an exit criterion for the Newton iterations we check if the distance between two following iterations (measured in $H^{1}$ and $L^{2}$ norm respectively for $u$ and $w$) is below a suitable tolerance, which we set as $tol=10^{-14}$. In accordance with experimental observations (see, e.g., [9]), the nonlinear dynamics shows a first quick propagation of the stimulus in the tissue and, after a plateau phase, a slow decrease of the electrical potential.

6.1 Spatial and temporal analysis

We now verify the validity of the estimates stated in Theorem 5.1. Due to the lack of an analytical expression for the solution of (2.1), we need to build a high-fidelity numerical solution $(\tilde{u},\tilde{w})$. In particular, we employ a reference fine mesh with $h_{ref}=4\cdot 10^{-3}$ and a time step $\tau_{ref}=2\cdot 10^{-3}$ to solve the Newton scheme (2.8)-(2.9), where $tol=10^{-15}$ is employed to make negligible the linearization error (see Remark 5.18). Employing $(\tilde{u},\tilde{w})$ it is possible to compute the error associated to different discrete solutions, obtained with different values of $h$ and $\tau$, and to assess the validity of the a posteriori error estimates introduced in Theorem 5.1 employing in particular the estimators defined in (5.18).

In Figure 2 we report the numerical verification of the upper bound (5.1) for two different choices of the discretization parameters $h$ and $\tau$. Each line is piecewise constant on every interval $(t_{n-1},t_{n})$. The red line represents the norm of the error on the interval $(0,t_{n})$ (see the left-hand side of (5.1) for its precise definition) computed with respect to the high-fidelity solution, whereas the blue line shows the sum of the estimators in each interval until $t_{n}$ (see the left-hand side of (5.1)). In this case the upper bound holds with constant $1$.

Moreover, in Figure 3 we investigate the convergence rates for both the a posteriori estimator and the error norm with respect to the mesh size $h$ and the timestep $\tau$. The results are obtained by linearly reducing both $h$ and $\tau$ at the same time. The convergence history reported in Figure 3 shows that the error decays with linear rate, as expected from the a priori estimate in Theorem 3.3, and the a posteriori estimator decays with the same (linear) rate.

6.2 Linearization analysis

We now numerically assess the validity of the a posteriori estimate concerning the linearization error. In order to reduce as much as possible the numerical error induced by spatial and temporal approximations, we perform the the numerical experiments with the same discretization parameters ($h_{ref}=4\cdot 10^{-3}$, $\tau_{ref}=2\cdot 10^{-3}$) employed to build the high-fidelity numerical solution. Selecting an instant $t_{n}$, we compute several iterations of the Newton scheme (2.8)-(2.9) until the convergence criterion is satisfied with $tol=10^{-15}$. The iterative scheme produces a sequence $\{u_{h,K_{n}}^{n},w_{h,K_{n}}^{n}\}_{k=0,\ldots,K}$. Then, for each $k$ we compute $\gamma_{k}^{n}$ and compare it with the linearization error. In Figure 4 we report the described comparison at $t_{n}=2.5$ and $t_{n}=10$.

We observe that for each $k=1,\ldots,K$ the estimator is above the error, and they decrease with the same rate.

7 Conclusions

We considered the numerical approximation of the monodomain model, a system of a parabolic semilinear reaction-diffusion equation coupled with a nonlinear ordinary differential equation. The monodomain model arises from the (simplified) mathematical description of the electrical activity of the heart. In particular, we derived a posteriori error estimators accounting for different sources of error (space/time discretization and linearization). Moreover, after obtaining an a priori error estimate, we showed reliability and efficiency (this latter under a suitable assumption) of the error indicators. Lastly, a set of numerical experiments assess the validity of the theoretical results.