key: cord-0705731-3n4gbtmo authors: Wang, Wansheng; Yi, Lijun; Xiao, Aiguo title: A Posteriori Error Estimates for Fully Discrete Finite Element Method for Generalized Diffusion Equation with Delay date: 2020-07-01 journal: J Sci Comput DOI: 10.1007/s10915-020-01262-5 sha: 52ba4a6657074605f67e64a394f45c739c5e0ca8 doc_id: 705731 cord_uid: 3n4gbtmo 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. with A being self-adjoint positive definite, where "∇" 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 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., τ (t) = t − η(t) = τ with τ > 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 θ -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 ξ 0 = 0 and introduce the points {ξ i } i≥1 such that Under the given assumptions, the solution of (1.1)-(1.3) on (ξ i , ξ i+1 ] follows from the solution on (ξ i−1 , ξ i ] and the solution for t ≥ 0 can be found by a "method of steps". Although u(x, t) will inherit its smoothness properties on each interval [ξ i , ξ i+1 ] from the smoothness properties of φ on [t * , 0] for smooth f and η, the solution u(x, t) may suffer regularity loss of time smoothness at points ξ i , i = 0, 1, 2, . . .; see, e.g., [5] . This implies that even if the initial function φ is sufficiently smooth, the solution to DDEs or PDDEs is generally not sufficiently smooth and has breaking points ξ i , i = 0, 1, . . .. 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 ξ i , i = 0, 1, . . .. 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 ξ 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 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 ξ i , i = 0, 1, . . ., 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. 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. 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 ω ⊂ R d , let L p (ω), p ≥ 1, be the Lebesgue spaces with the corresponding norms · L p (ω) , and W s, p (ω), s ≥ 0, p ≥ 1, be the standard Sobolev spaces with the norms · s, p,ω . For convenience, we denote by · s,ω and · ∞,ω the norms of the space H s (ω) = W s,2 (ω) and L ∞ (ω), respectively, and · ω for the usual norm in L 2 (ω). We denote by H 1 0 (ω) the space of functions in H 1 (ω) that vanish on the boundary of ω (boundary values are taken in the sense of traces), and we denote by ·, · ω either the inner product in L 2 (ω) or the duality pairing between H 1 0 (ω) and its duality H −1 (ω). The norm of the space H −1 (ω) will be denoted by · −1,ω . Whenever ω = , we omit the letter ω in the subscripts of · s, p,ω , · s,ω , · ∞,ω , · ω , ·, · ω , and · −1,ω . We also introduce the Bochner spaces L 2 (a, b; X ) which are defined by To state the abstract weak formulation of problem (1.1), we need to define two bilinear forms, a(·, ·) and b(·, ·), corresponding to the elliptic operator A and B, respectively. They are defined respectively on H 1 With the above notations and assumptions in hand, the weak formulation of problem (1.1)-(1.3) may be stated as follows: Suppose f ∈ L 2 (0, T ; H −1 ( )) and φ ∈ L 2 (t * , 0; H 1 0 ( )), and seek u : Here we assume that the bilinear form a(·, ·) is coercive and continuous on H 1 0 ( ), i.e., with α, β ∈ R + . We also assume that the bilinear form b(·, ·) is continuous on H 1 0 ( ), i.e., |b(ϕ, χ)| ≤ γ ∇ϕ ∇χ , ∀ϕ, χ ∈ H 1 0 ( ), (2.4) with γ ∈ R + . 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 η(t) is continuously differentiable and defineη Further, let us define η (0) (t) = t, η (1) (t) = η(t), and recursively We assume that there exists a positive integer M such that η (M) (T ) ∈ [t * , 0]. Wang et. al. [47] illustrated the existence of such a positive integer M for three lag functions η(t) = t −1, 0.5t −1, t − √ t + 1. For the three lag functions, we haveη = 1, 2, 2, respectively. Here we give another nontrivial example of a lag function η satisfying condition η (M) (T ) ∈ [t * , 0]: η(t) = pt − ln(t + 2), 0.5 < p < 1. Obviously, for any fixed T > 0, there exists a positive integer M such that η (M) (T ) ∈ [− ln 2, 0]. It is also easy to getη = 2 2 p−1 for this case. Now we are in position to give the stability result. Theorem 2.1 (Stability) Suppose f ∈ L 2 (0, T ; H −1 ( )) and φ ∈ L 2 (t * , 0; H 1 0 ( )), and the bilinear forms a(·, ·) and b(·, ·) satisfy (2.3) and (2.4), respectively. Then under the condition which, in view of (2.3) and (2.4), implies that Integrate from 0 to T to get and therefore (2.7) holds. This completes the proof. We note that for constant delay equation, i.e., η(t) = t − τ (in this case,η = 1), Tian in [41] obtained a similar condition to (2.6) for stability. For any 0 < h < 1, let T h be a shape regular, conforming triangulation of into triangles K with diameter h K ≤ h. We define V h by the usual finite element space of continuous, piecewise linear functions on T h : where P 1 (K ) denotes the space of polynomials of degree at most 1. We now set The assumptions on the mesh and discrete space above ensure the existence of an interpolation operator h : Further, the following approximation property is satisfied: where the constants C 1 and C depend only on the shape regularity of the family of triangulations 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 < · · · < 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, . . . , N , we define It is noteworthy that the method (2.14) itself needs the continuous approximation U h (t). Now we shall introduce quadratic time reconstructions, which are appropriate for Eq. (2.1) corresponding method (2.14), and continuous time approximation (2.15). For 1 ≤ n ≤ N , we define the continuous, piecewise quadratic delay-dependent reconstruction U h : Here we define∂ η U n h as∂ Then in order to estimate the interpolation error in approximating the delay term, we set and hence R U (t) = 0 because η is a strictly monotone increasing function on [0, T ]. We observe that the reconstruction introduced in the above subsection depends on the past values U h (η(t)), and W n h in (3.1) is formally an approximation of ∂ f /∂t − A∂u/∂t − B∂u(η(t))/∂t in the time slab I n . Since the latter is equal to ∂ 2 u/∂t 2 , it seems natural to try to replace W n h in (3.1) by a finite difference approximation of ∂ 2 u/∂t 2 . We introduce thus a continuous, piecewise multi-point delay-independent reconstruction U h of U h in time defined for all t ∈ I n , 1 ≤ n ≤ N , by An argument similar to [35] ensures that U h (t) vanishes on the boundary so that and η(t n−1/2 ) may not belong to the same interval I i , i ≤ 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 Considering the delay property of the problem, it is natural to use the multi-point reconstruction. For t * < 0, we can define and Nevertheless, we only consider the case of W 1 h = W 1 h for simplicity in this paper. Based on the reconstructions U h and 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. We first derive long-time a posteriori error estimates in the L 2 (η(T ), T ; H 1 ) norm. For this purpose, we set e := u − U h and e := u − U h for the reconstructions U h and U h , and need to prove the following lemma. Proof For ϕ h ∈ V 0 h and 1 ≤ n ≤ N , (2.14) and (2.15) can yields Substitute the above relation into (4.3) to obtain Then, by the definition of W n h , we give the desired identity (4.1). To derive (4.2) involving the reconstruction U h , from (3.5) we get Replace n by n − 1 such that (2.14) becomes Using the relation It suffices now to substitute (4.10) on the right-hand side of (4.7) to obtain (4.2). Now we state the following long-time a posteriori error estimates. and, when the multi-point reconstruction U h is employed, where the contributions S K ,n in space discretization are defined on each triangle K of T h and each time interval I n as (4.13) Here and in what follows, [·] denotes the jump of the bracketed quantity across an interval edge, [·] = 0 for an edge on the boundary ∂ and n is the unit edge normal. Proof We first show the error estimate (4.11) for the delay-dependent reconstruction U h . For t ∈ I n , 1 ≤ n ≤ N , we use the weak formulation (2.1) and (4.4) to get ∂ tẽ (t), ϕ + a(e(t), ϕ) + b(e(η(t)), ϕ) Choosing ϕ =ẽ in the above error equation and using (4.4) again, we have (4.15) By the definition ofR U [see (3.4) ], we obtain . (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 T h . Then, in view of the continuity of the bilinear forms b(·, ·), using the Cauchy-Schwarz inequality and the Poincaré inequality, we obtain together with the coercivity of the bilinear form a(·, ·) and the relation we use the Young's inequality and Proposition 2.2 to obtain An application of the Young's inequality yields Since α 2 > γ 2η , we can choose = 8α α 2 −γ 2η . Integrate from t n−1 to t n to obtain lead to the upper bound estimate in (4.11). We now turn to establish the lower bound of errors. In view of and therefore 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 U h . For t ∈ I n , 2 ≤ n ≤ N , using (4.2) for ∂ t U h , hê + a(U h , hê ) and (3.6), we have ∂ tê (t),ê + a(e(t),ê) + b(e(η(t)),ê) Similar to the case of the delay-dependent reconstruction, it holds that and therefore d dt We choose = 10α α 2 −γ 2η and integrate from t n−1 to t n to obtain where C = max 10αC 2 . Summing over n = 2 to N , and using the factê(t n ) = e(t n ) for n ≥ 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. Remark 4. 1 We first observe that it follows from (3.1) that (4.31) Then in view of the upper bound of the error can be given by 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) 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. Proof Similar to (4.21), from (4.18) we use the Young's inequality and Proposition 2.2 to obtain An application of the Young's inequality yields Now choosing = 6 α and integrating from t n−1 to t n , we obtain α , 2β . Summing over n = 1 to N , notingẽ(t n ) = e(t n ) Using (4.24), by induction, we have where 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 , 6γ 2η α }. We now turn to establish the lower bound. It follows from (4.25) that As a consequence, it holds that (4.43) which gives the lower bound and completes the proof. (4.44) and, if e(t) = 0 for t ∈ [t * , 0], then If the multi-point reconstruction U h is explored, then we have the following estimate of CNFE scheme for general delay problems. (2.1) and (2.14), respectively. Then, there is a constant C such that Proof We use the Young's inequality and Proposition 2.2 to obtain the following inequality similar to (4.28), Using the Young's inequality again yields Now choosing = 7 α , integrating from t n−1 to t n , summing over n = 2 to N , and notinĝ e(t n ) = e(t n ) for 2 ≤ n ≤ N , we get where Then the upper bound in Theorem 4.4 is obtained with C = C(T ) max{1, C 3 , 7γ 2η α }. As for the lower bound, similar to the case of the delay-dependent reconstruction, we have and therefore This gives the lower bound and completes the proof. (4.53) As already mentioned in Remark 3.1, in this paper we also consider the reconstruction W 1 h = W 1 h , although the reconstruction (3.7) is natural. Due to the presence of the term e(t 1 ) 2 and t 1 η(t 1 ) ∇e(t) 2 dt 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 e(t 1 ) 2 and t 1 η(t 1 ) ∇e(t) 2 dt 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. 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 Here the forcing term f , and the initial and boundary value conditions are chosen such that the exact solution is u( 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 E K where the contributions S K ,n in spatial discretization have been introduced in Theorem 4.2. The time error estimator E T is defined by respectively. In this case, we test the long time a posteriori error estimates (4.12) . For this purpose, we consider linear proportional delay η(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 Err T the square of the error at time T , i.e., by Err η1 the square of the L 2 (η(T ), T ; V )-norm of the error, i.e., The true errors Err T , Err 1η , and the a posteriori error estimators and their temporal convergence orders are listed in Tables 1 and 2 , respectively. From these numerical results, we observe that the true error Err T and the a posteriori error estimators E W 2 , E U , E f , E η , are of optimal order 4 (since they estimate the square of the error), and the space error estimator E K does not depend on the time stepsize k. The a posteriori quantities K is of optimal order 2. We also note that since the true error Err 1η is badly affected by the space mesh diameter h, there is a continuing decreasing in its temporal convergence order. To clearly illustrate the temporal convergence behaviours of the true errors Err T , Err 1η , and the a posteriori error estimators E K , E W 1 , E W 2 , E U , E f , E η , 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 To measure the contribution of space discretization error and time discretization error, we also define the following effectivity indices in space and in time 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 E K behaves as the true error. To clearly illustrate the effectivity indices ei K U and ei T U , 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 K U and ei T U presented in Table 4 for CNFE method seem to be asymptotically constant (around 15.4, 15.4 and 0.3, respectively). In this case, we consider nonlinear delay η(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 Err 1 the square of the L 2 (t 1 , T ; V )-norm of the error, i.e., The true errors Err T , Err 1 , and the a posteriori error estimators E K , E W 1 , E W 2 , E U , E f , E η , 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. In view of the a posteriori error estimates (4.46), the lower and upper estimators are defined by α 8 E U and E K + E T , respectively. We are interested in computing the effectivity indices , 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 T U of the upper time error estimator 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 . 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 η(t) = 0.5t − 1 and nonlinear delay η(t) = 0.8t − ln(t + 2). For both delay problems, these experiments exactly verify the theoretical results, and the effectivity index ei T U of the upper time error estimator E T stays close to a constant. We also found that the upper space error estimator 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. A posteriori error estimates for the Crank-Nicolson method for parabolic equations A report on the use of delay differential equations in numerical modelling in the biosciences Discontinuous solutions of neutral delay differential equations A posteriori error control for fully discrete Crank-Nicolson schemes Numerical Methods for Delay Differential Equations Recent trends in the numerical solution of retarded functional differential equations. Acta Numer Convergence and stability analysis of the θ -method for delayed diffusion mathematical models The Mathematical Theory of Finite Element Methods Collocation Methods for Volterra Integral and Related Functional Equations A time delay dynamical model for outbreak of 2019-nCoV and the parameter identification A delay differential equation solver based on a continuous Runge-Kutta method with defect control Dissipativity of θ -methods for nonlinear delay differential equations of neutral type Numerical solutions of diffusion methematical models with delay Convergence of two implicit numerical schemes for diffusion mathematical models with delay Diffusion and hereditary effects in a class of population models On the asymptotic stability properties of Runge-Kutta methods for delay differential equations Asymptotic stability barriers for natural Runge-Kutta processes for delay equations Open issues in devising software for the numerical solution of implicit delay differential equations Computing breaking points in implicit delay differential equations Asymptotic expansions for regularized state-dependent neutral delay equations On the stability of predictor-corrector methods for parabolic equations with delay Stability analysis of numerical methods for systems of neutral delay-differential equations Stability and error analysis of one-leg methods for nonlinear delay differential equations Stability analysis of numerical methods for delay differential equations The stability of a class of Runge-Kutta methods for delay differential equations Asymptotic stability analysis of θ -methods for functional differential equations Numerical solution of neutral functional differential equations by Adams methods in divided difference form Introduction to the Theory and Applications of Functional Differential Equations Stability of Numerical Methods for Delay Differential Equations High order contractive Runge-Kutta methods for Volterra functional differential equations B-convergence theory of Runge-Kutta methods for stiff Volterra functional differential equations with infinite integration interval Numerical Analysis for Stiff Ordinary and Functional Differential Equations Convergence and asymptotic stability of Galerkin methods for linear parabolic equations with delays The stability of θ -methods in the numerical solution of delay differential equations An anisotropic error estimator for the Crank-Nicolson method: application to a parabolic problem Good behavior with respect to the stiffness in the numerical integration of retarded functional differential equations On the Crank-Nicolson anisotropic a posteriori error analysis for parabolic integro-differential equations An adaptive algorithm for the Crank-Nicolson scheme applied to a timedependent convection-diffusion problem Analytic-numerical solutions of diffusion mathematical models with delays Finite element interpolation of non-smooth functions satisfying boundary conditions Asymptotic stability of numerical methods for linear delay parabolic differential equations A posteriori error estimates for finite element discretizations of the heat equation Stability analysis of θ -methods for nonlinear neutral functional differential equations Nonlinear stability of general linear methods for neutral delay differential equations Stability of continuous Runge-Kutta-type methods for nonlinear neutral delay-differential equations Preserving stability implicit Euler method for nonlinear Volterra and neutral functional differential equations in Banach space A posteriori error analysis for the Crank-Nicolson-Galerkin method for the reaction-diffusion equations with delay The tracking of derivative discontinuities in systems of delay differential equations Theory and Applications of Partial Functional Differential Equations Modeling and prediction for the trend of outbreak of NCP based on a time-delay dynamic system Non-linear stability and D-convergence of Runge-Kutta methods for delay differential equations Dissipativity and exponentially asymptotic stability of the solutions for nonlinear neutral functional-differential equations Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations The authors would like to thank the referee for comments and suggestions that led to improvements in the presentation of this paper. This work was supported by the National Natural Science Foundation of China (Grant Nos. 11771060, 11771298, 11671343).