A Posteriori Error Estimates for Fully Discrete Finite Element Method for Generalized Diffusion Equation with Delay

Abstract

In this paper, we derive several a posteriori error estimators for generalized diffusion equation with delay in a convex polygonal domain. The Crank–Nicolson method for time discretization is used and a continuous, piecewise linear finite element space is employed for the space discretization. The a posteriori error estimators corresponding to space discretization are derived by using the interpolation estimates. Two different continuous, piecewise quadratic reconstructions are used to obtain the error due to the time discretization. To estimate the error in the approximation of the delay term, linear approximations of the delay term are used in a crucial way. As a consequence, a posteriori upper and lower error bounds for fully discrete approximation are derived for the first time. In particular, long-time a posteriori error estimates are obtained for stable systems. Numerical experiments are presented which confirm our theoretical results.

1 Introduction

In this paper we study a posteriori error estimates of fully discrete finite element method for parabolic delay differential equations (PDDEs). As an example, we consider the numerical methods for the generalized diffusion equation with delay

$$\begin{aligned} \partial _t u(x,t)+{\mathcal A}u(x,t)+{\mathcal B}u(x,\eta (t))= & {} f(x,t),\quad (x,t)\in \Omega \times [0,T], \end{aligned}$$

(1.1)

$$\begin{aligned} u(x,t)= & {} 0,~~~\qquad (x,t)\in \partial \Omega \times [0,T], \end{aligned}$$

(1.2)

$$\begin{aligned} u(x,t)= & {} \phi (x,t),\quad (x,t)\in \Omega \times [t_{*},0], \end{aligned}$$

(1.3)

where $\Omega \subset {\mathbb R}^d$, $d\in \mathbb N$, is a bounded convex polygonal domain with boundary $\partial \Omega$, $\partial _tu(x,t)=\frac{\partial u}{\partial t}(x,t)$, $T\le \infty$ is a constant, $\eta (t)=t-\tau (t)<t$ is a sufficiently smooth “lag”  function on the interval [0, T] and increases strictly monotonically, the initial function $\phi$ and the forcing term f(x, t) are assumed to be smooth, and $t_{*}=\eta (0)$ is a constant. More precisely, we shall derive a posteriori error estimators for the Crank–Nicolson finite element (CNFE) solutions to Eqs. (1.1)–(1.3). Here operators $\mathcal A$ and $\mathcal B$ are of the form, respectively,

$$\begin{aligned} \mathcal Au=-\nabla \cdot (A(x)\nabla u),\qquad \mathcal Bu=-\nabla \cdot (B(x)\nabla u), \end{aligned}$$

with $\mathcal A$ being self-adjoint positive definite, where “$\nabla$”  denotes the spatial gradient and the coefficient matrices A and B are assumed to be smooth.

The equation of the form (1.1), to the best of our knowledge, is one of the two typical examples of PDDEs (see, for example, [21, 39]). Equations of this type arise in several applications such as control theory, population dynamics, heat conduction in materials with thermal memory, biosciences, and so on (see [2, 15, 28, 39, 49], and references therein). Especially, delay systems have recently been used to model the outbreak of Corona Virus Disease 2019 ( COVID-19) in the world [10, 50].

For the numerical solutions of (1.1)–(1.3), the stability of predictor–corrector methods was first investigated in [21]. Later on, several researchers have discussed the stability and a priori error analysis of numerical methods for (1.1)–(1.3) with a constant delay, i.e., $\tau (t)=t-\eta (t)=\tau$ with $\tau >0$ being a real constant. The authors of [13, 14] considered the consistence, stability and convergence of an explicit difference scheme and two implicit difference schemes for the Eqs. (1.1)–(1.3) with a constant delay, respectively. Blanco-Cocom and Ávila-Vales in [7] studied the stability and convergence of the $\theta$-methods together with finite difference methods for constant delay reaction–diffusion equations. Using finite element methods for spatial discretization, Liang in [33] investigated the convergence and asymptotic stability of backward Euler and Crank–Nicolson schemes for (1.1)–(1.3) with a constant delay. It should be pointed out that applying the process of semi-discretization with respect to the spatial variable x to (1.1)–(1.3) will leads to delay differential equations (DDEs), and stability and convergence of numerical methods for DDEs have been investigated by many researchers; see, e.g., [6, 12, 16, 17, 22,23,24,25,26,27, 30, 31, 34, 36, 43,44,45,46, 51, 52]. For recent monographs relevant to the numerical solutions of DDEs, we refer the readers to [5, 9, 29, 32].

We write $\xi _0=0$ and introduce the points $\{\xi _i\}_{i\ge 1}$ such that

$$\begin{aligned} \eta (\xi _{i+1})=\xi _i, \qquad i=0,1,2,\ldots . \end{aligned}$$

Under the given assumptions, the solution of (1.1)–(1.3) on $(\xi _i,\xi _{i+1}]$ follows from the solution on $(\xi _{i-1},\xi _{i}]$ and the solution for $t\ge 0$ can be found by a “method of steps”. Although u(x, t) will inherit its smoothness properties on each interval $[\xi _i,\xi _{i+1}]$ from the smoothness properties of $\phi$ on $[t_{*},0]$ for smooth f and $\eta$, the solution u(x, t) may suffer regularity loss of time smoothness at points $\xi _i$, $i=0,1,2,\ldots$; see, e.g., [5]. This implies that even if the initial function $\phi$ is sufficiently smooth, the solution to DDEs or PDDEs is generally not sufficiently smooth and has breaking points $\xi _i$, $i=0,1,\ldots$. Several approaches have been proposed to detect and approximate these breaking points (see, e.g., [3, 11, 18, 19, 31, 48]) or to resolve the low regularity of the solution at these breaking points by regularization (see, e.g., [20]). These breaking points requires us to provide error estimates with minimal regularity of the solutions. A posteriori estimators of numerical methods can be also used as the indicator in the adaptive time grid refinement which is important in numerically solving DDEs or PDDEs because of the existence of the breaking points $\xi _i$, $i=0,1,\ldots$. Recently, Wang et. al. [47] have derived a posteriori error bounds for the Crank–Nicolson–Galerkin type time semi-discretizations for reaction–diffusion equations with delay. In that paper, numerical studies were reported for several test cases where the Crank–Nicolson–Galerkin type time discretization methods together with finite difference method were explored. From the numerical results presented in [47], we observed that although the errors at the point of low regularity have influence on the error estimators, the Crank–Nicolson method attains its expected second order convergence rate since $\xi _1$ is the only low regularity point for second order time methods. Therefore, in this paper we still consider the Crank–Nicolson method in time discretization but with finite element methods in space discretization for (1.1)–(1.3). To the best of our knowledge no article is available in the literature concerning a posteriori error analysis for fully discrete finite element approximations for PDDEs.

In the absence of the delay term (when $\mathcal B=0$), much work has been conducted on a posteriori error estimates for linear parabolic problems on isotropic meshes (see, e.g., [4, 42]) or anisotropic meshes (see, e.g., [35, 38]) during the last decades. In particular, several researchers have investigated a posteriori error analysis for the Crank–Nicolson method for parabolic problems in recent years; see, e.g., [1, 4, 35, 37, 38, 42]. It is natural to ask whether the a posteriori error analysis for the parabolic problem can be carried over to PDDEs which may be thought of as a perturbation of the parabolic problem. Such an extension is particularly important not only because the solution to PDDEs has some weak discontinuity points $\xi _i$, $i=0,1,\ldots$, but also because a posteriori error analysis has received little attention in the literature. Due to the low regularity property of the solution and the presence of the delay term in (1.1), the extension of those ideas from the simple heat equation to the linear PDDEs is of increased difficulty. To obtain a posteriori error estimators for the CNFE method for linear PDDEs, in this paper, we introduce new quadratic reconstructions based on interpolation approximations.

The paper is organized as follows. We start in Sect. 2 by introducing some definitions and notations relative to the problem, and by presenting the time and space discretization of the problem. In this section, the stability of generalized diffusion equations with delay is also discussed. Section 3 is devoted to quadratic reconstructions for CNFE method for this class of equations. The a posteriori error estimates for the two reconstructions are presented in Sect 4. Especially, the long-time a posteriori error bounds are derived for CNFE method for stable systems in this section. A numerical study is carried out for several test cases in Sect. 5. Section 6 contains a few concluding remarks.

2 CNFE for Generalized Diffusion Equations with Delay

In this section we consider the CNFE approximation of generalized diffusion equations with delay (1.1)–(1.3). The continuous time approximation will be obtained by linear interpolation.

2.1 Functional Setting and Abstract Formulation

We first include some preliminaries on the functional setting and some basic results to be used in our analysis.

For a given Lebesgue measurable set $\omega \subset {\mathbb R}^d$, let $L^p(\omega )$, $p\ge 1$, be the Lebesgue spaces with the corresponding norms $\Vert \cdot \Vert _{L^p(\omega )}$, and $W^{s,p}(\omega )$, $s\ge 0$, $p\ge 1$, be the standard Sobolev spaces with the norms $\Vert \cdot \Vert _{s,p,\omega }$. For convenience, we denote by $\Vert \cdot \Vert _{s,\omega }$ and $\Vert \cdot \Vert _{\infty ,\omega }$ the norms of the space $H^s(\omega )=W^{s,2}(\omega )$ and $L^\infty (\omega )$, respectively, and $\Vert \cdot \Vert _{\omega }$ for the usual norm in $L^2(\omega )$. We denote by $H^1_0(\omega )$ the space of functions in $H^1(\omega )$ that vanish on the boundary of $\omega$ (boundary values are taken in the sense of traces), and we denote by $\langle \cdot ,\cdot \rangle _\omega$ either the inner product in $L^2(\omega )$ or the duality pairing between $H^1_0(\omega )$ and its duality $H^{-1}(\omega )$. The norm of the space $H^{-1}(\omega )$ will be denoted by $\Vert \cdot \Vert _{-1,\omega }$. Whenever $\omega =\Omega$, we omit the letter $\omega$ in the subscripts of $\Vert \cdot \Vert _{s,p,\omega }$, $\Vert \cdot \Vert _{s,\omega }$, $\Vert \cdot \Vert _{\infty ,\omega }$, $\Vert \cdot \Vert _{\omega }$, $\langle \cdot ,\cdot \rangle _\omega$, and $\Vert \cdot \Vert _{-1,\omega }$. We also introduce the Bochner spaces $L^2(a,b;X)$ which are defined by

$$\begin{aligned} L^2(a,b;X):=\left\{ \varphi (t)\in X, t\in [a,b]\left| \int ^b_a\Vert \varphi \Vert ^2_Xdt<\infty \right. \right\} , \end{aligned}$$

with X being $H^{-1}(\Omega )$, $H^1_0(\Omega )$, or $H^1(\Omega )$.

To state the abstract weak formulation of problem (1.1), we need to define two bilinear forms, $a(\cdot ,\cdot )$ and $b(\cdot ,\cdot )$, corresponding to the elliptic operator $\mathcal A$ and $\mathcal B$, respectively. They are defined respectively on $H^1_0(\Omega )\times H^1_0(\Omega )$ by

$$\begin{aligned} a(\varphi ,\chi ):=\langle A(x)\nabla \varphi ,\nabla \chi \rangle ,\qquad b(\varphi ,\chi ):=\langle B(x)\nabla \varphi ,\nabla \chi \rangle ,\qquad \forall \varphi ,\chi \in H^1_0(\Omega ). \end{aligned}$$

With the above notations and assumptions in hand, the weak formulation of problem (1.1)–(1.3) may be stated as follows: Suppose $f\in L^2(0,T;H^{-1}(\Omega ))$ and $\phi \in L^2(t_{*},0;H^1_0(\Omega ))$, and seek $u: [0,T]\rightarrow H^1_0(\Omega )$ such that

$$\begin{aligned} \langle u_t,\varphi \rangle +a(u(t),\varphi )+b(u(\eta (t)),\varphi )= & {} \langle f,\varphi \rangle , ~~~\forall \varphi \in H^1_0(\Omega ),~t\in [0,T], \end{aligned}$$

(2.1)

$$\begin{aligned} u(\cdot ,t)= & {} \phi (\cdot ,t),\quad t\in [t_{*},0]. \end{aligned}$$

(2.2)

Here we assume that the bilinear form $a(\cdot ,\cdot )$ is coercive and continuous on $H^1_0(\Omega )$, i.e.,

$$\begin{aligned} a(\varphi ,\varphi )\ge \alpha \Vert \nabla \varphi \Vert ^2,\qquad |a(\varphi ,\chi )|\le \beta \Vert \nabla \varphi \Vert \Vert \nabla \chi \Vert ,\qquad \forall \varphi ,\chi \in H^1_0(\Omega ) \end{aligned}$$

(2.3)

with $\alpha ,\beta \in \mathbb R^+$. We also assume that the bilinear form $b(\cdot ,\cdot )$ is continuous on $H^1_0(\Omega )$, i.e.,

$$\begin{aligned} |b(\varphi ,\chi )|\le \gamma \Vert \nabla \varphi \Vert \Vert \nabla \chi \Vert ,\qquad \forall \varphi ,\chi \in H^1_0(\Omega ), \end{aligned}$$

(2.4)

with $\gamma \in \mathbb R^+$.

2.2 Stability of Generalized Diffusion Equations with Delay

Since we will derive long-time a posteriori error estimates for CNFE method for stable generalized diffusion equations with delay, in this subsection we discuss shortly the stability of this class of equations. For this purpose, we assume that the function $\eta (t)$ is continuously differentiable and define

$$\begin{aligned} \bar{\eta }=\sup \limits _{s\ge 0}\frac{1}{\eta ^\prime (s)}. \end{aligned}$$

(2.5)

Further, let us define $\eta ^{(0)}(t)=t$, $\eta ^{(1)}(t)=\eta (t)$, and recursively

$$\begin{aligned} \eta ^{(\mu )}(t)=\eta (\eta ^{(\mu -1)}(t)),\qquad \mu \in \mathbb N. \end{aligned}$$

We assume that there exists a positive integer M such that $\eta ^{(M)}(T)\in [t_{*},0]$. Wang et. al. [47] illustrated the existence of such a positive integer M for three lag functions $\eta (t)=t-1, 0.5t-1, t-\sqrt{t+1}$. For the three lag functions, we have $\bar{\eta }=1,~2,~2$, respectively. Here we give another nontrivial example of a lag function $\eta$ satisfying condition $\eta ^{(M)}(T)\in [t_{*},0]$: $\eta (t)=pt-\ln (t+2)$, $0.5<p<1$. Obviously, for any fixed $T>0$, there exists a positive integer M such that $\eta ^{(M)}(T)\in [-\ln 2,0]$. It is also easy to get $\bar{\eta }=\frac{2}{2p-1}$ for this case.

Now we are in position to give the stability result.

Theorem 2.1

(Stability) Suppose $f\in L^2(0,T;H^{-1}(\Omega ))$ and $\phi \in L^2(t_{*},0;H^1_0(\Omega ))$, and the bilinear forms $a(\cdot ,\cdot )$ and $b(\cdot ,\cdot )$ satisfy (2.3) and (2.4), respectively. Then under the condition

$$\begin{aligned} \alpha ^2>\gamma ^2\bar{\eta }, \end{aligned}$$

(2.6)

the solution u(t) to (2.1)–(2.2) satisfies

$$\begin{aligned}&\Vert u(T)\Vert ^2+\alpha \int ^T_{\eta (T)}\Vert \nabla u(t)\Vert ^2dt\nonumber \\&\qquad \le \Vert u(0)\Vert ^2+\frac{\alpha }{\bar{\eta }}\int ^0_{t_{*}}\Vert \nabla \phi (t)\Vert ^2dt+\frac{1}{C_\alpha }\int ^T_0\Vert f\Vert _{-1}^2dt, \end{aligned}$$

(2.7)

where $C_\alpha =\frac{\alpha ^2-\gamma ^2\bar{\eta }}{\alpha }$.

Proof

Setting $\varphi =u(t)$ in (2.1), we obtain

$$\begin{aligned} \langle u_t,u(t)\rangle +a(u(t),u(t))+b(u(\eta (t)),u(t))=\langle f,u(t)\rangle , \end{aligned}$$

(2.8)

which, in view of (2.3) and (2.4), implies that

$$\begin{aligned} \frac{1}{2}\frac{d}{dt}\Vert u(t)\Vert ^2+\alpha \Vert \nabla u(t)\Vert ^2\le & {} \gamma \Vert \nabla u(\eta (t))\Vert \Vert \nabla u(t)\Vert +\Vert f\Vert _{-1}\Vert \nabla u(t)\Vert \nonumber \\\le & {} \frac{\alpha }{2}\Vert \nabla u(t)\Vert ^2+\frac{\alpha }{2\bar{\eta }}\Vert \nabla u(\eta (t))\Vert ^2+\frac{1}{2C_\alpha }\Vert f\Vert _{-1}^2.~~~~ \end{aligned}$$

(2.9)

Integrate from 0 to T to get

$$\begin{aligned}&\Vert u(T)\Vert ^2 + \alpha \int ^T_0\Vert \nabla u(t)\Vert ^2dt\nonumber \\&\qquad \le \Vert u(0)\Vert ^2+\frac{\alpha }{\bar{\eta }}\int ^T_0\Vert \nabla u(\eta (t))\Vert ^2dt+\frac{1}{C_\alpha }\int ^T_0\Vert f\Vert _{-1}^2dt. \end{aligned}$$

(2.10)

Using (2.5), we have

$$\begin{aligned} \frac{\alpha }{\bar{\eta }}\int _{0}^{T}\Vert \nabla u(\eta (t))\Vert ^2dt\le \alpha \int _{\eta (0)}^{\eta (T)}\Vert \nabla u(t)\Vert ^2dt, \end{aligned}$$

(2.11)

and therefore (2.7) holds. This completes the proof. $\square$

We note that for constant delay equation, i.e., $\eta (t)=t-\tau$ (in this case, $\bar{\eta }=1$), Tian in [41] obtained a similar condition to (2.6) for stability.

2.3 CNFE Method for Generalized Diffusion Equations with Delay

For any $0<h<1$, let $\mathcal T_h$ be a shape regular, conforming triangulation of $\overline{\Omega }$ into triangles K with diameter $h_K\le h$. We define $V_h$ by the usual finite element space of continuous, piecewise linear functions on $\mathcal T_h$:

$$\begin{aligned} V_{h}=\left\{ v_h\in H^1(\Omega )|v_{h|K}\in \mathcal P_1(K), ~\forall K\in \mathcal T_h\right\} , \end{aligned}$$

where $\mathcal P_1(K)$ denotes the space of polynomials of degree at most 1. We now set

$$\begin{aligned} V^0_{h}=V_h\cap H^1_0(\Omega ). \end{aligned}$$

The assumptions on the mesh and discrete space above ensure the existence of an interpolation operator $\Pi _h:H^1_0(\Omega )\rightarrow V^0_h$, satisfying the following approximation estimates.

Proposition 2.2

[8, 40] Let $\Pi _h:H^1_0(\Omega )\rightarrow V^0_h$ be the usual linear Lagrange interpolation operator. Then we have

$$\begin{aligned} \Vert \nabla \Pi _hv\Vert \le C_1\Vert \nabla v\Vert . \end{aligned}$$

(2.12)

Further, the following approximation property is satisfied:

$$\begin{aligned} h^{-1}_K\Vert v-\Pi _hv\Vert _K+h^{-1/2}_K\Vert v-\Pi _hv\Vert _{\partial K}\le C_{\Pi }\Vert \nabla v\Vert _K, \end{aligned}$$

(2.13)

where the constants $C_1$ and $C_{\Pi }$ depend only on the shape regularity of the family of triangulations $\mathcal T_h$.

In order to describe the time discretization corresponding to (2.1)–(2.2), we introduce a time mesh $0=t_0<t_1<\cdots <t_N=T$ which is a partition of the interval [0, T] with $I_n:=[t_{n-1},t_n]$ and we denote the time steps by $k_n=t_n-t_{n-1}$. For $n=1,2,\ldots ,N$, we define

$$\begin{aligned} t_{n-1/2}:=\frac{t_n+t_{n-1}}{2},~~v^{n-1/2}:=\frac{v^n+v^{n-1}}{2},~~{\hbox {and}}~~\bar{\partial }v^n:=\frac{v^n-v^{n-1}}{k_n}. \end{aligned}$$

Let $U_h(t)=\Pi _h\phi (t)$ for $t\in [t_{*},0]$. Then our CNFE approximation can be stated as follows: For the given $U_h(t)$, $t\in [t_{*},0]$, find $U^n_h\in V^0_h$ for all $n=1,2,\ldots ,N$ such that

$$\begin{aligned} \langle \bar{\partial }U^n_h,\varphi _h\rangle +a(U^{n-1/2}_h,\varphi _h)+b\left( U_h(\eta (t_{n-1/2})),\varphi _h\right) =\langle f(t_{n-1/2}),\varphi _h\rangle ,~ \forall \varphi _h\in V^0_h,\nonumber \\ \end{aligned}$$

(2.14)

where the continuous time approximation $U_h(t)~(t\in [0,T])$ is defined by linearly interpolating the nodal values $U^{n}_h$ and $U^{n-1}_h$ for $t\in I_n$, $1\le n\le N$, i.e.

$$\begin{aligned} U_h(t):= \frac{t-t_{n-1}}{k_n}U^{n}_h+\frac{t_{n}-t}{k_n}U^{n-1}_h=U^{n-1/2}_h+(t-t_{n-1/2})\bar{\partial }U^{n}_h. \end{aligned}$$

(2.15)

It is noteworthy that the method (2.14) itself needs the continuous approximation $U_h(t)$.

3 Quadratic Reconstructions for Generalized Diffusion Equations with Delay

Now we shall introduce quadratic time reconstructions, which are appropriate for Eq. (2.1) corresponding method (2.14), and continuous time approximation (2.15).

3.1 Crank–Nicolson Delay-Dependent Reconstruction

For $1\le n\le N$, we define the continuous, piecewise quadratic delay-dependent reconstruction ${\widetilde{U}}_h: I_n\rightarrow V^0_h$ of $U_h$ in time as

$$\begin{aligned} {\widetilde{U}}_h(t):= & {} U_h(t)+\frac{1}{2}(t-t_{n-1})(t-t_{n}){\widetilde{W}}^{n}_h\nonumber \\= & {} U^{n-1/2}_h+(t-t_{n-1/2})\bar{\partial }U^{n}_h+\frac{1}{2}(t-t_{n-1})(t-t_{n}){\widetilde{W}}^{n}_h, \end{aligned}$$

(3.1)

where ${\widetilde{W}}^{n}_h\in V^0_h$ is defined by

$$\begin{aligned} \langle {\widetilde{W}}^{n}_h,\varphi _h\rangle =\langle \bar{\partial }f^{n},\varphi _h\rangle -a(\bar{\partial }U^{n}_h,\varphi _h)-\langle \eta ^\prime (t_{n-1/2})B\nabla \bar{\partial }_\eta U^{n}_h,\nabla \varphi _h\rangle . \end{aligned}$$

(3.2)

Here we define $\bar{\partial }_\eta U^{n}_h$ as

$$\begin{aligned} \bar{\partial }_\eta U^{n}_h:=\frac{U_h(\eta (t_{n}))-U_h(\eta (t_{n-1}))}{\eta (t_{n})-\eta (t_{n-1})}. \end{aligned}$$

(3.3)

For $t\in [t_{*},0]$, we let ${\widetilde{U}}_h(t)=U_h(t)$. Note that even if $\eta (t_{n-1/2})\in [t_{i-1},t_i]$, $i\le n-1$, $\eta (t)$ for $t\in [t_{n-1},t_n]$ may not belong to $[t_{i-1},t_i]$. Then in order to estimate the interpolation error in approximating the delay term, we set

$$\begin{aligned} \widetilde{R}_U(t):=U_h(\eta (t))-U_h(\eta (t_{n-1/2}))-(t-t_{n-1/2})\eta ^\prime (t_{n-1/2})\bar{\partial }_\eta U^{n}_h,~ t\in [t_{n-1},t_n].~~~~ \end{aligned}$$

(3.4)

Obviously, if $\eta (t_{n})$ and $\eta (t_{n-1})$ belong to the same small interval $[t_{i-1},t_i]$, $i\le n-1$, then we have $\eta (t)\in [t_{i-1},t_i]$ for $t\in [t_{n-1},t_n]$ and hence $\widetilde{R}_U(t)=0$ because $\eta$ is a strictly monotone increasing function on [0, T].

3.2 Multi-point Delay-Independent Reconstruction

We observe that the reconstruction introduced in the above subsection depends on the past values $U_h(\eta (t))$, and ${\widetilde{W}}^{n}_h$ in (3.1) is formally an approximation of $\partial f/\partial t-\mathcal A \partial u/\partial t-\mathcal B\partial u(\eta (t))/\partial t$ in the time slab $I_n$. Since the latter is equal to $\partial ^2 u/\partial t^2$, it seems natural to try to replace ${\widetilde{W}}^{n}_h$ in (3.1) by a finite difference approximation of $\partial ^2 u/\partial t^2$. We introduce thus a continuous, piecewise multi-point delay-independent reconstruction ${\widehat{U}}_h$ of $U_h$ in time defined for all $t\in I_n$, $1\le n\le N$, by

$$\begin{aligned} {\widehat{U}}_h(t):= & {} U_h(t)+\frac{1}{2}(t-t_{n-1})(t-t_{n}){\widehat{W}}^{n}_h\nonumber \\= & {} U^{n-1/2}_h+(t-t_{n-1/2})\bar{\partial }U^{n}_h+\frac{1}{2}(t-t_{n-1})(t-t_{n}){\widehat{W}}^{n}_h, \end{aligned}$$

(3.5)

where

$$\begin{aligned} {\widehat{W}}^{n}_h:=\left\{ \begin{array}{ll} {\widetilde{W}}^1_h &{}\quad {\hbox {for}}~~~ n=1,\\ \bar{\partial }^2 U^{n}_h:=\frac{2(\bar{\partial }U^{n}_h-\bar{\partial }U^{n-1}_h)}{k_{n}+k_{n-1}},&{}\quad {\hbox {for}}~~~ 2\le n\le N.\end{array}\right. \end{aligned}$$

An argument similar to [35] ensures that $\widehat{U}_h(t)$ vanishes on the boundary so that ${\widehat{U}}_h(t),~ \widehat{W}^{n}_h\in V^0_h$. We note that $U_h(t),~{\widetilde{U}}_h(t)$ and ${\widehat{U}}_h(t)$ coincide at $t_i$, $i=1,2,\ldots ,N$. For $t\in [t_{n-3/2},t_{n-1/2}]$, $\eta (t)$, $\eta (t_{n-3/2})$ and $\eta (t_{n-1/2})$ may not belong to the same interval $I_i$, $i\le n-1$. Then, to estimate the time discretization error due to quadrature approximation of the delay term in the context of multi-point reconstruction, we define

$$\begin{aligned} {\widehat{R}}_U(t):=U_h(\eta (t))-l_{n-1/2}(t)U_h(\eta (t_{n-1/2}))-l_{n-3/2}(t)U_h(\eta (t_{n-3/2})), \end{aligned}$$

(3.6)

where

$$\begin{aligned} l_{n-1/2}(t)=1+\frac{2(t-t_{n-1/2})}{k_n+k_{n-1}},\qquad l_{n-3/2}(t)=-\frac{2(t-t_{n-1/2})}{k_n+k_{n-1}}. \end{aligned}$$

Remark 3.1

Considering the delay property of the problem, it is natural to use the multi-point reconstruction. For $t_{*}<0$, we can define $k_0=-t_{*}$, $U^{-1}_h=I_h\phi (t_{*})$ and

$$\begin{aligned} {\widehat{W}}^{1}_h:= \bar{\partial }^2 U^{1}_h:=\frac{2((U^1_h-U^{0}_h)/k_1-(U^{0}_h-U^{-1}_h)/k_0)}{k_{1}+k_{0}}. \end{aligned}$$

(3.7)

Nevertheless, we only consider the case of $\widehat{W}^{1}_h={\widetilde{W}}^1_h$ for simplicity in this paper.

4 A Posteriori Error Estimates

Based on the reconstructions ${\widetilde{U}}_h$ and ${\widehat{U}}_h$, in this section, we derive four optimal order a posteriori upper bounds for the error $e:=u-U_h$ in the $L^2(H^1)$-norm, where u and $U_h$ satisfy (2.1) and (2.14), respectively.

4.1 Long-Time Error Estimates for Stable Systems

We first derive long-time a posteriori error estimates in the $L^2(\eta (T),T;H^1)$ norm. For this purpose, we set $\widetilde{e}:=u-{\widetilde{U}}_h$ and ${\widehat{e}}:=u-{\widehat{U}}_h$ for the reconstructions ${\widetilde{U}}_h$ and ${\widehat{U}}_h$, and need to prove the following lemma.

Lemma 4.1

Let $U_h$ be defined by (2.14). Then

• 1. For the reconstruction ${\widetilde{U}}_h$ defined by (3.1), we have, for any $t\in I_n$ with $1\le n\le N$, and for all $\varphi _h\in V^0_h$,

  $$\begin{aligned}&\left\langle \partial _t {\widetilde{U}}_h,\varphi _h\right\rangle +a(U_h,\varphi _h)+b\left( U_h(\eta (t_{n-1/2})),\varphi _h\right) \nonumber \\&\qquad =\langle {\tilde{f}},\varphi _h\rangle -(t-t_{n-1/2})b\left( \eta ^\prime (t_{n-1/2})\bar{\partial }_\eta U^n_h,\varphi _h\right) , \end{aligned}$$

  (4.1)

  where $\bar{\partial }_\eta U^{n}_h$ has been given in (3.3) and ${\tilde{f}}(t)=f(t_{n-1/2})+(t-t_{n-1/2})\partial f^{n}$, $t\in I_n$.

• 2. For the multi-point reconstruction ${\widehat{U}}_h$ defined by (3.5), we have, for any $t\in I_n$ with $2\le n\le N$, and for all $\varphi _h\in V^0_h$,

  $$\begin{aligned}&\left\langle \partial _t {\widehat{U}}_h,\varphi _h\right\rangle +a(U_h,\varphi _h)+l_{n-1/2}(t)b(U_h(\eta (t_{n-1/2})),\varphi _h)\nonumber \\&\qquad +\,l_{n-3/2}(t)b(U_h(\eta (t_{n-3/2})),\varphi _h), \nonumber \\&\qquad =\,\langle \hat{f},\varphi _h\rangle +\frac{k_{n-1}(t-t_{n-1/2})}{2}a(\widehat{W}^{n}_h,\varphi _h), \end{aligned}$$

  (4.2)

  where

  $$\begin{aligned} {\hat{f}}=f(t_{n-1/2})+(t-t_{n-1/2})\frac{2(f(t_{n-1/2})-f(t_{n-3/2}))}{k_{n-1}+k_n}.\end{aligned}$$

Proof

For $\varphi _h\in V^0_h$ and $1\le n\le N$, (2.14) and (2.15) can yields

$$\begin{aligned}&\langle \bar{\partial }U^{n}_h,\varphi _h\rangle +a(U_h,\varphi _h)+b\left( U_h(\eta (t_{n-1/2})),\varphi _h\right) \nonumber \\&\qquad =\langle f(t_{n-1/2}),\varphi _h\rangle +(t-t_{n-1/2})a(\bar{\partial }U^{n}_h,\varphi _h). \end{aligned}$$

(4.3)

It follows from (3.1) that

$$\begin{aligned} \partial _t {\widetilde{U}}_h(t)=\bar{\partial }U_h^{n}+\left( t-t_{n-1/2}\right) {\widetilde{W}}^{n}_h. \end{aligned}$$

(4.4)

Substitute the above relation into (4.3) to obtain

$$\begin{aligned}&\left\langle \partial _t \widetilde{U}_h,\varphi _h\right\rangle +a(U_h,\varphi _h)+b\left( U_h(\eta (t_{n-1/2})),\varphi _h\right) \nonumber \\&\qquad =\langle f(t_{n-1/2}),\varphi _h\rangle +(t-t_{n-1/2})\left( a(\bar{\partial }U^{n}_h,\varphi _h)+\langle \widetilde{W}^{n}_h,\varphi _h\rangle \right) . \end{aligned}$$

(4.5)

Then, by the definition of ${\widetilde{W}}^{n}_h$, we give the desired identity (4.1).

To derive (4.2) involving the reconstruction ${\widehat{U}}_h$, from (3.5) we get

$$\begin{aligned} \partial _t {\widehat{U}}_h(t)=\bar{\partial }U_h^{n}+\left( t-t_{n-1/2}\right) {\widehat{W}}^{n}_h, \quad t\in I_n,~~ n=2,3,\ldots ,N. \end{aligned}$$

(4.6)

For $\varphi _h\in V^0_h$ and $2\le n\le N$, substitute the above equation into (4.3) to obtain

$$\begin{aligned}&\left\langle \partial _t \widehat{U}_h(t),\varphi _h\right\rangle +a(U_h,\varphi _h)+b\left( U_h(\eta (t_{n-1/2})),\varphi _h\right) \nonumber \\&\qquad =\langle f(t_{n-1/2}),\varphi _h\rangle +(t-t_{n-1/2})\left( a(\bar{\partial }U^{n}_h,\varphi _h)+\langle {\widehat{W}}^{n}_h,\varphi _h\rangle \right) . \end{aligned}$$

(4.7)

Replace n by $n-1$ such that (2.14) becomes

$$\begin{aligned} \langle \bar{\partial }U^{n-1}_h,\varphi _h\rangle +a(U^{n-3/2}_h,\varphi _h)+b\left( U_h(\eta (t_{n-3/2})),\varphi _h\right)= & {} \langle f(t_{n-3/2}),\varphi _h\rangle ,\nonumber \\&\quad \forall \varphi _h\in V^0_h. \end{aligned}$$

(4.8)

Subtract (4.8) from (2.14) to obtain for all $\varphi _h\in V^0_h$

$$\begin{aligned}&\left\langle {\widehat{W}}^{n}_h,\varphi _h\right\rangle +a\left( \frac{U^n_h-U^{n-2}_h}{k_{n-1}+k_n},\varphi _h\right) \nonumber \\&\qquad +\,\frac{2}{k_{n-1}+k_n}b(U_h(\eta (t_{n-1/2}))-U_h(\eta (t_{n-3/2})),\varphi _h)\nonumber \\&\qquad =\,\left\langle \frac{2(f(t_{n-1/2})-f(t_{n-3/2}))}{k_{n-1}+k_n},\varphi _h\right\rangle . \end{aligned}$$

(4.9)

Using the relation

$$\begin{aligned} \frac{U^n_h-U^{n-2}_h}{k_{n-1}+k_n}=\bar{\partial }U^n_h-\frac{k_{n-1}}{2}{\widehat{W}}^{n}_h, \end{aligned}$$

from (4.9) we obtain

$$\begin{aligned}&\left\langle \widehat{W}^{n}_h,\varphi _h\right\rangle +a\left( \bar{\partial }U^n_h,\varphi _h\right) +\frac{2}{k_{n-1}+k_n}b(U_h(\eta (t_{n-1/2}))-U_h(\eta (t_{n-3/2})),\varphi _h)\nonumber \\&\qquad =\,\left\langle \frac{2(f(t_{n-1/2})-f(t_{n-3/2}))}{k_{n-1}+k_n},\varphi _h\right\rangle +\frac{k_{n-1}}{2}a(\widehat{W}^{n}_h,\varphi _h). \end{aligned}$$

(4.10)

It suffices now to substitute (4.10) on the right-hand side of (4.7) to obtain (4.2). $\square$

Now we state the following long-time a posteriori error estimates.

Theorem 4.2

(Long-time a posteriori error estimates) Let u and $U_h$ satisfy (2.1) and (2.14), respectively. Then under the stability condition (2.6), there is a constant C independent of the final time T such that, when the delay-dependent reconstruction ${\widetilde{U}}_h$ is employed,

$$\begin{aligned}&\Vert e(T)\Vert ^2+\frac{C_\alpha }{4}\int _{\eta (T)}^{T}\Vert \nabla {\widetilde{U}}_h(t)-\nabla U_h(t)\Vert ^2dt\nonumber \\&\quad \le \Vert e(T)\Vert ^2+\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\nonumber \\&\quad \le \Vert e(0)\Vert ^2+\alpha \int _{\eta (0)}^{0}\Vert \nabla e(t)\Vert ^2dt\nonumber \\&\qquad +\,C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt\right. +\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\nonumber \\&\qquad +\,\int _{t_{n-1}}^{t_{n}}\Vert \nabla (\widetilde{U}_h(t)-U_h(t))\Vert ^2_{K}dt\left. +\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widetilde{R}}_U(t)\Vert ^2_{K}dt\right\} , \end{aligned}$$

(4.11)

and, when the multi-point reconstruction ${\widehat{U}}_h$ is employed,

$$\begin{aligned}&\Vert e(T)\Vert ^2 + \frac{C_\alpha }{4}\int _{\eta (T)}^{T}\Vert \nabla {\widehat{U}}_h(t)-\nabla U_h(t)\Vert ^2dt\nonumber \\&\quad \le \Vert e(T)\Vert ^2+\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\nonumber \\&\quad \le \Vert e(t_1)\Vert ^2+\alpha \int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt+C\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt\right. \nonumber \\&\qquad +\,\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{1}{48} k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\hat{f}}\Vert ^2_{K}dt\nonumber \\&\qquad +\,\int _{t_{n-1}}^{t_{n}}\Vert \nabla (\widehat{U}_h(t)-U_h(t))\Vert ^2_{K}dt\left. +\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt\right\} , \end{aligned}$$

(4.12)

where the contributions $S_{K,n}$ in space discretization are defined on each triangle K of $\mathcal T_h$ and each time interval $I_n$ as

$$\begin{aligned} S_{K,n}:= & {} h_K\left\| f-\bar{\partial }U^{n}_h+\nabla \cdot (A\nabla U_h)+\nabla \cdot (B\nabla U_h(\eta (t)))\right\| _{K}\nonumber \\&+\,h^{1/2}_{K}\Vert [A\nabla U_h\cdot {\varvec{n}}]\Vert _{\partial K}+h^{1/2}_K\Vert [B\nabla U_h(\eta (t))\cdot {\varvec{n}}]\Vert _{\partial K}. \end{aligned}$$

(4.13)

Here and in what follows, $[\cdot ]$ denotes the jump of the bracketed quantity across an interval edge, $[\cdot ]=0$ for an edge on the boundary $\partial \Omega$ and $\varvec{n}$ is the unit edge normal.

Proof

We first show the error estimate (4.11) for the delay-dependent reconstruction ${\widetilde{U}}_h$. For $t\in I_n$, $1\le n\le N$, we use the weak formulation (2.1) and (4.4) to get

$$\begin{aligned}&\left\langle \partial _t {\tilde{e}}(t),\varphi \right\rangle +a(e(t),\varphi )+b(e(\eta (t)),\varphi )\nonumber \\&\qquad =\langle f-\bar{\partial }U^{n}_h,\varphi \rangle -a(U_h,\varphi )-b(U_h(\eta (t)),\varphi )-\left( t-t_{n-1/2}\right) \langle \widetilde{W}^{n}_h,\varphi \rangle .~~~ \end{aligned}$$

(4.14)

Choosing $\varphi ={\tilde{e}}$ in the above error equation and using (4.4) again, we have

$$\begin{aligned}&\left\langle \partial _t {\tilde{e}}(t),{\tilde{e}} \right\rangle +a(e(t),{\tilde{e}})+b(e(\eta (t)),{\tilde{e}})\nonumber \\&\quad =\langle f-\bar{\partial }U^{n}_h,{\tilde{e}}-\Pi _h{\tilde{e}}\rangle -a(U_h,{\tilde{e}}-\Pi _h{\tilde{e}})-b(U_h(\eta (t)),{\tilde{e}}-\Pi _h{\tilde{e}})\nonumber \\&\qquad -\,\left( t-t_{n-1/2}\right) \langle {\widetilde{W}}^{n}_h,{\tilde{e}}-\Pi _h{\tilde{e}}\rangle +\langle f,\Pi _h{\tilde{e}}\rangle -\left\langle \partial _t {\widetilde{U}}_h,\Pi _h{\tilde{e}} \right\rangle \nonumber \\&\qquad -\,a(U_h,\Pi _h{\tilde{e}})-b(U_h(\eta (t),\Pi _h{\tilde{e}}). \end{aligned}$$

(4.15)

Using (4.1) for $\langle \partial _t \widetilde{U}_h,\Pi _h{\tilde{e}} \rangle +a(U_h,\Pi _h{\tilde{e}})$ yields

$$\begin{aligned}&\left\langle \partial _t {\tilde{e}}(t),{\tilde{e}} \right\rangle +a(e(t),{\tilde{e}})+b(e(\eta (t)),{\tilde{e}})\nonumber \\&\quad =\langle f-\bar{\partial }U^{n}_h, {\tilde{e}}-\Pi _h{\tilde{e}}\rangle -a(U_h,{\tilde{e}}-\Pi _h{\tilde{e}})-b(U_h(\eta (t)),{\tilde{e}}-\Pi _h{\tilde{e}})\nonumber \\&\qquad -\,\left( t-t_{n-1/2}\right) \langle {\widetilde{W}}^{n}_h,{\tilde{e}}-\Pi _h{\tilde{e}}\rangle +\langle f-{\tilde{f}}, \Pi _h{\tilde{e}} \rangle -b(U_h(\eta (t)),\Pi _h{\tilde{e}})\nonumber \\&\qquad +\,(t-t_{n-1/2})b\left( \eta ^\prime (t_{n-1/2})\bar{\partial }_\eta U^n_h,\Pi _h{\tilde{e}}\right) +b(U_h(\eta (t_{n-1/2})),\Pi _h{\tilde{e}}). \end{aligned}$$

(4.16)

By the definition of $\tilde{R}_U$ [see (3.4)], we obtain

$$\begin{aligned}&\left\langle \partial _t {\tilde{e}}(t),{\tilde{e}} \right\rangle +a(e(t),{\tilde{e}})+b(e(\eta (t)),{\tilde{e}})\nonumber \\&\quad =\langle f-\bar{\partial }U^{n}_h,{\tilde{e}}-\Pi _h{\tilde{e}}\rangle -a(U_h,{\tilde{e}}-\Pi _h{\tilde{e}})-b(U_h(\eta (t)),{\tilde{e}}-\Pi _h{\tilde{e}})\nonumber \\&\qquad -\,\left( t-t_{n-1/2}\right) \langle {\widetilde{W}}^{n}_h,\tilde{e}-\Pi _h{\tilde{e}}\rangle +\langle f-{\tilde{f}},\Pi _h{\tilde{e}} \rangle -b({\widetilde{R}}_U(t),\Pi _h{\tilde{e}}). \end{aligned}$$

(4.17)

Now, we integrate the second and third terms on the right-hand side of (4.17) by parts on each of the elements K of $\mathcal T_h$. Then, in view of the continuity of the bilinear forms $b(\cdot ,\cdot )$, using the Cauchy-Schwarz inequality and the Poincaré inequality, we obtain

$$\begin{aligned}&\frac{1}{2}\frac{d}{dt}\Vert {\tilde{e}}(t)\Vert ^2 + a(e(t),{\tilde{e}})\le \gamma \Vert \nabla e(\eta (t))\Vert \Vert \nabla {\tilde{e}}\Vert \nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ \left\| f-\bar{\partial }U^{n}_h+\nabla \cdot (A\nabla U_h)+\nabla \cdot (B\nabla U_h(\eta (t)))\right\| _{K}\Vert {\tilde{e}}-\Pi _h{\tilde{e}}\Vert _{K}\right. \nonumber \\&\qquad +\,\frac{1}{2}\Vert [A\nabla U_h\cdot {\varvec{n}}]\Vert _{\partial K}\Vert {\tilde{e}}-\Pi _h{\tilde{e}}\Vert _{\partial K}+\frac{1}{2}\Vert [B\nabla U_h(\eta (t))\cdot {\varvec{n}}]\Vert _{\partial K}\Vert {\tilde{e}}-\Pi _h{\tilde{e}}\Vert _{\partial K}\nonumber \\&\qquad +\,\left| t-t_{n-1/2}\right| \Vert {\widetilde{W}}^{n}_h\Vert _{K}\Vert {\tilde{e}}-\Pi _h{\tilde{e}}\Vert _{K}+C_Ph_K\Vert f-{\tilde{f}}\Vert _{K}\Vert \nabla \Pi _h{\tilde{e}}\Vert _{K}\nonumber \\&\qquad +\,\left. \gamma \Vert \nabla {\widetilde{R}}_U(t)\Vert _{K}\Vert \nabla \Pi _h\tilde{e}\Vert _{K}\right\} , \end{aligned}$$

(4.18)

where $C_P$ is the Poincaré constant. Now, based on the fact

$$\begin{aligned} a(e,{\tilde{e}})=\frac{1}{2}\{a(e,e)+a({\tilde{e}},{\tilde{e}})-a(e-{\tilde{e}}, e-{\tilde{e}})\} \end{aligned}$$

(4.19)

together with the coercivity of the bilinear form $a(\cdot ,\cdot )$ and the relation

$$\begin{aligned} |a(e-{\tilde{e}}, e-{\tilde{e}})|\le \beta \Vert \nabla (e-\tilde{e})\Vert ^2=\beta \Vert \nabla ({\widetilde{U}}_h-U_h)\Vert ^2, \end{aligned}$$

(4.20)

we use the Young’s inequality and Proposition 2.2 to obtain

$$\begin{aligned}&\frac{1}{2}\frac{d}{dt}\Vert {\tilde{e}}(t)\Vert ^2+\frac{\alpha }{2}\Vert \nabla e\Vert ^2+\frac{\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\nonumber \\&\quad \le \frac{\gamma ^2\bar{\eta }}{2\alpha }\Vert \nabla \tilde{e}\Vert ^2+\frac{\alpha }{2\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2 +\sum \limits _{K\in \mathcal T_h}\left\{ C_\Pi S_{K,n}\Vert \nabla {\tilde{e}}\Vert _{K}\right. \nonumber \\&\qquad +\,C_\Pi h_K\left| t-t_{n-1/2}\right| \Vert {\widetilde{W}}^{n}_h\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K}+C_PC_1h_K\Vert f-{\tilde{f}}\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K}\nonumber \\&\qquad +\,\beta \Vert \nabla ({\widetilde{U}}_h(t)-U_h(t))\Vert ^2_{K}+\left. C_1\gamma \Vert {\widetilde{R}}_U(t)\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K} \right\} . \end{aligned}$$

(4.21)

An application of the Young’s inequality yields

$$\begin{aligned}&\frac{d}{dt}\Vert {\tilde{e}}(t)\Vert ^2+\alpha \Vert \nabla e\Vert ^2+\alpha \Vert \nabla {\tilde{e}}\Vert ^2\le \frac{\gamma ^2\bar{\eta }}{\alpha }\Vert \nabla {\tilde{e}}\Vert ^2+\frac{\alpha }{\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ \epsilon C^2_{\Pi } S^2_{K,n}\right. +\epsilon C^2_{\Pi }h^2_{K}\left( t-t_{n-1/2}\right) ^2 \Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\frac{4}{\epsilon }\Vert \nabla {\tilde{e}}\Vert ^2_{K}\nonumber \\&\qquad +\,\left. \epsilon C_P^2C_1^2h_K^2\Vert f-\tilde{f}\Vert ^2_{K}+2\beta \Vert \nabla ({\widetilde{U}}_h(t)-U_h(t))\Vert ^2_{K}+\epsilon \gamma ^2C_1^2\Vert \nabla {\widetilde{R}}_U(t)\Vert ^2_{K}\right\} .\nonumber \\ \end{aligned}$$

(4.22)

Since $\alpha ^2>\gamma ^2\bar{\eta }$, we can choose $\epsilon =\frac{8\alpha }{\alpha ^2-\gamma ^2\bar{\eta }}$. Integrate from $t_{n-1}$ to $t_{n}$ to obtain

$$\begin{aligned}&\Vert {\tilde{e}}(t_{n})\Vert ^2 +\int _{t_{n-1}}^{t_{n}}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\le \Vert {\tilde{e}}(t_{n-1})\Vert ^2+\int _{t_{n-1}}^{t_{n}}\frac{\alpha }{\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2dt\nonumber \\&\qquad +\,C\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla (\widetilde{U}_h(t)-U_h(t))\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} , \end{aligned}$$

(4.23)

where $C=\max \left\{ \frac{8\alpha C^2_\Pi }{\alpha ^2-\gamma ^2\bar{\eta }}, \frac{8\alpha C^2_P C^2_1}{\alpha ^2-\gamma ^2\bar{\eta }}, \frac{8\alpha \gamma ^2 C^2_1}{\alpha ^2-\gamma ^2\bar{\eta }},2\beta \right\}$. Summing over $n=1$ to N, noting ${\tilde{e}}(t_n)=e(t_n)$ for $n\ge 0$, and

$$\begin{aligned} \int _{0}^{T}\frac{\alpha }{\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2dt\le \alpha \int _{\eta (0)}^{\eta (T)}\Vert \nabla e(t)\Vert ^2dt, \end{aligned}$$

(4.24)

lead to the upper bound estimate in (4.11).

We now turn to establish the lower bound of errors. In view of

$$\begin{aligned} \Vert \nabla \tilde{U}(t)-\nabla U(t)\Vert \le \Vert \nabla e(t)\Vert +\Vert \nabla \tilde{e}(t)\Vert \end{aligned}$$

(4.25)

and $\frac{C_\alpha }{2}\le \alpha$, we have

$$\begin{aligned} \frac{C_{\alpha }}{2}\Vert \nabla \tilde{U}(t)-\nabla U(t)\Vert ^2\le 2\left( \alpha \Vert \nabla e(t)\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla \tilde{e}(t)\Vert ^2\right) , \end{aligned}$$

and therefore

$$\begin{aligned}&\Vert e(T)\Vert ^2+\frac{C_\alpha }{4} \int _{\eta (T)}^{T}\Vert \nabla \tilde{U}_h(t)-\nabla U_h(t)\Vert ^2dt \nonumber \\&\qquad \le \Vert e(T)\Vert ^2+\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e(t)\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla \tilde{e}(t)\Vert ^2\right) dt, \end{aligned}$$

(4.26)

which gives the lower bound and completes the proof of (4.11).

Similar argument can be used to estimate the errors for the multi-point reconstruction ${\widehat{U}}_h$. For $t\in I_n$, $2\le n\le N$, using (4.2) for $\langle \partial _t \widehat{U}_h,\Pi _h{\hat{e}} \rangle +a(U_h,\Pi _h{\hat{e}})$ and (3.6), we have

$$\begin{aligned}&\left\langle \partial _t {\hat{e}}(t),{\hat{e}} \right\rangle +a(e(t),{\hat{e}})+b(e(\eta (t)),{\hat{e}})\nonumber \\&\quad =\langle f-\bar{\partial }U^{n}_h,{\hat{e}}-\Pi _h{\hat{e}}\rangle -a(U_h,{\hat{e}}-\Pi _h{\hat{e}})-b(U_h(\eta (t)),{\hat{e}}-\Pi _h{\hat{e}})\nonumber \\&\qquad -\,\left( t-t_{n-1/2}\right) \langle {\widehat{W}}^{n}_h,{\hat{e}}-\Pi _h{\hat{e}}\rangle +\langle f-{\hat{f}},\Pi _h{\hat{e}}\rangle \nonumber \\&\qquad -\,\frac{k_{n-1}(t-t_{n-1/2})}{2}a({\widehat{W}}^{n}_h,\Pi _h\hat{e})-b({\widehat{R}}_U(t),\Pi _h{\hat{e}}). \end{aligned}$$

(4.27)

Similar to the case of the delay-dependent reconstruction, it holds that

$$\begin{aligned}&\frac{1}{2}\frac{d}{dt}\Vert {\hat{e}}(t)\Vert ^2+\frac{\alpha }{2}\Vert \nabla e\Vert ^2+\frac{\alpha }{2}\Vert \nabla {\hat{e}}\Vert ^2\le \frac{\gamma ^2\bar{\eta }}{2\alpha }\Vert \nabla {\hat{e}}\Vert ^2+\frac{\alpha }{2\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ C_\Pi S_{K,n}\Vert \nabla {\hat{e}}\Vert _K+C_\Pi h_K\left| t-t_{n-1/2}\right| \Vert {\widehat{W}}^{n}_h\Vert _{K}\Vert \nabla {\hat{e}}\Vert _K\right. \nonumber \\&\qquad +\,\frac{C_1\beta k_{n-1}|t-t_{n-1/2}|}{2}\Vert \nabla {\widehat{W}}^{n}_h\Vert _K\Vert \nabla {\hat{e}}\Vert _K+C_PC_1h_K\Vert f-{\hat{f}}\Vert _{K}\Vert \nabla {\hat{e}}\Vert _{K}\nonumber \\&\qquad +\,\left. \beta \Vert \nabla ({\widetilde{U}}_h(t)-U_h(t))\Vert ^2_{K}+\gamma C_1\Vert \nabla {\widehat{R}}_{ U}(t)\Vert _{K}\Vert \nabla {\hat{e}}\Vert _{K}\right\} , \end{aligned}$$

(4.28)

and therefore

$$\begin{aligned}&\frac{d}{dt}\Vert {\hat{e}}(t)\Vert ^2+\alpha \Vert \nabla e\Vert ^2+\alpha \Vert \nabla {\hat{e}}\Vert ^2\le \frac{\gamma ^2\bar{\eta }}{\alpha }\Vert \nabla {\hat{e}}\Vert ^2+\frac{\alpha }{\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ \epsilon C^2_\Pi S^2_{K,n}+\,\epsilon C_\Pi ^2h^2_{K}\left( t-t_{n-1/2}\right) ^2 \Vert {\widehat{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\frac{5}{\epsilon }\Vert \nabla {\hat{e}}\Vert ^2_{K}+\frac{\epsilon C^2_1\beta ^2}{4} k^2_{n-1}(t-t_{n-1/2})^2\Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_K+\epsilon C_P^2C_1^2h_K^2\Vert f-{\hat{f}}\Vert ^2_{K}\nonumber \\&\qquad +\,\left. 2\beta \Vert \nabla ({\widehat{U}}_h(t)-U_h(t))\Vert ^2_{K}+\epsilon \gamma ^2C^2_1\Vert \nabla {\widehat{R}}_{ U}(t)\Vert ^2_{K}\right\} . \end{aligned}$$

(4.29)

We choose $\epsilon =\frac{10\alpha }{\alpha ^2-\gamma ^2\bar{\eta }}$ and integrate from $t_{n-1}$ to $t_{n}$ to obtain

$$\begin{aligned}&\Vert {\hat{e}}(t_{n})\Vert ^2+\int _{t_{n-1}}^{t_{n}}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\hat{e}}\Vert ^2\right) dt\le \Vert {\hat{e}}(t_{n-1})\Vert ^2+\int _{t_{n-1}}^{t_{n}}\frac{\alpha }{\bar{\eta }}\Vert \nabla e(\eta (t))\Vert ^2dt\nonumber \\&\qquad +\,C\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt\right. +\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{1}{48} k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}\nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla (\widehat{U}_h(t)-U_h(t))\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\hat{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt \right\} ,\nonumber \\ \end{aligned}$$

(4.30)

where $C=\max \left\{ \frac{10\alpha C^2_\Pi }{\alpha ^2-\gamma ^2\bar{\eta }}, \frac{10\alpha C^2_P C^2_1}{\alpha ^2-\gamma ^2\bar{\eta }}, \frac{10\alpha \gamma ^2 C^2_1}{\alpha ^2-\gamma ^2\bar{\eta }},2\beta ,\frac{10\alpha C^2_1\beta ^2}{\alpha ^2-\gamma ^2\bar{\eta }}\right\}$. Summing over $n=2$ to N, and using the fact ${\hat{e}}(t_n)=e(t_n)$ for $n\ge 0$, and (4.24), we obtain the upper bound in (4.12). The lower bound in (4.12) can be obtained with the same argument as for the delay-dependent reconstruction. Thus, we complete the proof. $\square$

Remark 4.1

We first observe that it follows from (3.1) that

$$\begin{aligned} \Vert \nabla (\widetilde{U}_h(t)-U_h(t))\Vert ^2=\frac{1}{4}(t-t_{n-1})^2(t-t_{n})^2\Vert \nabla {\widetilde{W}}^{n}_h\Vert ^2,\quad \forall t\in I_n. \end{aligned}$$

(4.31)

Then in view of

$$\begin{aligned} \int _{I_n}(t-t_{n-1})^2(t-t_{n})^2=\frac{1}{30}k^5_n, \end{aligned}$$

(4.32)

the upper bound of the error can be given by

$$\begin{aligned}&\Vert e(T)\Vert ^2 +\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\le \Vert e(0)\Vert ^2+\alpha \int _{\eta (0)}^{0}\Vert \nabla e(t)\Vert ^2dt\nonumber \\&\qquad +\,C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt\right. +\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\nonumber \\&\qquad +\,\frac{1}{30}k^5_n\Vert \nabla \widetilde{W}^{n}_h\Vert ^2_{K}\left. +\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt\right\} . \end{aligned}$$

(4.33)

Further, if $e(t)=0$ for $t\in [t_*,0]$, then

$$\begin{aligned}&\Vert e(T)\Vert ^2+\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\nonumber \\&\quad \le C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}\right\} \nonumber \\&\qquad +\,C\left\{ \int _{0}^{T}h^2\Vert f-\tilde{f}\Vert ^2dt+\frac{1}{30}\sum \limits ^N_{n=1}k^5_n\Vert \nabla \widetilde{W}^{n}_h\Vert ^2+\int _{0}^{T}\Vert \nabla \widetilde{R}_U(t)\Vert ^2dt\right\} .~~~~ \end{aligned}$$

(4.34)

Some observations analogous to the delay-dependent reconstruction can be presented for the multi-point reconstruction. In view of (3.5) and (4.32), we have the following alternative estimate under the assumption (2.6)

$$\begin{aligned}&\Vert e(T)\Vert ^2 +\int _{\eta (T)}^{T}\left( \alpha \Vert \nabla e\Vert ^2+\frac{C_\alpha }{2}\Vert \nabla {\hat{e}}\Vert ^2\right) dt\le \Vert e(t_1)\Vert ^2+\alpha \int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt\nonumber \\&\qquad +\,C\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt\right. +\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}\nonumber \\&\qquad +\,\left( \frac{k_{n-1}^2k_n^3}{48}+\frac{k^5_n}{30}\right) \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}\nonumber \\&\qquad +\,\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\hat{f}\Vert ^2_{K}dt+\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla \widehat{G}_U(t)\Vert ^2_{K}dt\right\} . \end{aligned}$$

(4.35)

Remark 4.2

Under the condition (2.6), the constant C in Theorem 4.2 is independent of the final time T, and therefore the a posteriori error estimators can serve as long-time a posteriori error estimators. To the best of our knowledge, this is the first result on long-time a posteriori error estimates for numerical methods for DDEs and PDDEs. For more general case, we have the following estimates where the constant C will depend on the final time T.

4.2 A Posteriori Error Estimates on Finite Time Interval

If the condition (2.6) doesn’t hold, then we have the following a posteriori error estimates for CNFE scheme on a finite time interval.

Theorem 4.3

Let u and $U_h$ satisfy (2.1) and (2.14), respectively. Then, there is a constant C such that

$$\begin{aligned}&\Vert e(T)\Vert ^2+\frac{\alpha }{7} \int _{0}^{T}\Vert \nabla \tilde{U}_h(t)-\nabla U_h(t)\Vert ^2dt \nonumber \\&\quad \le \Vert e(T)\Vert ^2+\alpha \int _{0}^{T}\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{6}\Vert \nabla {\tilde{e}}(t)\Vert ^{2}\right) dt\nonumber \\&\quad \le C\left( \Vert e(0)\Vert ^2+\int _{\eta (0)}^{0}\Vert \nabla e(t)\Vert ^2dt+\int _{0}^{T}\Vert \nabla U_h(t)-\nabla {\widetilde{U}}_h(t)\Vert ^2dt\right) \nonumber \\&\qquad +\,C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S^2_{K,n}dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\tilde{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} . \end{aligned}$$

(4.36)

Proof

Similar to (4.21), from (4.18) we use the Young’s inequality and Proposition 2.2 to obtain

$$\begin{aligned}&\frac{1}{2}\frac{d}{dt}\Vert {\tilde{e}}(t)\Vert ^2+\frac{\alpha }{2}\Vert \nabla e\Vert ^2+\frac{\alpha }{2}\Vert \nabla {\tilde{e}}\Vert ^2\le \frac{\alpha }{12}\Vert \nabla {\tilde{e}}\Vert ^2+\frac{3\gamma ^2}{\alpha }\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ C_\Pi S_{K,n}\Vert \nabla {\tilde{e}}\Vert _{K}\right. +C_\Pi h_K\left| t-t_{n-1/2}\right| \Vert {\widetilde{W}}^{n}_h\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K}\nonumber \\&\qquad +\,C_PC_1h_K\Vert f-{\tilde{f}}\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K}+\beta \Vert \nabla ({\widetilde{U}}_h(t)-U_h(t))\Vert ^2_{K}\nonumber \\&\qquad +\,\left. C_1\gamma \Vert {\widetilde{R}}_U(t)\Vert _{K}\Vert \nabla {\tilde{e}}\Vert _{K} \right\} . \end{aligned}$$

(4.37)

An application of the Young’s inequality yields

$$\begin{aligned}&\frac{d}{dt}\Vert {\tilde{e}}(t)\Vert ^2+\alpha \Vert \nabla e\Vert ^2+\alpha \Vert \nabla {\tilde{e}}\Vert ^2\le \frac{\alpha }{6}\Vert \nabla {\tilde{e}}\Vert ^2+\frac{6\gamma ^2}{\alpha }\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\qquad +\,\sum \limits _{K\in \mathcal T_h}\left\{ \epsilon C_\Pi ^2 S_{K,n}^2\right. +\epsilon C_\Pi ^2h^2_{K}\left( t-t_{n-1/2}\right) ^2 \Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\frac{4}{\epsilon }\Vert \nabla {\tilde{e}}\Vert ^2_{K}\nonumber \\&\qquad +\,\epsilon C_P^2C_1^2h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}+2\beta \Vert \nabla ({\widetilde{U}}_h(t)-U_h(t))\Vert ^2_{K}\nonumber \\&\qquad +\,\left. \epsilon \gamma ^2 C_1^2\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}\right\} . \end{aligned}$$

(4.38)

Now choosing $\epsilon =\frac{6}{\alpha }$ and integrating from $t_{n-1}$ to $t_{n}$, we obtain

$$\begin{aligned}&\Vert {\tilde{e}}(t_{n})\Vert ^2+\alpha \int _{t_{n-1}}^{t_{n}}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\le \Vert {\tilde{e}}(t_{n-1})\Vert ^2+\frac{6\gamma ^2}{\alpha }\int _{t_{n-1}}^{t_{n}}\Vert \nabla e(\eta (t))\Vert ^2dt\nonumber \\&\qquad +\,C_2\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla (\widetilde{U}_h(t)-U_h(t))\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} , \end{aligned}$$

(4.39)

where $C_2=\max \left\{ \frac{6 C^2_\Pi }{\alpha }, \frac{6 C^2_P C^2_1}{\alpha }, \frac{6\gamma ^2 C^2_1}{\alpha },2\beta \right\}$. Summing over $n=1$ to N, noting ${\tilde{e}}(t_n)=e(t_n)$ for $1\le n\le N$, and ${\tilde{U}}_h(t)=U(t)$ for $t\in [t_{*},0]$, we get

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{0}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\le \Vert e(0)\Vert ^2+\frac{6\gamma ^2}{\alpha }\int _{0}^{T}\Vert \nabla e(\eta (t))\Vert ^2dt\nonumber \\&\qquad +\,C_2\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widetilde{R}}_U(t)\Vert ^2_{K}dt \right\} +C_2\int _{0}^{T}\Vert \nabla U_h(t)-\nabla \widetilde{U}_h(t)\Vert ^2dt. \end{aligned}$$

(4.40)

Using (4.24), by induction, we have

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{0}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\le \Vert e(0)\Vert ^2+\frac{6\gamma ^2\bar{\eta }}{\alpha ^2}\alpha \int _{\eta (0)}^{\eta (T)}\Vert \nabla e(t)\Vert ^2dt\nonumber \\&\qquad +\,C_2\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\tilde{f}}\Vert ^2_{K}dt\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widetilde{R}}_U(t)\Vert ^2_{K}dt \right\} +C_2\int _{0}^{T}\Vert \nabla U_h(t)-\nabla {\widetilde{U}}_h(t)\Vert ^2dt\nonumber \\&\quad \le C(T)\left( \Vert e(0)\Vert ^2+\frac{6\gamma ^2\bar{\eta }}{\alpha }\int _{\eta (0)}^{0}\Vert \nabla e(t)\Vert ^2dt+C_2\int _{0}^{T}\Vert \nabla U_h(t)-\nabla {\widetilde{U}}_h(t)\Vert ^2dt\right) \nonumber \\&\qquad +\,C(T)C_2\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\tilde{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} , \end{aligned}$$

(4.41)

where

$$\begin{aligned} C(T)=\left\{ \begin{array}{ll} M,\quad &{}\quad {\hbox {if}}\quad 6\gamma ^2\bar{\eta }=\alpha ^2,\\ \frac{\left( 6\gamma ^2\bar{\eta }\right) ^M-\alpha ^{2M}}{6\gamma ^2\bar{\eta }\alpha ^{2M-2}-\alpha ^{2M}},&{}\quad {\hbox {if}}\quad 6\gamma ^2\bar{\eta }\not =\alpha ^2. \end{array} \right. \end{aligned}$$

Note that M, which depends on T, has been defined in Sect. 2.2. Then the upper bound in Theorem 4.3 is proved with $C=C(T)\max \{1,C_2,\frac{6\gamma ^2\bar{\eta }}{\alpha }\}$.

We now turn to establish the lower bound. It follows from (4.25) that

$$\begin{aligned} \Vert \nabla \tilde{U}_h(t)-\nabla U_h(t)\Vert ^2&\le 7\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{6}\Vert \nabla {\tilde{e}}(t)\Vert ^{2}\right) . \end{aligned}$$

(4.42)

As a consequence, it holds that

$$\begin{aligned}&\Vert e(T)\Vert ^2+\frac{\alpha }{7} \int _{0}^{T}\Vert \nabla \tilde{U}_h(t)-\nabla U_h(t)\Vert ^2dt \nonumber \\&\qquad \le \Vert e(T)\Vert ^2+\alpha \int _{0}^{T}\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{6}\Vert \nabla {\tilde{e}}(t)\Vert ^{2}\right) dt, \end{aligned}$$

(4.43)

which gives the lower bound and completes the proof. $\square$

Remark 4.3

In view of (4.31) and (4.32), $\int _{0}^{T}\Vert \nabla \tilde{U}_h(t)-\nabla U_h(t)\Vert ^2dt$ in the lower and upper bounds of Theorem 4.3 can be replaced by $\frac{1}{30}\sum \limits _{n=1}^Nk^5_n\Vert \nabla {\widetilde{W}}^n_h\Vert .$ Moreover, the upper bound can be also estimated by

$$\begin{aligned}&\Vert e(T)\Vert ^2 +\alpha \int _{0}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\nonumber \\&\quad \le C\left( \Vert e(0)\Vert ^2+\int _{\eta (0)}^{0}\Vert \nabla e(t)\Vert ^2dt+\frac{1}{30}\sum \limits _{n=1}^Nk^5_n\Vert \nabla {\widetilde{W}}^n_h\Vert \right) \nonumber \\&\qquad +\,C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\tilde{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} , \end{aligned}$$

(4.44)

and, if $e(t)=0$ for $t\in [t_*,0]$, then

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{0}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\tilde{e}}\Vert ^2\right) dt\nonumber \\&\quad \le \frac{C}{30}\sum \limits _{n=1}^Nk^5_n\Vert \nabla {\widetilde{W}}^n_h\Vert +C\sum \limits ^N_{n=1}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widetilde{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\tilde{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla \widetilde{R}_U(t)\Vert ^2_{K}dt \right\} . \end{aligned}$$

(4.45)

If the multi-point reconstruction ${\widehat{U}}_h$ is explored, then we have the following estimate of CNFE scheme for general delay problems.

Theorem 4.4

Let u and $U_h$ satisfy (2.1) and (2.14), respectively. Then, there is a constant C such that

$$\begin{aligned}&\Vert e(T)\Vert ^2 +\frac{\alpha }{8} \int _{t_1}^{T}\Vert \nabla \widehat{U}_h(t)-\nabla U_h(t)\Vert ^2dt \nonumber \\&\quad \le \Vert e(T)\Vert ^2+\alpha \int _{t_1}^{T}\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{7}\Vert \nabla {\hat{e}}(t)\Vert ^{2}\right) dt\nonumber \\&\quad \le C\left( \Vert e(\cdot ,t_1)\Vert ^2+\int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt+\int _{t_1}^{T}\Vert \nabla U_h(t)-\nabla {\widehat{U}}_h(t)\Vert ^2dt\right) \nonumber \\&\qquad +\,C\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{1}{48} k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\hat{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt \right\} .~~~~ \end{aligned}$$

(4.46)

Proof

We use the Young’s inequality and Proposition 2.2 to obtain the following inequality similar to (4.28),

$$\begin{aligned}&\frac{1}{2}\frac{d}{dt}\Vert {\hat{e}}(t)\Vert ^2 +\frac{\alpha }{2}\Vert \nabla e\Vert ^2+\frac{\alpha }{2}\Vert \nabla {\hat{e}}\Vert ^2\le \frac{\alpha }{14}\Vert \nabla {\hat{e}}\Vert ^2+\frac{7\gamma ^2}{2\alpha }\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\quad +\sum \limits _{K\in \mathcal T_h}\left\{ C_\Pi S_{K,n}\Vert \nabla {\hat{e}}\Vert _{K} \right. +C_\Pi \left| t-t_{n-1/2}\right| \Vert {\widehat{W}}^{n}_h\Vert _{K}\Vert \nabla {\hat{e}}\Vert _{K} \nonumber \\&\quad +\,C_PC_1h_K\Vert f-{\hat{f}}\Vert _{K}\Vert \nabla {\hat{e}}\Vert _{K}+\frac{C_1\beta k_{n-1}|t-t_{n-1/2}|}{2}\Vert \nabla {\widehat{W}}^{n}_h\Vert _K\Vert \nabla {\hat{e}}\Vert _K\nonumber \\&\quad +\left. \beta \Vert \nabla ({\widehat{U}}_h(t)-U_h(t))\Vert ^2_{K}+C_1\gamma \Vert {\widehat{R}}_U(t)\Vert _{K}\Vert \nabla {\hat{e}}\Vert _{K} \right\} . \end{aligned}$$

(4.47)

Using the Young’s inequality again yields

$$\begin{aligned}&\frac{d}{dt}\Vert {\hat{e}}(t)\Vert ^2+\alpha \Vert \nabla e\Vert ^2+\alpha \Vert \nabla {\hat{e}}\Vert ^2\le \frac{\alpha }{7}\Vert \nabla {\hat{e}}\Vert ^2+\frac{7\gamma ^2}{\alpha }\Vert \nabla e(\eta (t))\Vert ^2\nonumber \\&\quad +\,\sum \limits _{K\in \mathcal T_h}\left\{ \epsilon C_\Pi ^2 S_{K,n}^2\right. +\epsilon C_\Pi ^2h^2_{K}\left( t-t_{n-1/2}\right) ^2 \Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{5}{\epsilon }\Vert \nabla {\hat{e}}\Vert ^2_{K}\nonumber \\&\quad +\,\epsilon C_P^2C_1^2h_K^2\Vert f-{\hat{f}}\Vert ^2_{K}+\frac{\epsilon C_1^2\beta ^2}{4} k^2_{n-1}(t-t_{n-1/2})^2\Vert \nabla {\widehat{W}}^{n}_h\Vert _K\nonumber \\&\quad +\,\left. 2\beta \Vert \nabla ({\widehat{U}}_h(t)-U_h(t))\Vert ^2_{K}+\epsilon \gamma ^2 C_1^2\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}\right\} . \end{aligned}$$

(4.48)

Now choosing $\epsilon =\frac{7}{\alpha }$, integrating from $t_{n-1}$ to $t_{n}$, summing over $n=2$ to N, and noting $\hat{e}(t_n)=e(t_n)$ for $2\le n\le N$, we get

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{t_1}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{7}\Vert \nabla {\hat{e}}\Vert ^2\right) dt\nonumber \\&\quad \le \Vert e(0)\Vert ^2+\frac{7\gamma ^2}{\alpha }\int _{t_1}^{T}\Vert \nabla e(\eta (t))\Vert ^2dt+C_3\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt\right. \nonumber \\&\qquad +\,\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{1}{48} k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}+\int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-{\hat{f}}\Vert ^2_{K}dt\nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt \right\} +C_3\int _{t_1}^{T}\Vert \nabla U_h(t)-\nabla \widehat{U}_h(t)\Vert ^2dt, \end{aligned}$$

(4.49)

where $C_3=\max \left\{ \frac{7 C^2_\Pi }{\alpha }, \frac{7 C^2_P C^2_1}{\alpha }, \frac{7 \gamma ^2C^2_1}{\alpha },\frac{7 \beta ^2C^2_1}{\alpha }, 2\beta \right\}$. By induction, it follows that

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{t_1}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{6}\Vert \nabla {\hat{e}}\Vert ^2\right) dt\nonumber \\&\quad \le C(T)\left( \Vert e(t_1)\Vert ^2+\int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt+C_3\int _{t_1}^{T}\Vert \nabla U_h(t)-\nabla {\widehat{U}}_h(t)\Vert ^2dt\right) \nonumber \\&\qquad +\,C(T)C_3\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt+\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}\right. \nonumber \\&\qquad +\,\frac{1}{48} k_{n-1}^2k_n^3 \Vert \nabla \widehat{W}^{n}_h\Vert ^2_{K}+\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\hat{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt \right\} ,~~~~ \end{aligned}$$

(4.50)

where

$$\begin{aligned} C(T)=\left\{ \begin{array}{ll} M-1,\qquad &{}{\hbox {if}}\quad 7\gamma ^2\bar{\eta }=\alpha ^2,\\ \frac{\left( 7\gamma ^2\bar{\eta }\right) ^{M-1}-\alpha ^{2M-2}}{7\gamma ^2\bar{\eta }\alpha ^{2M-4}-\alpha ^{2M-2}},&{}{\hbox {if}}\quad 7\gamma ^2\bar{\eta }\not =\alpha ^2. \end{array} \right. \end{aligned}$$

Then the upper bound in Theorem 4.4 is obtained with $C=C(T)\max \{1,C_3, \frac{7\gamma ^2\bar{\eta }}{\alpha }\}$.

As for the lower bound, similar to the case of the delay-dependent reconstruction, we have

$$\begin{aligned} \int _{t_1}^{T}\Vert \nabla \widehat{U}_h(t)-\nabla U_h(t)\Vert ^2dt&\le 8\int _{t_1}^{T}\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{7}\Vert \nabla \hat{e}(t)\Vert ^{2}\right) dt, \end{aligned}$$

(4.51)

and therefore

$$\begin{aligned}&\Vert e(T)\Vert ^2+\frac{\alpha }{8} \int _{t_1}^{T}\Vert \nabla \widehat{U}_h(t)-\nabla U_h(t)\Vert ^2dt\nonumber \\&\quad \le \Vert e(T)\Vert ^2+\alpha \int _{t_1}^{T}\left( \Vert \nabla e(t)\Vert ^{2}+\frac{1}{7}\Vert \nabla {\hat{e}}(t)\Vert ^{2}\right) dt. \end{aligned}$$

(4.52)

This gives the lower bound and completes the proof. $\square$

Remark 4.4

We first make an analogous observation to Remark 4.3, that is, $\int _{0}^{T}\Vert \nabla \widehat{U}_h(t)-\nabla U_h(t)\Vert ^2dt$ in the lower and upper bounds of Theorem 4.4 can be replaced by $\frac{1}{30}\sum \nolimits _{n=1}^Nk^5_n\Vert \nabla {\widehat{W}}^n_h\Vert$. Moreover, the upper bound can also be estimated by

$$\begin{aligned}&\Vert e(T)\Vert ^2+\alpha \int _{t_1}^{T}\left( \Vert \nabla e\Vert ^2+\frac{1}{7}\Vert \nabla {\hat{e}}\Vert ^2\right) dt\nonumber \\&\quad \le C\left( \Vert e(t_1)\Vert ^2+\int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt+\frac{1}{30}\sum \limits _{n=1}^Nk^5_n\Vert \nabla \widehat{W}^n_h\Vert \right) \nonumber \\&\qquad +\,C\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}\left\{ \int _{t_{n-1}}^{t_{n}}S_{K,n}^2dt\right. +\frac{1}{12} k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}+\frac{1}{48}k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K}\nonumber \\&\qquad +\,\left. \int _{t_{n-1}}^{t_{n}}h_K^2\Vert f-\hat{f}\Vert ^2_{K}dt+\int _{t_{n-1}}^{t_{n}}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt \right\} . \end{aligned}$$

(4.53)

Remark 4.5

As already mentioned in Remark 3.1, in this paper we also consider the reconstruction $\widehat{W}^{1}_h={\widetilde{W}}^1_h$, although the reconstruction (3.7) is natural. Due to the presence of the term $\Vert e(t_1)\Vert ^2$ and $\int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt$ on the right-hand sides of (4.12), (4.35) and (4.46), it seems that the estimators derived in (4.12), (4.35) and (4.46) are not sufficient to provide three meaningful a posteriori error estimators. Nevertheless, the optimality of the a posteriori error estimates in (4.12), (4.35) and (4.46) are justified by using the estimates of $\Vert e(t_1)\Vert ^2$ and $\int _{\eta (t_1)}^{t_1}\Vert \nabla e(t)\Vert ^2dt$ from (4.11) and (4.36) in which the error bounds are of optimal order in the $L^2(0,t_1;H^1)$-norm.

5 Numerical Examples

We now illustrate the theoretical results on the error estimators for which the error comes either from the time discretization, or from the space discretization, or from both of them. Let us consider applying the CNFE method (2.14) to the following generalized diffusion equation with a delay term

$$\begin{aligned} \frac{\partial }{\partial t}u(x,t)=a\frac{\partial ^2}{\partial x^2}u(x,t)+b\frac{\partial ^2}{\partial x^2} u(x,\eta (t))+f(x,t),\quad t\in [0,T], \quad x\in [0,1]. \end{aligned}$$

(5.1)

Here the forcing term f, and the initial and boundary value conditions are chosen such that the exact solution is $u(x,t)=x(1-x)\sin t$.

Since it is natural to use the multi-point reconstruction of the numerical solution to DDEs for a posteriori error estimates, in the following numerical examples we will mainly consider the multi-point reconstruction (3.5). In our implementation, we use the uniform time partition and uniform space partition with constant stepsize k and h, respectively. Moreover, all the integral between $t_{n-1}$ and $t_n$ are approximated by the Gauss-Legendre quadrature formula with three nodes, where the integral, as a polynomial of degree 4, is integrated exactly by this formula.

Let us define the space error estimator $\mathcal E_K$

$$\begin{aligned} \mathcal E_K:=\sum \limits _{n=2}^N\sum \limits _{K\in \mathcal T_h}\int ^{t_n}_{t_{n-1}}S^2_{K,n}dt, \end{aligned}$$

where the contributions $S_{K,n}$ in spatial discretization have been introduced in Theorem 4.2. The time error estimator $\mathcal E_T$ is defined by

$$\begin{aligned} \mathcal E_T:= & {} \mathcal E_{W1}+\mathcal E_{W2}+\mathcal E_U+\mathcal E_f+\mathcal E_\eta , \end{aligned}$$

where the contributions $\mathcal E_{W1}$, $\mathcal E_{W2}$, $\mathcal E_U$, $\mathcal E_f$ and $\mathcal E_\eta$ are computed by

$$\begin{aligned} \mathcal E_{W1}:= & {} \frac{1}{12}\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K},\quad \mathcal E_{W2}:=\frac{1}{48}\sum \limits ^N_{n=2}\sum \limits _{K\in \mathcal T_h}k_{n-1}^2k_n^3 \Vert \nabla {\widehat{W}}^{n}_h\Vert ^2_{K},\\ \mathcal E_U:= & {} \sum \limits _{n=2}^N \sum \limits _{K\in \mathcal T_h}\int _{t_{n-1}}^{t_n}\Vert \nabla ({\widehat{U}}_h(t)-U_h(t))\Vert ^2_{K}dt,\\ {\mathcal {E}}_f:= & {} \sum \limits _{n=2}^N \sum \limits _{K\in \mathcal T_h}\int _{t_{n-1}}^{t_n}h_K^2\Vert f-{\hat{f}}\Vert ^2_{K}dt,\quad {\mathcal {E}}_\eta :=\sum \limits _{n=2}^N \sum \limits _{K\in \mathcal T_h}\int _{t_{n-1}}^{t_n}\Vert \nabla {\widehat{R}}_U(t)\Vert ^2_{K}dt, \end{aligned}$$

respectively.

5.1 Example 1: Linear Proportional Delay

In this case, we test the long time a posteriori error estimates (4.12). For this purpose, we consider linear proportional delay $\eta (t)=0.5t-1$ and take $a=10$, $b=5$, and $T=20$ in problem (5.1). It is easy to verify the condition (2.6) which implies that the system is stable.

Let us denote by ${\hbox {Err}}_T$ the square of the error at time T, i.e.,

$$\begin{aligned} {\hbox {Err}}_T=\Vert e(T)\Vert ^2, \end{aligned}$$

by ${\hbox {Err}}_{\eta 1}$ the square of the $L^2(\eta (T),T;V)$-norm of the error, i.e.,

$$\begin{aligned} {\hbox {Err}}_{1\eta }=\alpha \int _{\eta (T)}^{T}\Vert \nabla e\Vert ^2dt, \end{aligned}$$

The true errors ${\hbox {Err}}_T$, ${\hbox {Err}}_{1\eta }$, and the a posteriori error estimators $\mathcal E_K$, $\mathcal E_{W1}$, $\mathcal E_{W2}$, $\mathcal E_U$, ${\mathcal {E}}_f$, ${\mathcal {E}}_\eta$, and their temporal convergence orders are listed in Tables 1 and 2, respectively. From these numerical results, we observe that the true error ${\hbox {Err}}_T$ and the a posteriori error estimators $\mathcal E_{W2}$, $\mathcal E_U$, ${\mathcal {E}}_f$, ${\mathcal {E}}_\eta$, are of optimal order 4 (since they estimate the square of the error), and the space error estimator $\mathcal E_K$ does not depend on the time stepsize k. The a posteriori quantities $\mathcal E_{W1}:=\frac{1}{12}\sum ^N_{n=2}\sum _{K\in \mathcal T_h}k_n^3 h^2_{K}\Vert {\widehat{W}}^{n}_h\Vert ^2_{K}$ is of optimal order 2. We also note that since the true error ${\hbox {Err}}_{1\eta }$ is badly affected by the space mesh diameter h, there is a continuing decreasing in its temporal convergence order.

Table 1 Example 1: The errors and their orders of CNFE method (2.14) for (5.1) with $\eta (t)=0.5t-1$, where $T=20$

Table 2 Example 1: The errors and their orders of CNFE method (2.14) for (5.1) with $\eta (t)=0.5t-1$, where $T=20$

To clearly illustrate the temporal convergence behaviours of the true errors ${\hbox {Err}}_T$, ${\hbox {Err}}_{1\eta }$, and the a posteriori error estimators $\mathcal E_K$, $\mathcal E_{W1}$, $\mathcal E_{W2}$, $\mathcal E_U$, ${\mathcal {E}}_f$, ${\mathcal {E}}_\eta$, we also present them in Fig. 1 against 1/k in a log–log scale.

In order to measure the quality of our estimators, the estimated error is compared to the true error by introducing the so-called effecitvity index. We define the effectivity index $ei_U$ of the upper error bound as

$$\begin{aligned} ei_U:=\left( \frac{\mathcal E_K+\mathcal E_T}{{\hbox {Err}}_T+{\hbox {Err}}_{1\eta }}\right) ^{\frac{1}{2}}. \end{aligned}$$

To measure the contribution of space discretization error and time discretization error, we also define the following effectivity indices in space and in time

$$\begin{aligned} ei^K_U:=\left( \frac{{\mathcal {E}}_K}{{\hbox {Err}}_T+{\hbox {Err}}_{1\eta }}\right) ^{\frac{1}{2}} \quad {\hbox {and}}\quad ei^T_U:=\left( \frac{\mathcal E_T}{{\hbox {Err}}_T+{\hbox {Err}}_{1\eta }}\right) ^{\frac{1}{2}}. \end{aligned}$$

respectively.

Referring to Table 3, we observe that for a fixed space mesh diameter h, the computing error is mainly due to the space discretizaiton and the space error estimator $\mathcal E_K$ behaves as the true error. To clearly illustrate the effectivity indices $ei_U^K$ and $ei_U^T$, we change simultaneously the space mesh diameter h and the time stepsize k with h being proportional to $k^2$, which is the natural choice because of the error behaviour $O(h+k^2)$. The numerical results are presented in Table 4. The correct temporal and spatial convergence orders are observed for the true errors. We can see that the effectivity indices $ei_U$, $ei_U^K$ and $ei_U^T$ presented in Table 4 for CNFE method seem to be asymptotically constant (around 15.4, 15.4 and 0.3, respectively).

Fig. 1

Log–log plot of the errors of CNFE method (2.14) for (5.1) with $\eta (t)=0.5t-1$ in Example 1, where $T=20$ and $h=0.000125$

Table 3 Example 1: The effecitvity indices of error estimators of CNFE method (2.14) for (5.1) with $\eta (t)=0.5t-1$, where $T=20$

Table 4 Example 1: The effecitvity indices of error estimators of CNFE method (2.14) for (5.1) with $\eta (t)=0.5t-1$, where $T=20$

5.2 Example 2: Nonlinear Delay

In this case, we consider nonlinear delay $\eta (t)=0.8t-\ln (t+2)$ and take $a=1$ and $b=2$ in problem (5.1). It is easy to verify that the condition (2.6) is not satisfied. Thus we consider the a posteriori error estimates (4.46) on a finite time interval [0, T] with $T=1$. On account of the nonlinearity of the delay, the error behaviours are certainly much more complicated compared to the previous example. We denote by ${\hbox {Err}}_1$ the square of the $L^2(t_1,T;V)$-norm of the error, i.e.,

$$\begin{aligned} {\hbox {Err}}_1=\alpha \int _{t_1}^{T}\Vert \nabla e\Vert ^2dt. \end{aligned}$$

The true errors ${\hbox {Err}}_T$, ${\hbox {Err}}_{1}$, and the a posteriori error estimators $\mathcal E_K$, $\mathcal E_{W1}$, $\mathcal E_{W2}$, $\mathcal E_U$, ${\mathcal {E}}_f$, ${\mathcal {E}}_\eta$, and their temporal convergence orders are listed in Tables 5 and 6, respectively. The temporal convergence behaviours of these quantities are also presented in Fig. 2 against 1/k in a log–log scale. From these numerical results, we still observe the correct convergence orders for all these quantities, though the convergence orders are slightly affected by the nonlinearity of the delay and the space mesh diameter h.

Table 5 Example 2: The errors and their orders of CNFE method (2.14) for (5.1) with $\eta (t)=0.8t-\ln (t+2)$, where $T=1$

Table 6 Example 2: The errors and their orders of CNFE method (2.14) for (5.1) with $\eta (t)=0.8t-\ln (t+2)$, where $T=1$

Fig. 2

Log–log plot of the errors of CNFE method (2.14) for (5.1) with $\eta (t)=0.8t-\ln (t+2)$ in Example 2, where $T=1$, $h=0.000125$

In view of the a posteriori error estimates (4.46), the lower and upper estimators are defined by $\displaystyle \frac{\alpha }{8} \mathcal E_U$ and $\mathcal E_K+\mathcal E_T,$ respectively. We are interested in computing the effectivity indices $ei_L$ and $ei_U$, defined as

$$\begin{aligned} ei_L:=\left( \frac{\hbox {Lower~estimator}}{\frac{8}{7}{\hbox {Err}}_1}\right) ^{\frac{1}{2}} \quad \hbox {and}\quad ei_U:=\left( \frac{\hbox {Upper~estimator}}{{\hbox {Err}}_T+\frac{8}{7}{\hbox {Err}}_1}\right) ^{\frac{1}{2}}, \end{aligned}$$

respectively. These effectivity indices are presented in Table 7 for a fixed space mesh diameter h and in Table 8 for changeable h with simultaneous decrease of the time stepsize k. From Table 8, we find that the effectivity index $ei_L$ of the lower estimator of the CNFE method is around 0.18 and the effectivity index $ei_U^T$ of the upper time error estimator $\mathcal E_T$ is around 0.66 when the space mesh diameter h and the time stepsize k are changed simultaneously with h being proportional to $k^2$.

Table 7 Example 2: The effecitvity indices of error estimators of CNFE method (2.14) for (5.1) with $\eta (t)=0.8t-\ln (t+2)$, where $T=1$

Table 8 Example 2: The effecitvity indices of error estimators of CNFE method (2.14) for (5.1) with $\eta (t)=0.8t-\ln (t+2)$, where $T=1$

6 Concluding Remarks

Seeing that DDEs and PDDEs and their variants play a pivotal role in the modeling of several physical and biological phenomena, and that the solution to this class of equations is generally not sufficiently smooth and has breaking points, a posteriori error estimates are extremely important for numerically solving them. In this work we presented the a posteriori error analysis for the CNFE method (2.14) for generalized diffusion equation with delay. We first obtained the sufficient condition for the stability of this class of equations. Based on this stability condition, we proved long-time a posteriori error estimates for stable PDDEs for the first time, after we introduced two different continuous and piecewise quadratic reconstructions to estimate the time discretization error. By means of these reconstructions, we further derived a posteriori upper and lower error bounds for the CNFE method (2.14) on a finite time interval for general systems. These lower bounds are optimal in a certain sense. Although many researchers have investigated the stability and a priori error estimates of some numerical methods for PDDEs and their variant, it is the first time to explore a posteriori error analysis for fully discrete numerical methods for such types of equations. It is worth emphasizing that the techniques exploited in this paper can be easily extended to other types of linear PDDEs, even PDDEs with several delays.

It is noteworthy that these error bounds depend only upon the discretization parameters and the data of problems, and thus they are computable. We carried out various numerical experiments for the CNFE method (2.14) for generalized diffusion equation with different delays, e.g., linear proportional delay $\eta (t)=0.5t-1$ and nonlinear delay $\eta (t)=0.8t-\ln (t+2)$. For both delay problems, these experiments exactly verify the theoretical results, and the effectivity index $ei_U^T$ of the upper time error estimator $\mathcal E_T$ stays close to a constant. We also found that the upper space error estimator $\mathcal E_K$ may slightly affect the effectivity of the upper error estimator when the CNFE method with linear element was applied to general systems with nonlinear delay. Studying a posteriori error estimates of high-order finite element methods for partial functional differential equations will be our future work.