An easily computable error estimator in space and time for the wave equation

We propose a cheaper version of \textit{a posteriori} error estimator from arXiv:1707.00057 for the linear second-order wave equation discretized by the Newmark scheme in time and by the finite element method in space. The new estimator preserves all the properties of the previous one (reliability, optimality on smooth solutions and quasi-uniform meshes) but no longer requires an extra computation of the Laplacian of the discrete solution on each time step.

that the work of Adjerid et al. is related to some space-time Galerkin discretizations, which are different from the time marching schemes considered in the present work as well as in other papers cited above.
In deriving our a posteriori estimates, we follow first the approach of [15]. First of all, we recognize that the Newmark method can be reinterpreted as the Crank-Nicolson discretization of the reformulation of the governing equation as the first-order system, as in [5]. We then use the techniques stemming from a posteriori error analysis for the Crank-Nicolson discretization of the heat equation in [17], based on a piecewise quadratic polynomial in time reconstruction of the numerical solution. Finally, in a departure from [15], we replace the second derivatives in space (Laplacian of the discrete solution) in the error estimate with the forth derivatives in time by reusing the governing equation. This leads to the new a posteriori error estimate in time and also allows us to easily recover the error estimates in space that turn out to be the same as those of [15]. The resulting estimate is referred to as the 5-point estimator since it contains the fourth order finite differences in time and thus involves the discrete solution at 5 points in time at each time step. On the other hand, the estimate [15] involves only 3 points in time at each time step and will be thus referred to as the 3-point estimator.
Like in the case of the 3-point estimator, we are able to prove that the new 5-point estimator is reliable on general regular meshes in space and non-uniform meshes in time (with constants depending on the regularity of meshes in both space and time). Moreover, the 5-point estimator is proved to be of optimal order at least on sufficiently smooth solutions, quasi-uniform meshes in space and uniform meshes in time, again reproducing the results known for the 3-point estimator. Numerical experiments demonstrate that 3-point and 5-point error estimators produce very similar results in the majority of test cases. Both turn out to be of optimal order in space and time, even in situations not accessible to the current theory (non quasi-uniform meshes, not constant time steps). It should be therefore possible to use the new estimator for mesh adaptation in space and time. In fact, the best strategy in practice may be to combine both estimators to take benefit from the strengths of each of them: the relative cheapness of the 5-point one, and the better numerical behavior of the 3-point estimator under abrupt changes of the mesh.
The outline of the paper is as follows. We present the governing equations and the discretization in Section 2. Since our work is based on techniques from [15], Section 3 recalls the a posteriori bounds in time and space from there. In Section 4, the 5-point a posteriori error estimator for the fully discrete wave problem is derived. Numerical experiments on several test cases are presented in Section 5.

The Newmark scheme for the wave equation
We consider initial-boundary-value problem for the wave equation. Let Ω be a bounded domain in R 2 with boundary ∂Ω and T > 0 be a given final time. Let u = u(x, t) : Ω × [0, T ] → R be the solution to where f, u 0 , v 0 are given functions. Note that if we introduce the auxiliary unknown v = ∂u ∂t then model (2.1) can be rewritten as the following first-order in time system The above problem (2.1) has the following weak formulation [11]: for given f ∈ L 2 (0, T ; L 2 (Ω)), u 0 ∈ H 1 0 (Ω) and v 0 ∈ L 2 (Ω) find a function where ·, · denotes the duality pairing between H −1 (Ω) and H 1 0 (Ω) and the parentheses (·, ·) stand for the inner product in L 2 (Ω). Following Chapter 7, Section 2, Theorem 5 of [11], we observe that in fact Higher regularity results with more regular data are also available in [11].
Let us now discretize (2.1) or, equivalently, (2.2) in space using the finite element method and in time using an appropriate marching scheme. We thus introduce a regular mesh T h on Ω with triangles where E h represents the internal edges of the mesh T h and the standard finite element space V h ⊂ H 1 0 (Ω) of piecewise polynomials of degree k ≥ 1: Let us also introduce a subdivision of the time interval [0, T ] 0 = t 0 < t 1 < · · · < t N = T, with non-uniform time steps τ k = t n+1 − t k for n = 0, . . . , N − 1 and τ = max The Newmark scheme [18,19] with coefficients β = 1/4, γ = 1/2 as applied to the wave equation (2.1): given and then compute u n+1 h ∈ V h for n = 1, . . . , N − 1 from equation where f n is an abbreviation for f (·, t k ). Following [5,15], we observe that this scheme is equivalent to the Crank-Nicolson discretization of the governing equation written in the form (2.2): Note that the additional unknowns v h k are the approximations are not present in the Newmark scheme (2.5)-(2.6). If needed, they can be recovered on each time step by the following easy computation From now on, we shall use the following notations We apply this notations to all quantities indexed by a superscript, so that, for example, f n+1/2 = (f n+1 + f n )/2. We also denote u(x, t k ), v(x, t k ) by u n , v n so that, for example, u n+1/2 = u n+1 + u n /2 = (u(x, t n+1 ) + u(x, t k )) /2.
We shall measure the error in the following norm Here and in what follows, we use the notations u(t) and ∂u ∂t (t) as a shorthand for, respectively, u(·, t) and ∂u ∂t (·, t). The norms and semi-norms in Sobolev spaces H k (Ω) are denoted, respectively, by · H k (Ω) and |·| H k (Ω) .
We call (2.11) the energy norm referring to the underlying physics of the studied phenomenon. Indeed, the first term in (2.11) may be assimilated to the kinetic energy and the second one to the potential energy.

The 3-point time error estimator
The aim of this section is to recall a posteriori bounds in time and space from [15] for the error measured in the norm (2.11). Their derivation is based on the following piecewise quadratic (in time) 3-point reconstruction of the discrete solution.
Definition 3.1. Let u n h be the discrete solution given by the scheme (2.6). Then, the piecewise quadratic reconstructionũ hτ (t) : [0, T ] → V h is constructed as the continuous in time function that is equal on [t k , t n+1 ], n ≥ 1, to the quadratic polynomial in t that coincides with u n+1 h (respectively u n h , u n−1 h ) at time t n+1 (respectively t k , t n−1 ). Moreover,ũ hτ (t) is defined on [t 0 , t 1 ] as the quadratic polynomial in t that coincides with u 2 . The quadratic reconstructionsũ hτ ,ṽ hτ are thus based on three points in time (normally looking backwards in time, with the exemption of the initial time slab [t 0 , t 1 ]). This is also the case for the time error estimator (3.3), recalled in the following Theorem and therefore referred to as the 3-point estimator.
Theorem 3.2. The following a posteriori error estimate holds between the solution u of the wave equation (2.1) and the discrete solution u n h given by (2.5) and (2.6) for all t k , 0 ≤ n ≤ N with v n h given by (2.9): where the space indicator is defined by here C 1 , C 2 are constants depending only on the mesh regularity, [·] stands for a jump on an edge E ∈ E h , and u hτ ,ṽ hτ are given by Definition 3.1. The error indicator in time for k = 1, . . . , N − 1 is

4)
and η T (t 0 ) = 5 12 We also recall an optimality result for the 3-point time error estimator. We introduce to this end the , ). Suppose that mesh T h is quasi-uniform, the mesh in time is uniform (t k = kτ ), and the initial approximations are chosen as Then, the 3-point time error estimator η T (t k ) defined by (3.3, 3.5) is of order τ 2 , i.e.
with a positive constant C depending only on u, f , and the regularity of mesh T h .
Remark 3.4. Note that the particular choice for the approximation of initial conditions in (3.7) using the H 1 0 -orthogonal projection (3.6) is crucial to obtain the optimal order of the 3-point time error estimator, as confirmed both theoretically and numerically in [15].

The 5-point A P OSTERIORI error estimator
As already mentioned in the Introduction, the time error estimator (3.3) contains a finite element approximation to the Laplacian of u k h , i.e. z k h given by (3.4). This is unfortunate because z k h should be computed by solving an additional finite element problem that implies additional computational effort. Having in mind that the term ∂ 2 n f h − z n h in (3.3) is a discretization of ∂ 2 f /∂t 2 + ∆u = ∂ 4 u/∂t 4 at time t n our goal now is to avoid the second derivatives in space in the error estimates and replace them with the forth derivatives in time.
We introduce a "fourth order finite difference in time" ∂ 4 n defined by on any sequence {w n h } n=0,1,... ∈ V h . This can be rewritten as a composition of two second order finite difference operators where ∂ 2 w h is the standard finite difference (2.10) applied to w h , and∂ 2 n is a modified second order finite difference defined by∂ Note that a lower subscript "n" is lacking from ∂ 2 w h in (4.2) consistent with the fact that∂ 2 n is applied there to the sequence {∂ 2 n w h } n=0,1,... rather than to a single instance of ∂ 2 n w h . In full detail, (4.2) should be interpreted as Remark 4.1. In the case of constant time steps τ n = τ , (4.1) is reduced to It is thus indeed a standard finite difference approximation to the fourth derivative. In particular, it is exact on polynomials (in time) of degree up to 4. However, a standard fourth order finite difference in the general case of non constant time steps would be given by the divided differences Clearly, the formulas for ∂ 4 n w h and∂ 4 n w h , although similar, do not coincide in general, and consequently ∂ 4 n w h is not necessarily consistent with the fourth derivative in time of w h . Definition (4.1) may seem thus artificial and counter-intuitive. We shall see however that it arises naturally in the analysis of Newmark scheme, cf. forthcoming Lemma 4.2. Indeed, in order to "differentiate" in time the averaged quantitiesw n h defined by (4.4) and present in the scheme (2.6), cf. also (4.13), one needs to employ the modified second order finite differencê ∂ 2 n , which shall be composed further with ∂ 2 n to give rise to ∂ 4 n .
For any sequence {w n h } n=0,1,... ∈ V h , we denotē Consistently with the conventions above,w h will stand for the collection of any sequence {w n h } n=0,1,... . The following technical lemma establishes a connection between second order discrete derivatives∂ 2 n and ∂ 2 n .
where c and C are positive constants depending only on the mesh regularity in time, i.e. on max k≥0 Proof. We first note that relation (4.5) does not contain any derivatives in space and thus it should hold at any point x ∈ Ω. Consequently, it is sufficient to prove this Lemma assuming that w n h , ∂ 2 k w h , etc. are real numbers, i.e. replacing V h by R. This is the assumption adopted in this proof. We shall thus drop the sub-indexes h everywhere. Furthermore, it will be convenient to reinterpret w n in (4.2), (4.3) and (4.4) as the values of a real valued function w(t) at t = t n . We shall also use the notations likew n , ∂ 2 n w, and so on, where w is a continuous function on R, always assuming w n = w(t n ).
Observe that∂ 2 nw is a linear combination of 5 numbers {w n−3 , . . . , w n+1 }. Thus, it is enough to check equality (4.5) on any 5 continuous functions φ (k) (t), k = n − 3, . . . , n + 1, such that the vector of values of φ (k) at times t l , l = n − 3, . . . , n + 1, form a basis of R 5 . For fixed n, let us choose these functions as First we notice that for every linear function u(t) on [t n−3 , t n+1 ] we have∂ 2 nū = ∂ 2 n u = 0. Thus, we get immediately∂ 2 nφ(n−3) = ∂ 2 n φ (n−3) = 0 and∂ 2 kφ (n+1) = ∂ 2 k φ (n+1) = 0 so that (4.5) is fulfilled on functions φ (n−3) , φ (n+1) with any coefficients α k , k = n − 2, n − 1, n. Now we want to provide coefficients α k , k = n − 2, n − 1, n for which (4.5) is fulfilled on functions φ (n−2) , φ (n−1) and φ (n) . For brevity, we demonstrate the idea only for function φ (n) (t). Function φ (n) (t) is linear on [t n−3 , t n ] and thus From direct computations it is easy to show that where ∼ hides some factors that can be bounded by constants depending only on the mesh regularity. Thus we are able to establish expression for coefficient α n =∂ 2 kφ (n) ∂ 2 n φ (n) ≤ C. Similar reasoning for function φ (n−1) and φ (n−2) shows that α n−1 =∂ 2 nφ(n−1) The next step is to show boundedness from below of n k=n−2 α k . We will show it by applying equality (4.5) to second order polynomial function s(t) = t 2 2 . Using a Taylor expansion of s(t) aroundt n in the definition ofs n givess n = τ n (t 2 n +t n τ n−1 + 1 4 (τ 2 n + τ 2 n−1 )) + τ n−1 (t 2 n −t n τ n + 1 4 (τ 2 n + τ 2 n−1 )) 2(τ n + τ n−1 ) Substituting this into the definition of∂ 2 ns we obtain Using (4.5) and the fact that ∂ 2 n s = 1 for k = n − 2, n − 1, n we note that This implies n k=n−2 α k ≥ C. For all n ≥ 3 there exist coefficients β k , k = n − 2, n − 1, n such that n k=n−2 where coefficients α k , k = n − 2, n − 1, n are introduced in Lemma 4.2. Moreover where C is a positive constant depending only on the mesh regularity in time, i.e. on max k≥0 Proof. As in proof of Lemma 4.2, we assume V h = R, drop the sub-indexes h and interpret w n , s n as the values of continuous real valued functions w(t), s(t) at t = t n . Using (4.7) and notations (2.10) implies ∂ 2 k w = ∂ k s. Now, we are able to rewrite (4.8) in terms of s n only n k=n−2 As in the proof of Lemma 4.2 we take into account the fact that equation (4.9) should hold for every 5 numbers {s n−3 , . . . , s n+1 } and therefore it's enough to check equality (4.9) on 5 linearly independent piecewise linear functions φ (k) introduced by (4.6). Using the reasoning as in Lemma 4.2 leads to desired result (4.8).
We can now prove an a posteriori error estimate involving ∂ 4 n u h . Since the latter is computed through 5 points in time {t n−3 , . . . , t n+1 }, we shall refer to this approach as the 5-point estimator. For the same reason, this estimator is only applicable from time t 4 . The error at first 3 time steps should be thus measured differently, for example using the 3-point estimator from Theorem 3.2.
Theorem 4.4. The following a posteriori error estimate holds between the solution u of the wave equation (2.1) and the discrete solution u n h given by (2.5)-(2.6) for all t n , 4 ≤ n ≤ N with v n h given by (2.9): where the space error indicator is defined by (3.2) and the time error indicator iŝ with additional higher order termsη The constant C > 0 depends only on the mesh regularity in time, i.e. on max k≥0 Proof. We note first of all that it is sufficient to prove the Theorem for the final time, i.e. n = N because the statement for the general case n < N will follow by resetting the final time N to n. Introducing the L 2 -orthogonal projection P h : we can rewrite scheme (2.6) as for n = 0, . . . , N − 1 wheref n h is defined through averaging (4.4) from f n h = P h f (t n , ·). Taking a linear combination of instances of (4.13) at steps n, n − 1, n − 2 with appropriate coefficients gives (4.14) Using the definition of operator∂ 2 n and re-introducing v n h by (2.7) leads tô with coefficients α k , β k introduced in Lemmas 4.2 and 4.3. Moreover, by Lemma 4.2 γ = n k=n−2 α k −1 is positive and bounded so that with γ k = γβ k that are all uniformly bounded on regular meshes in time. Similarly, Thus, The rest of the proof follows closely that of Theorem 3.2, cf. [15]. We adopt the vector notation where v = ∂u/∂t. Note that the first equation in (2.2) implies that Similarly, Newmark scheme (2.7) and (2.8) can be rewritten as The a posteriori analysis relies on an appropriate residual equation for the quadratic reconstructionŨ hτ = ũ hτ v hτ . We have thus for t ∈ [t n , t n+1 ], n = 1, . . . , N − 1 17) so that, after some simplifications, Consider now (4.16) at time steps n and n − 1. Subtracting one from another and dividing by τ n−1/2 yields Introduce the error between reconstructionŨ hτ and solution U to problem (4.15): (4.20) or, component-wise Taking the difference between (4.19) and (4.15) we obtain the residual differential equation for the error valid for t ∈ [t n , t n+1 ], n = 1, . . . , andĨ h : H 1 0 (Ω) → V h is a Clément-type interpolation operator which is also a projection [10,20]. Noting that (A∇E, ∇E) = 0 and Integrating (4.21) in time from t 3 to some t * ≥ t 3 yields and assume that t * ∈ [t 3 , t N ] is the point in time where Z attains its maximum and t * ∈ (t n , t n+1 ] for some n. We have for the first and second terms in (4.22) This follows from an integration by parts with respect to time and the estimates on operators Π h andĨ h , cf. [15], and gives rise to the space part of the error estimate (4.10). Indeed, we can summarize the bounds above as The third term in (4.22) is responsible for the time estimator. It can be written as Recalling that Z(t * ) is the maximum of Z(t) and using the estimate Ĩ h E v L 2 (Ω) ≤ C E v L 2 (Ω) we continue as Noting tm+1 tm |p m |dt ≤ 1 12 we can finally bound III as Summing together the estimates on the terms I, II, III, and recalling Z(t * ) ≥ Z(t N ) yields (4.10) at the final time t N .
Remark 4.5. The termsη h.o.t T (t k ) in (4.10) are (at least formally) of higher order thanη T (t k ). We propose therefore to ignoreη h.o.t T (t k ) in practice together with the integral of f −f τ , and to useη T (t k ) as the indicator of error due to the discretization in time. The following Theorem shows that the latter is indeed of optimal order τ 2 , at least for sufficiently smooth solutions, on quasi-uniform meshes in space and uniform meshes in time.
Proof. The result follows from Theorem 3.3 by using (4.14) and Lemma 4.2 . Remark 4.7. Note, that as in the case for 3-point error estimator, the approximation of initial conditions is crucial for the optimal rate of our time error estimator.

A toy model: a second order ordinary differential equation
Let us consider first the following ordinary differential equation with a constant A > 0. This problem serves as simplification of the wave equation in which we get rid of the space variable. The Newmark scheme reduces in this case to and the error becomes The 3-point and the 5-point a posteriori error estimates are then defined as follows: We define the following effectivity indices in order to measure the quality of the 3-point and the 5-point estimators We consider problem (5.1) with the exact solution u = cos( √ At), and the final time T = 1. The results of simulations with constant time steps τ n = τ = T /N are presented in Table 1. We observe that 3-point and 5-point estimators are divided by about 100 when the time step τ is divided by 10. The true error e(t N ) also behaves as O(τ 2 ) and hence both time error estimators behave as the true error.
In order to check the behavior of time error estimators for variable time step we take the previous example with the following time step ∀n : 0 ≤ n ≤ N τ n = 0.1τ * , if mod(n, 2) = 0, τ * , if mod(n, 2) = 1, (5.6)   where τ * is a given fixed value, see Table 2. As in the case of constant time step we have the equivalence between the true error and both estimated errors. We have plotted on Figure  τ kηT (t k ) compared to the error e(t n ). Table 3 contains the results for even more non-uniform time step ∀n : 0 ≤ n ≤ N τ n = 0.01τ * , if mod(n, 2) = 0, τ * , if mod(n, 2) = 1, (5.7) on otherwise the same test case. Note that in case when A = 100 and τ * = 0.001 the 5-points error estimator significantly over-predicts the true error, while the 3-point estimator remains very close to it. This effect is consistent with Theorem 3.2. Indeed, the constants the 5-point error estimator may depend on the meshes regularity in time.
Our conclusion is that both 3-point and 5-point a posteriori error estimators are reliable for the toy model (5.1), although not asymptotically exact. The effectivity indices range from 1 to around 8 so that the the estimates can be rather pessimistic with respect to the true error. Moreover, Figure 1 suggests that the effectivity of the error estimator can deteriorate in the long time simulations.

The error estimator for the wave equation on Delaunay mesh
We now report numerical results for initial boundary-value problem for the wave equation (2.1) using piecewise linear finite elements in space) and Newmark scheme with non-uniform time steps and study the behavior of the 3-point time error estimator (3.3) and the 5-point time error estimator (4.11). All the computations are done with the help of FreeFEM++ [16]. In practice, we compute the space estimator (3.2) as follows:  with I h denoting the nodal interpolator to piecewise linear functions. The quality of our error estimators in space and time is determined by the effectivity indices Thus, u is a Gaussian function, whose center moves from point (0.3, 0.3) at t = 0 to point (0.7, 0.7) at t = 1. The transport velocity 0.8t(1, 1) T is peaking at t = 1. We choose non-uniform time step for n = 1, . . . , N − 1 as The initial conditions are computed with the orthogonal projections as in (3.7), cf. Remarks 3.4 and 4.7. Unstructured Delaunay meshes in space are used in all the experiments. Numerical results are reported in Table 4. Note that this case is chosen so that the non-uniform time step is required, see Figure 2.
Referring to Table 4, we observe that when setting initial time step as τ 2 ∼ O(h) the error is divided by 2 each time h is divided by 2, consistent with e ∼ O(τ 2 + h). The space error estimator and the two time error estimators behave similarly and thus provide a good representation of the true error. Both effectivity indices tend to a constant value.
We therefore conclude that our space and time error estimators are reliable in the regime of non-uniform time steps and Delaunay space meshes. They separate well the two sources of the error and can be thus used for the mesh adaptation in space and time. In particularly, 3-point and 5-point time estimators become more and more close to each other when h and τ tend to 0.
Finally, we report in Table 5 the computational times of our simulations using either 3-point or 5-point error estimators for the Newmark scheme on uniform meshes in time (with time step τ ) and unstructured quasi-uniform Delaunay meshes in space with maximum mesh size h. We have used FreeFEM++ to implement the algorithm and run it on a modern laptop computer with Intel Core I7 processor and 16GB of memory. The reported CPU times correspond to the whole computation, including the construction of the mesh, setting up the initial conditions, and factorizing the matrices, which is done only once before entering into the time marching loop (we have used the default UMFPACK direct sparse solver). We observe that the advantage of the 5-point estimator over the 3-point one which grows with refining the mesh and is up to 20% in our experiments. Some preliminary results for a fully adaptive algorithm showing the behavior of the estimators from [15] and from this article in more realistic settings are available in the Ph.D. thesis of the first author [14].