# A posteriori error estimates
for leap-frog and cosine methods
for second order evolution problems

Emmanuil H. Georgoulis  Department of Mathematics, University of Leicester,
Leicester LE1 7RH, England UK, and School of Applied Mathematical and Physical
Sciences, National Technical University of Athens, Athens 15780, Greece.
<EMAIL_ADDRESS>Omar Lakkis  Department of Mathematics,
University of Sussex, Brighton BN1 9QH, England UK.
http://www.maths.sussex.ac.uk/Staff/OL Charalambos G. Makridakis  Department
of Mathematics, University of Sussex, Brighton BN1 9QH, England UK.
<EMAIL_ADDRESS>Juha M. Virtanen  Department of Mathematics,
University of Leicester, Leicester LE1 7RH, England UK<EMAIL_ADDRESS>
###### Abstract

We consider second order explicit and implicit two-step time-discrete schemes
for wave-type equations. We derive optimal order a posteriori estimates
controlling the time discretization error. Our analysis, has been motivated by
the need to provide a posteriori estimates for the popular _leap-frog_ method
(also known as _Verlet_ ’s method in molecular dynamics literature); it is
extended, however, to general cosine-type second order methods. The estimators
are based on a novel reconstruction of the time-dependent component of the
approximation. Numerical experiments confirm similarity of convergence rates
of the proposed estimators and of the theoretical convergence rate of the true
error.

##  1 Introduction

This work is concerned with second order explicit and implicit two-step time-
discrete schemes for wave-type equations. Our objective is to derive optimal
order a posteriori estimates controlling the time-discretization error. To the
best of our knowledge, error control for wave equations, discretized by
popular methods is limited so far to first order schemes [9, 13]. Despite the
importance of such wave-type problems, the lack of error control for time-
discretizations used extensively in applications is probably due to the two-
step character of these methods, and the associated technical issues. Our
analysis, has been motivated by the need to provide a posteriori estimates for
the _leap-frog_ method or, as often termed, _Verlet_ ’s method in the
molecular dynamics literature. It extends, however to general cosine-type
second order methods [5, 6].

Adaptivity and a posteriori error control for parabolic problems have been
developed in, e.g., [12, 25, 21, 15, 19, 8, 10, 17]. In particular, as far as
time discretization is concerned, all implicit one-step methods can be treated
within the framework developed in [2, 20, 3, 4, 18]. Although some of these
results apply (directly or after appropriate modifications) to the wave
equation also, when written as a first order system and discretised by
implicit Runge-Kutta or Galerkin schemes, this framework does not cover
popular two-step implicit or explicit time-discretisation methods. The recent
results in [9, 13] cover only first order time discrete schemes; see also [1]
for certain estimators to standard implicit time-stepping finite element
approximations of the wave equation. For earlier works on adaptivity for wave
equations from various perspectives we refer, e.g., to [16, 7, 23, 24].

#### Model problem and notation.

Let $(H,\langle{\cdot},{\cdot}\rangle)$ be a Hilbert space and
$\mathcal{A}:[0,T]\to D(\mathcal{A})$, positive definite, self-adjoint, linear
operator on $D(\mathcal{A})$, the domain of $\mathcal{A}$, which is assumed to
be dense in $H$. For time $t\in(0,T]$, we consider the linear second order
hyperbolic problem: find $u:[0,T]\to D(\mathcal{A})$, such that

$\displaystyle u^{\prime\prime}(t)+\mathcal{A}u(t)$
$\displaystyle=f(t)\quad\text{for }0<t\leq T,$ (1) $\displaystyle u(0)$
$\displaystyle=u_{0},$ $\displaystyle u^{\prime}(0)$ $\displaystyle=v_{0},$

where $f:[0,T]\to H$, $u_{0},v_{0}\in H$.

#### Leap-frog time-discrete schemes.

We shall be concerned with the popular leap-frog time-discrete scheme for (1).
We consider a subdivision of the time interval $(0,T]$ into disjoint
subintervals $(t^{n},t^{n+1}]$, $n=0,\dots,N-1$, with $t^{0}=0$ and $t^{N}=T$,
and we define $k_{n}:=t^{n+1}-t^{n}$, the time-step. For simplicity of the
presentation, we shall assume that $k_{n}=k$ is constant, although this is not
a restriction of the analysis below. Despite being two-step, the schemes
considered herein can be formulated for variable time steps also, with their
consistency and stability properties then being influenced accordingly, cf.
[22, 11]; the study of such extensions is out of the scope of this work. We
shall use the notation $t^{n+1/2}:=(t^{n+1}+t^{n})/2$.

The time-discrete leap-frog scheme (or Verlet’s method in the terminology of
initial value problems or of molecular dynamics) for the wave problem (1) is
defined by finding approximations $U^{n+1}\in D(\mathcal{A})$ of the exact
values $u^{n+1}:=u(t^{n+1})$, such that:

$\partial^{2}U^{n+1}+\mathcal{A}U^{n}=f^{n},\quad n=1,\dots,N-1,$ (2)

where $f^{n}:=f(t^{n})\in H$,

$\partial^{2}U^{n+1}:=\frac{\partial U^{n+1}-\partial
U^{n}}{k}=\frac{U^{n+1}-2U^{n}+U^{n-1}}{k^{2}},$ (3)

with

$\partial U^{n+1}:=\frac{U^{n+1}-U^{n}}{k},$

assuming knowledge of $U^{0}$ and $U^{1}$. We set $U^{0}:=u_{0}$ and we define
$U^{1}$ by

$\frac{\partial U^{1}-v_{0}}{k}+\frac{1}{2}\mathcal{A}U^{0}=\frac{1}{2}f^{0},$
(4)

where $f^{0}:=f(0)$ and $\partial U^{1}:=(U^{1}-U^{0})/k$. This is a widely
used and remarkable method in many ways: it is the only two-step explicit
scheme for second order problems which is second order accurate, it has
important conservation and geometric properties, as it is symplectic, and it
is very natural and simple to formulate and implement. We refer to the review
article [14] for a thorough discussion. Explicit schemes, such as (2) are
suited for the discretization of wave-type partial differential equations,
since their implementation requires a mild CFL-type condition of the form
$k/h\leq C$ (in contrast to parabolic problems,) where $h$ stands for the
space-discretization parameter.

#### Cosine methods.

The leap-frog scheme is a member of a general class of two-step methods for
second order evolution problems, which are based on the approximation of
cosine and are used extensively in practical computations. In a two-step
cosine method, for $n=1,\dots,N-1,$ we seek approximations $U^{n+1}$ such that

$\partial^{2}U^{n+1}+\big{[}\,q_{1}\,\mathcal{A}U^{n+1}-2p_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}=\big{[}\,q_{1}\,f^{n+1}-2p_{1}f^{n}+q_{1}\,f^{n-1}\,\big{]},\
$ (5)

where we assume that $p_{1}=q_{1}-\frac{1}{2}$ for second order accuracy; we
refer to [5, 6] for a detailed discussion and analysis of general multi-step
cosine schemes. In this case, the rational function
$r(x)=(1+p_{1}x^{2})/(1+q_{1}x^{2})$ is a second order approximation to the
cosine, in the sense that for $|x|$ sufficiently small,

$|r(x)-\cos(x)|\leq C\,x^{4}\,.$ (6)

When $q_{1}=0$ the above methods are explicit, and the condition
$p_{1}=q_{1}-\frac{1}{2}$ implies that the only explicit second order member
of this family is the leap-frog method (2).

In this work, we derive a posteriori error bounds in the
$\operatorname{L}^{\infty}($)-in-time-energy-in-space norm of the error. The
derived bounds are of optimal order, i.e., of the same order as the error
(which is known to be second order) for the class of the schemes considered
[5, 6]. This is verified by the numerical experiments presented herein. Our
approach is based on the following ingredients: first, we rewrite the schemes
as one-step system on staggered time grids. In turn, this can be seen as a
second order perturbation of the staggered midpoint method. Further, by
introducing appropriate interpolants, we arrive to a form which can be viewed
as perturbation of (1) written as a first-order system. Finally, we employ an
adaptation of the time reconstruction from [2], yielding the desired a
posteriori error estimates. An interesting observation is that our estimates
hold _without_ any additional time-step assumption, which at the fully
discrete level would correspond to a CFL-type restriction. Thus, in a
posteriori analysis, standard stability considerations of time discretisation
schemes might influence the behaviour of the estimator, but are not explicitly
required; the possible instability is sufficiently reflected by the behaviour
of the estimator, see Section 4. Although not done here, by employing space
reconstruction techniques it would be possible to derive error estimates for
fully discrete schemes in various norms, using ideas from [19, 17, 13].

The remaining of this work is organized as follows. In §2 we reformulate the
numerical methods appropriately – this is a crucial step in our approach. We
start with the leap-frog method and we continue with providing two alternative
reformulations of general cosine methods. In §3 we introduce appropriate time
reconstructions and we derive the error bounds. In §4 we present detailed
numerical experiments which yield experimental orders of convergence for the
estimators that are the same with those of the actual error. Finally, in §5 we
draw some concluding remarks.

##  2 Reformulation of the methods

It will be useful for the analysis to reformulate the methods as a system in
two staggered grids.

###  2.1 Leap-frog

Starting with the leap-frog method, we introduce the auxiliary variable

$V^{n+1/2}:=\partial U^{n+1},$ (7)

for $n=0,1,\dots,N-1$, and we set $V^{-1/2}:=2v_{0}-V^{1/2}$. (Note that,
then, $v_{0}=(V^{-1/2}+V^{1/2})/2$.) Also, we define $U^{-1}:=U^{0}-kV^{-1/2}$
and we observe that we have

$v_{0}=\frac{U^{1}-U^{-1}}{2k}.$

Further, we introduce the notation

$\partial V^{n+1/2}:=\frac{V^{n+1/2}-V^{n-1/2}}{k},\quad n=0,1,\dots,N-1,$ (8)

noting that the identity $\partial V^{1/2}=2(\partial U^{1}-v_{0})/k$ also
holds.

We can now write the method (2) as a system in the staggered form considered
in [14]:

$\begin{split}&\partial U^{n+1}-V^{n+1/2}=0,\\\ \quad&\partial
V^{n+1/2}+\mathcal{A}U^{n}=f^{n}\,,\end{split}$ (9)

for $n=0,1,\dots,N-1$.

Next, our goal is to recast (9) using globally defined piecewise linear
functions. We define $U:[-k,T]\to D(\mathcal{A})$ to be the piecewise linear
interpolant of the sequence $\\{U^{n}\\}_{n=-1}^{N}$, at the points
$\\{t^{n}\\}_{n=-1}^{N}$, with $t^{-1}:=-k$. In addition, let
$V:[-k/2,t^{N-1/2}]\to D(\mathcal{A})$ be the piecewise linear interpolant of
the sequence $\\{V^{n+1/2}\\}_{n=-1}^{N-1}$, at the points
$\\{t^{n+1/2}\\}_{n=-1}^{N-1}.$ Using the notation

$\begin{split}U^{n+1/2}:=&\ U(t^{n+1/2}),\\\ \quad V^{n}:=&\ V(t^{n}),\qquad
n=0,\dots,N-1\,.\end{split}$ (10)

we, then, have

$U^{n+1/2}=\frac{1}{2}(U^{n+1}+U^{n}),\quad
V^{n}=\frac{1}{2}(V^{n+1/2}+V^{n-1/2}),$ (11)

for $n=0,1,\dots,N-1$.

Hence, in view of (11), (9) implies

$\displaystyle\partial U^{n+1}-\frac{1}{2}(V^{n+1}+V^{n})$
$\displaystyle=-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2}),$ (12)
$\displaystyle\partial V^{n+1/2}+\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})$
$\displaystyle=f^{n}+\frac{1}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1}),$

for $n=0,\dots,N-1$. Upon defining the piecewise constant residuals

$R_{U}(t)|_{(t^{n-1/2},t^{n+1/2}]}\equiv
R_{U}^{n}:=\frac{1}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1}),$
$R_{V}(t)|_{(t^{n},t^{n+1}]}\equiv
R_{V}^{n+1/2}:=-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2})\,,$

it is easy to check that, given that the leap-frog method is second order (in
both $U^{n}$ and $V^{n+1/2}$), we have $R_{U}^{n}=O(k^{2})$ and
$R_{V}^{n+1/2}=O(k^{2}).$ Hence, (12) can be viewed as a second order
perturbation of the staggered mid-point method for (1) written as first order
system

$\displaystyle u^{\prime}-v=0\,,$ (13) $\displaystyle
v^{\prime}+\mathcal{A}u=f\,.$

In what follows, it will be useful to rewrite (12) as a perturbation of the
continuous system (13). To this end, we introduce two time interpolants onto
piecewise linear functions defined on the staggered grids: we define
$U_{1}:[0,T]\to D(\mathcal{A})$ to be the piecewise linear interpolant of the
sequence $\\{U^{n+1/2}\\}_{n=-1}^{N-1}$ and $V_{1}:[0,t^{N-1}]\to
D(\mathcal{A})$ to be the piecewise linear interpolant of the sequence
$\\{V^{n}\\}_{n=0}^{N-1}$. Then, (12) can be written as

$\displaystyle U^{\prime}-I_{0}V_{1}$ $\displaystyle=R_{V},$ (14)
$\displaystyle V^{\prime}+\mathcal{A}\tilde{I}_{0}U_{1}$
$\displaystyle=\tilde{I}_{0}f+R_{U},$

where we define the interpolators

$\begin{split}\tilde{I}_{0}:&\text{ piecewise constant midpoint interpolator
on }\\{(t^{n-1/2},t^{n+1/2}]\\}_{n=0}^{N-1},\\\ I_{0}&\text{ piecewise
constant midpoint interpolator on
}\\{(t^{n-1},t^{n}]\\}_{n=1}^{N-1}.\end{split}$ (15)

This formulation will be the starting point of our analysis in the next
section.

###  2.2 Cosine methods: Formulation 1

We shall see that cosine methods (5) can be reformulated in a similar way as a
staggered system. As in the leap-frog case we introduce the auxiliary variable

$V^{n+1/2}:=\partial U^{n+1},$ (16)

and we let

$\partial V^{n+1/2}:=\frac{V^{n+1/2}-V^{n-1/2}}{k},\quad n=0,1,\dots,N-1,$
(17)

Then the methods (5) can be rewritten in system form:

$\begin{split}&\partial U^{n+1}-V^{n+1/2}=\ 0,\\\ \quad&\partial
V^{n+1/2}+\big{[}\,q_{1}\,\mathcal{A}U^{n+1}-2p_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}\\\
=&\big{[}\,q_{1}\,f^{n+1}-2p_{1}f^{n}+q_{1}\,f^{n-1}\,\big{]}\,,\end{split}$
(18)

for $n=0,1,\dots,N-1$. Using the same notation and conventions as in the leap-
frog case, we observe, respectively,

$\begin{split}\big{[}\,q_{1}\,&\mathcal{A}U^{n+1}-2p_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}=\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})\\\
&-\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})+\big{[}\,q_{1}\,\mathcal{A}U^{n+1}-2p_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}\\\
=&\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})\\\
&-\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})+\mathcal{A}U^{n}+\big{[}\,q_{1}\,\mathcal{A}U^{n+1}-2q_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}\\\
=&\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})\\\
&-\frac{1}{4}\big{[}\mathcal{A}U^{n+1}-2\mathcal{A}U^{n}+\,\mathcal{A}U^{n-1}\,\big{]}+\big{[}\,q_{1}\,\mathcal{A}U^{n+1}-2q_{1}\mathcal{A}U^{n}+q_{1}\,\mathcal{A}U^{n-1}\,\big{]}\\\
=&\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})-\frac{(1-4q_{1})}{4}\big{[}\mathcal{A}U^{n+1}-2\mathcal{A}U^{n}+\,\mathcal{A}U^{n-1}\,\big{]}\,,\end{split}$
(19)

where we used the fact that $p_{1}=q_{1}-\frac{1}{2}.$ Therefore, as before we
conclude that

$\displaystyle\partial U^{n+1}-\frac{1}{2}(V^{n+1}+V^{n})$
$\displaystyle=-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2}),$ (20)
$\displaystyle\partial V^{n+1/2}+\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})$
$\displaystyle=\tilde{f}^{n}+\frac{(1-4q_{1})}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1}),$

where
$\tilde{f}^{n}=\big{[}\,q_{1}\,f^{n+1}-2p_{1}f^{n}+q_{1}\,f^{n-1}\,\big{]}$
and $n=0,\dots,N-1$. Let us now define

$\displaystyle R^{\cos}_{U}(t)|_{(t^{n-1/2},t^{n+1/2}]}\equiv
R_{U}^{\cos,n}:=$ $\displaystyle\
\frac{(1-4q_{1})}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1})$
$\displaystyle+q_{1}\,\big{[}\,\,f^{n+1}-2f^{n}+\,f^{n-1}\,\big{]},$
$R^{\cos}_{V}(t)|_{(t^{n},t^{n+1}]}\equiv
R_{V}^{\cos,n+1/2}:=-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2})\,.$

As in the leap-frog case, it is easy to check that given that the method is
second order, we have $R_{U}^{\cos,n}=O(k^{2})$ and
$R_{V}^{\cos,n+1/2}=O(k^{2}).$ Hence, (12) can be seen as a second order
perturbation of the staggered mid-point method for (13). Further, still using
the same notation for time interpolants as the leap-frog case, we obtain

$\displaystyle U^{\prime}-I_{0}V_{1}$ $\displaystyle=R^{\cos}_{V},$ (21)
$\displaystyle V^{\prime}+\mathcal{A}\tilde{I}_{0}U_{1}$
$\displaystyle={\tilde{I}_{0}\tilde{f}}+R^{\cos}_{U}\,.$

where $\tilde{I}_{0}\tilde{f}|_{(t^{n-1/2},t^{n+1/2}]}=\tilde{f}^{n}.$ It is
interesting to compare (21) to (14).

###  2.3 Cosine methods: Formulation 2

We briefly discuss an alternative formulation of cosine methods. This time we
let

$V^{n+1/2}:=\big{(}\,I+k^{2}q_{1}\mathcal{A}\,\big{)}\partial U^{n+1},$ (22)

and, as before,

$\partial V^{n+1/2}:=\frac{V^{n+1/2}-V^{n-1/2}}{k},\quad n=0,1,\dots,N-1.$
(23)

Then, using again $p_{1}=q_{1}-\frac{1}{2},$ we rewrite the methods (5) as

$\begin{split}\partial U^{n+1}-V^{n+1/2}=&-k^{2}q_{1}\mathcal{A}\,\partial
U^{n+1},\\\ \quad\partial
V^{n+1/2}+\mathcal{A}U^{n}=&\big{[}\,q_{1}\,f^{n+1}-2p_{1}f^{n}+q_{1}\,f^{n-1}\,\big{]}\,,\end{split}$
(24)

for $n=0,1,\dots,N-1$. Using the same notation and conventions as in the leap-
frog case, we finally conclude

$\displaystyle\partial
U^{n+1}-\frac{1}{2}(V^{n+1}+V^{n})=-k^{2}q_{1}\mathcal{A}\,\partial
U^{n+1}-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2}),$ (25)
$\displaystyle\partial
V^{n+1/2}+\frac{1}{2}\mathcal{A}(U^{n+1/2}+U^{n-1/2})=\tilde{f}^{n}+\frac{1}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1}),$

where
$\tilde{f}^{n}=\big{[}\,q_{1}\,f^{n+1}-2p_{1}f^{n}+q_{1}\,f^{n-1}\,\big{]}$
and $n=0,\dots,N-1$. Upon defining the perturbations as

$R^{\cos,2}_{U}(t)|_{(t^{n-1/2},t^{n+1/2}]}\equiv
R_{U}^{\cos,2,n}:=\frac{1}{4}\mathcal{A}(U^{n+1}-2U^{n}+U^{n-1})+q_{1}\,\big{[}\,\,f^{n+1}-2f^{n}+\,f^{n-1}\,\big{]},$
$R^{\cos,2}_{V}(t)|_{(t^{n},t^{n+1}]}\equiv
R_{V}^{\cos,2,n+1/2}:=-k^{2}q_{1}\mathcal{A}\,\partial
U^{n+1}-\frac{1}{4}(V^{n+3/2}-2V^{n+1/2}+V^{n-1/2})\,,$

it is easy to check, again, that the method is second order
($R^{\cos,2,n}_{U}=O(k^{2})$ and $R_{V}^{\cos,2,n+1/2}=O(k^{2}),$) and, thus,
(12) can be interpreted as a second order perturbation of the staggered mid-
point method for (13). Still, using the same notation as before, we have

$\displaystyle U^{\prime}-I_{0}V_{1}$ $\displaystyle=R^{\cos,2}_{V},$ (26)
$\displaystyle V^{\prime}+\mathcal{A}\tilde{I}_{0}U_{1}$
$\displaystyle=\tilde{I}_{0}\tilde{f}+R^{\cos,2}_{U}\,.$

##  3 A posteriori error bounds

We have seen that all above schemes can be written in the form,

$\displaystyle V^{\prime}+\mathcal{A}\tilde{I}_{0}U_{1}$
$\displaystyle=\tilde{I}_{0}(f+\rho_{U}),$ (27) $\displaystyle
U^{\prime}-I_{0}V_{1}$ $\displaystyle=I_{0}\rho_{V},$

where $\tilde{I}_{0},I_{0}$ are defined in (15) and $\tilde{I}_{0}\rho_{U}$
equal to $R_{U},R^{\cos}_{U}$, or $R^{\cos,2}_{U}$, for the leap-frog, the
first and second cosine method formulations, respectively; similarly,
$I_{0}\rho_{V}$ is equal to $R_{V},R^{\cos}_{V}$, or $R^{\cos,2}_{V}$ for each
of the respective 3 formulations; cf., (14), (21) and (26). It is possible, in
principle, to consider non-constant $\rho_{U},\rho_{V}$ on each time-step;
nevertheless the, easiest to implement, constant ones considered here suffice
to deliver optimal estimator convergence rates, as will be highlighted in the
numerical experiments below.

###  3.1 Reconstructions

We continue by defining appropriate time reconstructions, cf., [2]. To this
end, on each interval $(t^{n-1/2},t^{n+1/2}]$, for $n=0,\dots,N-1$, we define
the reconstruction $\hat{V}$ of $V$ by

$\hat{V}(t):=V^{n-1/2}+\int_{t_{n-1/2}}^{t}(-\mathcal{A}U_{1}+\tilde{I}_{1}f+\rho_{U})\,\mathrm{d}t,$

where $\tilde{I}_{1}$ is a piecewise linear interpolant on the mesh
$\\{(t^{n-1/2},t^{n+1/2}]\\}_{n=1}^{N-1}$, such that
$\tilde{I}_{1}f(t^{n})=f^{n}$. We observe that $\hat{V}(t^{n-1/2})=V^{n-1/2}$,
and

$\hat{V}(t^{n+1/2})=V^{n-1/2}+k(-\mathcal{A}U_{1}(t^{n})+\tilde{I}_{1}f(t^{n})+\rho_{U}(t^{n}))=V^{n+1/2},$

using the mid-point rule to evaluate the integral and the first equation in
(12).

Also, on each interval $(t^{n-1},t^{n}]$, for $n=1,\dots,N-1$, we define the
reconstruction $\hat{U}$ of $U$ by

$\hat{U}(t):=U^{n-1}+\int_{t_{n-1}}^{t}(V_{1}+\rho_{V})\,\mathrm{d}t.$

Again, we observe that $\hat{U}(t^{n-1})=U^{n-1}$ and that

$\hat{U}(t^{n})=U^{n-1}+k(V_{1}(t^{n-1/2})+\rho_{V}(t^{n-1/2}))=U^{n},$

using the mid-point rule. Notice that each of the above reconstructions is
similar in spirit to the Crank-Nicolson reconstruction of [2]; however, we
note that, although $\hat{U},\hat{V}$ are both globally continuous functions,
their derivatives jump alternatingly at the nodes of the staggered grid.

###  3.2 Error equation and estimators

Setting $\hat{e}_{U}:=u-\hat{U}$ and $\hat{e}_{V}:=u^{\prime}-\hat{V}$, we
deduce

$\displaystyle\hat{e}_{V}^{\prime}+\mathcal{A}\hat{e}_{U}$
$\displaystyle=\mathcal{R}_{1}+\mathcal{R}_{f}$ (28)
$\displaystyle\hat{e}_{U}^{\prime}-\hat{e}_{V}$
$\displaystyle=\mathcal{R}_{2},$

with

$\displaystyle\mathcal{R}_{1}$
$\displaystyle:=-\mathcal{A}(\hat{U}-U_{1})-\rho_{U},$ (29)
$\displaystyle\mathcal{R}_{2}$ $\displaystyle:=\hat{V}-V_{1}-\rho_{V},$
$\displaystyle\mathcal{R}_{f}$ $\displaystyle:=f-\tilde{I}_{1}f.$

For $\Phi=(\phi_{1},\phi_{2})$, $\Psi=(\psi_{1},\psi_{2})$ $\in
D(\mathcal{A})\times H$, we define the bilinear form

$\langle\\!\langle{\Phi},{\Psi}\rangle\\!\rangle:=\langle{\mathcal{A}^{1/2}\phi_{1}},{\mathcal{A}^{1/2}\psi_{1}}\rangle+\langle{\phi_{2}},{\psi_{2}}\rangle.$

It is evident that $\langle\\!\langle{\cdot},{\cdot}\rangle\\!\rangle$ is an
inner product on $[D(\mathcal{A}^{1/2})\times H]^{2}$. This is the standard
energy inner product and the induced norm, denoted by
$|\kern-0.7113pt\|{\cdot}|\kern-0.7113pt\|,$ i.e.,

$|\kern-0.7113pt\|{\Phi}|\kern-0.7113pt\|=(||\mathcal{A}^{1/2}\phi_{1}||^{2}+||\phi_{2}||^{2})^{1/2}\,,$

is the natural energy norm for (1).

Then, the a posteriori error estimates will follow by applying standard energy
arguments to the error equation (28). More specifically, in view of (28), we
have

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})}|\kern-0.7113pt\|^{2}$
$\displaystyle=\langle\\!\langle{(\hat{e}_{U}^{\prime},\hat{e}_{V}^{\prime})},{(\hat{e}_{U},\hat{e}_{V})}\rangle\\!\rangle$
$\displaystyle=\langle{\mathcal{A}\hat{e}_{U}^{\prime}},{\hat{e}_{U}}\rangle+\langle{\hat{e}_{V}^{\prime}},{\hat{e}_{V}}\rangle$
$\displaystyle=\langle{\mathcal{A}\hat{e}_{V}},{\hat{e}_{U}}\rangle+\langle{\mathcal{A}\mathcal{R}_{2}},{\hat{e}_{U}}\rangle-\langle{\mathcal{A}\hat{e}_{U}},{\hat{e}_{V}}\rangle+\langle{\mathcal{R}_{1}},{\hat{e}_{V}}\rangle+\langle{\mathcal{R}_{f}},{\hat{e}_{V}}\rangle$
$\displaystyle=\langle{\mathcal{A}\mathcal{R}_{2}},{\hat{e}_{U}}\rangle+\langle{\mathcal{R}_{1}},{\hat{e}_{V}}\rangle+\langle{\mathcal{R}_{f}},{\hat{e}_{V}}\rangle,$

using the self-adjointness of $\mathcal{A}$. Hence, using the Cauchy-Schwarz
inequality, we arrive to

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})}|\kern-0.7113pt\|^{2}$
$\displaystyle\leq|\kern-0.7113pt\|{(\mathcal{R}_{2},\mathcal{R}_{1}+\mathcal{R}_{f})}|\kern-0.7113pt\||\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})}|\kern-0.7113pt\|.$

Integrating between $0$ and $\tau$, with $0\leq\tau\leq t^{N}$ such that

$|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(\tau)}|\kern-0.7113pt\|{}=\sup_{t\in[0,t^{N}]}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(t)}|\kern-0.7113pt\|,$

we arrive to

$\frac{1}{2}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(\tau)}|\kern-0.7113pt\|^{2}\leq\frac{1}{2}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(0)}|\kern-0.7113pt\|^{2}+|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(\tau)}|\kern-0.7113pt\|\int_{0}^{\tau}|\kern-0.7113pt\|{(\mathcal{R}_{2},\mathcal{R}_{1}+\mathcal{R}_{f})}|\kern-0.7113pt\|\mathrm{d}t,$

which implies

$|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(\tau)}|\kern-0.7113pt\|^{2}\leq
2|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(0)}|\kern-0.7113pt\|^{2}+4\Big{(}\int_{0}^{\tau}|\kern-0.7113pt\|{(\mathcal{R}_{2},\mathcal{R}_{1}+\mathcal{R}_{f})}|\kern-0.7113pt\|\mathrm{d}t\Big{)}^{2}.$

This already gives the following a posteriori bound.

######  Theorem 1.

Let $u$ be the solution of (1), $\hat{e}_{U}:=u-\hat{U}$ and
$\hat{e}_{V}:=u^{\prime}-\hat{V}.$ Then, the following a posteriori error
estimate holds

$\sup_{t\in[0,t^{N}]}|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(t)}|\kern-0.7113pt\|^{2}\leq
2|\kern-0.7113pt\|{(\hat{e}_{U},\hat{e}_{V})(0)}|\kern-0.7113pt\|^{2}+4\Big{(}\int_{0}^{t^{N}}|\kern-0.7113pt\|{(\mathcal{R}_{2},\mathcal{R}_{1}+\mathcal{R}_{f})}|\kern-0.7113pt\|\mathrm{d}t\Big{)}^{2},$

where $\mathcal{R}_{2},\mathcal{R}_{1}$, and $\mathcal{R}_{f}$ are defined in
(29).

An immediate Corollary from Theorem 1 is an posteriori bound for the error
$\sup_{[0,t^{N}]}|\kern-0.7113pt\|{(u-U,u^{\prime}-V)}|\kern-0.7113pt\|$ which
can be trivially deduced through a triangle inequality.

######  Remark 2.

Notice that due to the two-mesh stagerring, the computation of the ‘last’ $V$
used in the above estimate, $V^{N-1/2},$ requires the computation of
$U^{N+1}.$ This can be obtained by advancing one more time step in the
computation before estimating.

##  4 Numerical Experiments

###  4.1 Fully discrete formulation.

Although the focus of the present work is in time-discretization, we shall
introduce a fully discrete version of the time-stepping schemes for the
numerical experiments bellow. To this end, we consider Let
$\Omega\subset\mathbb{R}^{d}$, $d=2,3$ be a domain with boundary
$\partial\Omega$. We consider the initial-boundary value problem for the wave
equation: find
$u\in\operatorname{L}^{\infty}(0,T;\operatorname{H}^{1}_{0}(\Omega))$ such
that,

$\displaystyle u_{tt}+\mathcal{A}u$ $\displaystyle=$ $\displaystyle
f\quad\;\text{ in }\;\Omega\times(0,T],$ (30) $\displaystyle u$
$\displaystyle=$ $\displaystyle u_{0}\quad\text{ in
}\;\Omega\times\\{0\\}\quad$ (31) $\displaystyle u_{t}$ $\displaystyle=$
$\displaystyle v_{0}\quad\text{ in }\;\Omega\times\\{0\\}\quad\text{}$ (32)
$\displaystyle u$ $\displaystyle=$ $\displaystyle g\quad\text{ on
}\;\partial\Omega\times(0,T],$ (33)

where, for simplicity, we take $\mathcal{A}=-c^{2}\Delta,c\neq 0,$ and
$g\in\operatorname{H}^{1/2}(\partial\Omega)$. Further, for each $n$, we
consider the standard, conforming finite element space
$S_{h}^{p}\subset\operatorname{H}^{1}_{0}(\Omega)$, based on a quasiuniform
triangulation of $\Omega$ consisting of finite elements of polynomial degree
$p$, with $h$ denoting the largest element diameter. Focusing on time-
discretization issues, we shall use the same spatial discretization for all
time-steps. The respective discrete spatial operator is denoted by
$\mathcal{A}_{h}^{p}$. The fully-discrete leap-frog method is then defined as
follows: for each $n=2,\dots,N$, find $U^{n+1}\in S_{h}^{p}$, such that

$\displaystyle U^{n+1}$
$\displaystyle=2U^{n}-U^{n-1}+k^{2}(\bar{f}^{n}-\mathcal{A}_{h}^{p}U^{n}),$
(34)

and $U^{1}\in S_{h}^{p}$ such that

$\displaystyle U^{1}$
$\displaystyle=U^{0}+kV^{0}+\frac{k^{2}}{2}(\mathcal{A}_{h}^{p}U^{0}-\bar{f}^{0});$
(35)

here $\bar{f}^{n}(\cdot):=\Pi f(\cdot,t_{n})$ for each $t_{n}$, where
$\Pi:\operatorname{L}^{2}(\Omega)\to S_{h}^{p}$ denotes a suitable
interpolation/projection operator onto the finite element space $S_{h}^{p}$ of
the source function $f$. We also set $U^{0}:=\Pi u_{0}$ and $V^{0}:=\Pi
v_{0}$. Note that $V$ can be calculated as above through $U$, or can be
computed as follows: find $V^{n+1/2}\in S_{h}^{p}$ such that

$\displaystyle V^{n+1/2}$
$\displaystyle=V^{n-1/2}+k(\bar{f}^{n}-\mathcal{A}_{h}^{p}U^{n});$ (36)

(36) can be used to overcome the difficulty of evaluating estimators defined
on staggered time mesh and depending on the term $V^{n+3/2}$.

To assess the time-error estimator, we replace $\mathcal{A}$ by its
approximation $\mathcal{A}_{h}^{p}$ in the a posteriori estimators discussed
above, and in the $((\cdot,\cdot))$-inner product and $|||\cdot|||$-norm. For
brevity, we introduce the notation:

$e_{R}:=(\hat{e}_{U},\hat{e}_{V})\quad\text{and}\quad
e_{L}:=e_{R}+(U-\hat{U},V-\hat{V})=(u-U,u_{t}-V).$

The objective is to study the performance of the a posteriori estimator

$\eta_{1}:={\left(2|||e_{R}(0)|||^{2}+4{\left(\int_{0}^{T}|||(R_{2},R_{1}+R_{f})|||\right)}^{2}\right)}^{1/2},$
(37)

from Theorem 1.

###  4.2 Specific tests

For $t\in[0,1]$ and $\Omega:=(0,1)^{2}$, we consider the model problem (30) -
(32) in the following set up: $\mathcal{A}:=-c^{2}\Delta$ and $f=0$, for $c>0$
constant. Then, the exact solution $u$ of problem (30) - (32) is given by:

$u(t,x):=\sum_{j,k=1}^{\infty}\sin(k\pi x)\sin(j\pi
x)\left(\alpha_{k,j}\cos(\xi_{k,j}\pi t)+\beta_{k,j}\sin(\xi_{k,j}\pi
t)\right)$ (38)

where $\alpha_{k,j}>0$, $\beta_{k,j}>0$ and $\xi_{k,j}:=c\sqrt{k^{2}+j^{2}}$.
To illustrate the estimator’s behaviour, we have chosen the following sets of
parameters as numerical examples:

$\displaystyle\left\\{\begin{aligned} \text{ }&c:=1.0,\\\ \text{
}&\alpha_{1,1}=\beta_{1,1}=15.0,\beta_{k,j}=\alpha_{k,j}=0\;\text{for
}\;k,j\neq 1,\\\ \end{aligned}\right.$ (39)
$\displaystyle\left\\{\begin{aligned} \text{ }&c:=1.0,\\\ \text{
}&\alpha_{3,3}=\beta_{3,3}=1.0,\beta_{k,j}=\alpha_{k,j}=0\;\text{for
}\;k,j\neq 3,\\\ \end{aligned}\right.$ (40)
$\displaystyle\left\\{\begin{aligned} \text{ }&c:=5.0,\\\ \text{
}&\alpha_{1,1}=\beta_{1,1}=15.0,\beta_{k,j}=\alpha_{k,j}=0\;\text{for
}\;k,j\neq 1,\\\ \end{aligned}\right.$ (41)

The solutions of (39) - (41) are all smooth, but (41) oscillates much faster
temporally, while (40) has greater space-dependence of the error. In the
numerical experiments, the C++ library FEniCS/dolfin 1.2.0 and PETSc/SuperLU
were used for the finite element formulation and the linear algebra
implementation. For each of the examples, we compute the solution of (34)
using finite element spaces of polynomial degree $p=2$, and time step size
$k=Ch^{r}/(p+1)^{2}$, $r=1,2$, for some constant $C>0$, with $h>0$ denoting
the diameter of the largest element in the mesh associated with $S_{h}^{p}$.
The sequences of meshsizes considered (with respective colouring in the
figures below) are: $h=1/2$ (cyan),$1/4$ (green), $1/(4\sqrt{2})$ (yellow),
$1/8$ (red), $1/10$ (purple), for Examples (39) and (41), and
$h=1/(4\sqrt{2})$ (cyan), $1/8$ (green), $1/10$ (yellow), $1/12$ (red), $1/14$
(purple) for Example (40). Note that the CFL-condition required by the leap-
frog method, is satisfied by the restriction on the time step

$k\leq C\frac{h}{(p+1)^{2}},$ (42)

for sufficiently small $C>0$ constant, which couples the temporal and spatial
discretization sizes. We also include an experiment to highlight the behaviour
of the estimator when the CFL condition is violated.

We monitor the evolution of the values and the experimental order of
convergence of the estimator $\eta_{1}$, and the errors $e_{R}$ and $e_{L}$,
as well as the effectivity index over time on a sequence of uniformly refined
meshes with mesh sizes given as per each example. We also monitor the energy
of the reconstructed solution:

$E_{\text{reconstruction}}:=(\hat{U},\hat{V}).$ (43)

We define _experimental order of convergence_ ($\operatorname{EOC}$) of a
given sequence of positive quantities $a(i)$ defined on a sequence of meshes
of size $h(i)$ by

$\operatorname{EOC}(a,i)=\frac{\log(a(i+1)/a(i))}{\log(h(i+1)/h(i))},$ (44)

the _inverse effectivity index_ ,

$\operatorname{IEI}(\|{e}\|_{\operatorname{L}^{\infty}(0,t_{m};|||\cdot|||)},\eta_{1})=\frac{\|{e}\|_{\operatorname{L}^{\infty}(0,t_{m};|||\cdot|||)}}{\eta_{1}}.$
(45)

The $\operatorname{IEI}$has the same information as the (standard) effectivity
index and has the advantage of relating directly to the inequality appearing
in Theorem 1. The results of numerical experiments on uniform meshes, depicted
in Figures 1 and 2, indicate that the error estimators are reliable and also
efficient provided the time steps are kept sufficiently small. In the last
experiment, Figure 4(a), the behaviour of the estimator is displayed, for the
case when the CFL condition (42) is violated; the estimator remains reliable
in this case also: the instability is reflected in the behaviour of estimator,
cf., Figure 4(a). Indeed, the inverse effectivity idex $\operatorname{IEI}$ of
the estimator oscillates around a constant value for all times in the course
of the numerically unstable behaviour.

##  5 Concluding remarks

An a posteriori error bound for error measured in $\operatorname{L}^{\infty}$
\- norm in time and energy norm in space for leap-frog and cosine-type time
semi-discretizations for linear second order evolution problems was presented
and studied numerically. The estimator was found to be reliable, with the same
convergence rate as the theoretical convergence rate of the error. In a fully
discrete setting, this estimator corresponds to the control of time
discretization error. The estimators were also found to be sharp on uniform
meshes provided that the time steps and, thus, the time-dependent part of the
error, is kept sufficiently small. Investigation into the suitability of the
proposed estimators within an adaptive algorithm remains a future challenge.

## References

  * [1] Slimane Adjerid, A posteriori finite element error estimation for second-order hyperbolic problems, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 4699–4719.
  * [2] Georgios Akrivis, Charalambos Makridakis, and Ricardo H. Nochetto, A posteriori error estimates for the Crank-Nicolson method for parabolic equations, Math. Comp., 75 (2006), pp. 511–531 (electronic).
  * [3] , Optimal order a posteriori error estimates for a class of Runge-Kutta and Galerkin methods, Numer. Math., 114 (2009), pp. 133–160.
  * [4] , Galerkin and Runge-Kutta methods: unified formulation, a posteriori error estimates and nodal superconvergence, Numer. Math., 118 (2011), pp. 429–456.
  * [5] Garth A. Baker, Vassilios A. Dougalis, and Steven M. Serbin, High order accurate two-step approximations for hyperbolic equations, RAIRO Anal. Numér., 13 (1979), pp. 201–226.
  * [6] , An approximation theorem for second-order evolution equations, Numer. Math., 35 (1980), pp. 127–142.
  * [7] Wolfgang Bangerth and Rolf Rannacher, Adaptive finite element techniques for the acoustic wave equation, J. Comput. Acoust., 9 (2001), pp. 575–591.
  * [8] A. Bergam, C. Bernardi, and Z. Mghazli, A posteriori analysis of the finite element discretization of some parabolic equations, Math. Comp., 74 (2005), pp. 1117–1138 (electronic).
  * [9] Christine Bernardi and Endre Süli, Time and space adaptivity for the second-order wave equation, Math. Models Methods Appl. Sci., 15 (2005), pp. 199–225.
  * [10] Christine Bernardi and Rüdiger Verfürth, A posteriori error analysis of the fully discretized time-dependent Stokes equations, M2AN Math. Model. Numer. Anal., 38 (2004), pp. 437–455.
  * [11] M. P. Calvo and J. M. Sanz-Serna, The development of variable-step symplectic integrators, with application to the two-body problem, SIAM J. Sci. Comput., 14 (1993), pp. 936–952.
  * [12] Kenneth Eriksson and Claes Johnson, Adaptive finite element methods for parabolic problems. II. Optimal error estimates in ${L}_{\infty}{L}_{2}$ and ${L}_{\infty}{L}_{\infty}$, SIAM J. Numer. Anal., 32 (1995), pp. 706–740.
  * [13] Emmanuil H. Georgoulis, Omar Lakkis, and Charalambos Makridakis, A posteriori $L^{\infty}(L^{2})$-error bounds for finite element approximations to the wave equation, IMA J. Numer. Anal., 33 (2013), pp. 1245–1264.
  * [14] Ernst Hairer, Christian Lubich, and Gerhard Wanner, Geometric numerical integration illustrated by the Störmer-Verlet method, Acta Numer., 12 (2003), pp. 399–450.
  * [15] Paul Houston and Endre Süli, Adaptive Lagrange-Galerkin methods for unsteady convection-diffusion problems, Math. Comp., 70 (2001), pp. 77–106.
  * [16] Claes Johnson, Discontinuous Galerkin finite element methods for second order hyperbolic problems, Comput. Methods Appl. Mech. Engrg., 107 (1993), pp. 117–129.
  * [17] Omar Lakkis and Charalambos Makridakis, Elliptic reconstruction and a posteriori error estimates for fully discrete linear parabolic problems, Math. Comp., 75 (2006), pp. 1627–1658 (electronic).
  * [18] Christian Lubich and Charalambos Makridakis, Interior a posteriori error estimates for time discrete approximations of parabolic problems, Numer. Math., 124 (2013), pp. 541–557.
  * [19] Charalambos Makridakis and Ricardo H. Nochetto, Elliptic reconstruction and a posteriori error estimates for parabolic problems, SIAM J. Numer. Anal., 41 (2003), pp. 1585–1594 (electronic).
  * [20] , A posteriori error analysis for higher order dissipative methods for evolution problems, Numer. Math., 104 (2006), pp. 489–514.
  * [21] Marco Picasso, Adaptive finite elements for a linear parabolic problem, Comput. Methods Appl. Mech. Engrg., 167 (1998), pp. 223–237.
  * [22] Robert D. Skeel, Variable step size destabilizes the Störmer/leapfrog/Verlet method, BIT, 33 (1993), pp. 172–175.
  * [23] E. Süli, A posteriori error analysis and global error control for adaptive finite volume approximations of hyperbolic problems, in Numerical analysis 1995 (Dundee, 1995), vol. 344 of Pitman Res. Notes Math. Ser., Longman, Harlow, 1996, pp. 169–190.
  * [24] Endre Süli, A posteriori error analysis and adaptivity for finite element approximations of hyperbolic problems, in An introduction to recent developments in theory and numerics for conservation laws (Freiburg/Littenweiler, 1997), vol. 5 of Lect. Notes Comput. Sci. Eng., Springer, Berlin, 1999, pp. 123–194.
  * [25] R. Verfürth, A posteriori error estimates for finite element discretizations of the heat equation, Calcolo, 40 (2003), pp. 195–212.

Fig. 1:  Examples (39) and (40).

(a) Example (39). Errors, estimator and $\operatorname{IEI}$ are plotted on
the top row and $\operatorname{EOC}$s and energy of the reconstructed solution
on the bottom row over time ($x$-axis). Results are computed on the sequence
of uniform meshes with mesh size $h$, fixed time step $k=0.4h/(p+1)^{2}$ and
$p=2$. The $\operatorname{IEI}$ behaviour indicates that the error is well
estimated by the estimator and the convergence rate of the estimator remains
near to $\operatorname{EOC}\approx 2$, i.e., to that of the errors $e_{L}$ and
$e_{R}$.  (b) Example (40). Errors, estimator and $\operatorname{IEI}$ are
plotted on the top row and $\operatorname{EOC}$s and energy of the
reconstructed solution on the bottom row over time ($x$-axis). Results are
computed on the sequence of uniform meshes with mesh size $h$, fixed time step
$k=0.4h/(p+1)^{2}$ and $p=2$. The $\operatorname{IEI}$ behaviour indicates
that the error is well estimated by the estimator and the convergence rate of
the estimator remains near to $\operatorname{EOC}\approx 2$, i.e., to that of
the errors $e_{L}$ and $e_{R}$.

Fig. 2:  Example (41).

(a) Example (41). Errors, estimator and $\operatorname{IEI}$ are depicted on
the top row and $\operatorname{EOC}$s and energy of the reconstructed solution
on the bottom row over time ($x$-axis). Results are computed on the sequence
of uniform meshes with mesh size $h$, time-step $k=0.1h/(p+1)^{2}$ and $p=2$.
The $\operatorname{IEI}$ behaviour indicates that the error is overestimated
by the estimator. The convergence rate of the estimator is slightly below $2$,
probably due to the somewhat coarse time-step for asymptotic convergence, cf.,
subfigure (b) below.  (b) Example (41). Errors, estimator and
$\operatorname{IEI}$ are depicted on the top row and $\operatorname{EOC}$s and
energy of the reconstructed solution on the bottom row over time ($x$-axis).
Results are computed on the sequence of uniform meshes with mesh size $h$,
time-step $k=0.4h^{2}/(p+1)^{2}$, and $p=2$. The effect of using slightly
smaller time-step for fine $h$ (resulting from the $h^{2}$-term) is evident
compared to subfigure (a) above, in that the $\operatorname{EOC}\geq 2$ for
all times.

Fig. 3:  Example (39), violation of the CFL condition.  (a) Example (39).
Errors, estimator and $\operatorname{IEI}$ are depicted on the top row and
$\operatorname{EOC}$s and energy of the reconstructed solution on the bottom
row over time ($x$-axis). Results are computed on the sequence of uniform
meshes with mesh size $h$ and time step $k=2.0h/(p+1)^{2}$ and $p=3$. The
$\operatorname{IEI}$ behaviour indicates that the error is overestimated by
the estimator but follows the error behaviour. The method is unstable due to
the violation of the CFL condition (cf., Figure 2(b) for comparison with a
stable approximation.)


# Attention-based Spatial-Temporal Graph Neural ODE for Traffic Prediction

Weiheng Zhong
Department of Civil and Environmental Engineering
University of Illinois at Urbana-Champaign
Champaign, USA
<EMAIL_ADDRESS>
& Hadi Meidani
Department of Civil and Environmental Engineering
University of Illinois at Urbana-Champaign
Champaign, USA
<EMAIL_ADDRESS>
& Jane Macfarlane
Sustainable Energy Systems Group
Lawrence Berkeley National Laboratory
Berkeley, USA
<EMAIL_ADDRESS>

###### Abstract

Traffic forecasting is an important issue in intelligent traffic systems
(ITS). Graph neural networks (GNNs) are effective deep learning models to
capture the complex spatio-temporal dependency of traffic data, achieving
ideal prediction performance. In this paper, we propose attention-based graph
neural ODE (ASTGODE) that explicitly learns the dynamics of the traffic
system, which makes the prediction of our machine learning model more
explainable. Our model aggregates traffic patterns of different periods and
has satisfactory performance on two real-world traffic data sets. The results
show that our model achieves the highest accuracy of the root mean square
error metric among all the existing GNN models in our experiments.

_Keywords_ Graph neural ODE $\cdot$ Attention mechanism $\cdot$ Traffic
dynamics

##  1 Introduction

Accurate traffic prediction is of significant importance for the intelligent
transportation system, especially on the highway where a massive traffic flow
of high driving speed is moving. Traffic forecasting is a long-standing
challenge due to the complexity of spatio-temporal dependencies of traffic
data. On the one hand, traffic time series have strong nonlinear temporal
dynamics due to the varying traffic demand in different time points. On the
other hand, the correlation between different locations in a network changes
under different traffic conditions. Traffic accidents also cause non-
stationary traffic dynamics, making the traffic dynamics more unpredictable.

Graph neural networks (GNNs) attracted tremendous attention as a popular class
of deep learning models for traffic prediction. Unlike fully-connected neural
networks and convolution neural networks, GNNs can deal with graph-structured
data. GNNs outperform most traditional methods by aggregating traffic
information of neighbor locations and achieving satisfactory predictions of
future traffic conditions. There are many excellent GNN models, such as DCRNN
Li et al. (2017), STGCN Yu et al. (2017), ASTGCN Guo et al. (2019),
GraphWaveNet Sun et al. (2018), STSGCN Song et al. (2020), Graph Autoformer
Zhong et al. (2022), Etc. These models have a specific framework to aggregate
the spatial and temporal traffic features to predict future traffic data.

Neural ODE Chen et al. (2018) appears as an inspiration for us to predict the
traffic future based on the dynamic of the traffic system. Spatial-temporal
neural ODE Zhou et al. (2021) combines neural ODE and convolution neural
network for urban area traffic prediction, leveraging spatial-temporal data
from multiple sources. Spatial-temporal graph neural ODE (STGODE) Fang et al.
(2021) combines graph neural network and neural ODE to extract a more extended
range of spatio-temporal correlation without being affected by the over-
smoothing problem. However, to the best of our knowledge, there is not much
literature that explicitly uses a deep learning model to learn the dynamics of
the traffic system and predicts the future traffic data through the system’s
evolution.

In this paper, we focus on using neural ODE to imitate the dynamics of the
traffic system by performing supervise learning on the output of neural ODE at
each time step. We utilize attention mechanism Feng et al. (2017) to model the
dynamics of the traffic system, using the spatial and temporal attention to
estimate the trend of traffic under different traffic conditions. We also
shows the advantage of aggregating traffic trend of different periods for
higher prediction accuracy. The adjoint training Chen et al. (2018) effects
are also discussed in this paper.

##  2 Problem Setup

In this study, we define a traffic network as an undirected graph $G=(V,E,S)$,
where $V$ is a set of all vertices and $E$ is a set of all edges of the
network; $S\in\mathbb{R}^{N\times N}$ defines the adjacent matrix of the graph
and $N=|V|$ is the number of all vertices.

We use sensors in the vertices to collect measurements of the traffic features
with a fixed sampling frequency within a period. The number of traffic
features each sensor collects is $F$. Then the collected traffic data sequence
$\bm{X}=\\{x_{1},x_{2},...,x_{t},x_{t+1},...,x_{T}\\}$, where
$x_{t}\in\mathbb{R}^{N\times F}$ denotes the record of all the traffic
features of all nodes at time $t$.

We denote $T_{h}$ as the time length of one hour, $T_{d}$ as the time length
of one day, and $T_{w}$ as the time length of one week. Then we define four
types of time series segments of the length of $T_{h}$:

  * •

The predicted segment $\bm{\chi}^{t_{p}}$ is defined as
$\\{x_{t_{p}+1},x_{t_{p}+2},...,x_{t_{p}+T_{h}}\\}$, where $x_{t_{p}+1}$ is
the starting point of the predicted segment. We will use historical traffic
data segments to predict this traffic segment.

  * •

The recent segment $\bm{\chi}^{t_{r}}$ is defined as
$\\{x_{t_{p}-T_{h}+1},x_{t_{p}-T_{h}+2},...,x_{t_{p}}\\}$, where
$t_{r}=t_{p}-T_{h}$. It is the time segment that is temporally closest to the
predicted segment, which provides us with most important information for the
traffic condition of next one hour.

  * •

The daily period segment $\bm{\chi}^{t_{d}}$ is defined as
$\\{x_{t_{p}-T_{d}+1},x_{t_{p}-T_{d}+2},...,x_{t_{p}-T_{d}+T_{h}}\\}$, where
$t_{d}=t_{p}-T_{d}$. This is the traffic data of the same time in a day as the
predicted segment and we use it to capture the daily traffic pattern of the
system such as the rush hour in the morning and evening.

  * •

The weekly period segment $\bm{\chi}^{t_{w}}$ is defined as
$\\{x_{t_{p}-T_{w}+1},x_{t_{p}-T_{w}+2},...,x_{t_{p}-T_{w}+T_{h}}\\}$, where
$t_{w}=t_{p}-T_{w}$. This is the traffic data of the same time in a week as
the predicted segment and we use it to capture the weekly traffic pattern of
the system. Using this time segment we aims at detecting the dynamic of the
traffic system within the period of one week.

We intend to build a deep learning model using the historical traffic data
$\\{\bm{\chi}^{t_{r}},\bm{\chi}^{t_{d}},\bm{\chi}^{t_{w}}\\}$ as inputs to
predict future traffic data $\bm{\chi}^{t_{p}}$. An example of constructing
input data and output data is shown in Figure 1.

Figure 1: An example of the input and output traffic segments of our model is
shown. The recent traffic segment is defined as the traffic conditions of the
previous hour. The daily traffic segment is the traffic conditions of
yesterday at the same time slot as the predicted traffic segment. Similarly,
The weekly traffic segment is the traffic conditions of last week at the same
time slot of the same day as the predicted traffic segment.

##  3 Methodology

###  3.1 Model framework

Our model consists of three independent Neural ODE blocks to capture the
dynamic of the traffic system of different periods. We use a fully-connected
network as our encoder to map the traffic feature data of each node in each
time point to a higher-dimensional space, which is employed to strengthen the
expressiveness of the model.

Then we use Neural ODE block to calculate the hidden state of the future time
segments by learning the traffic pattern of different time periods. We
implement three times "Forward Pass" Chen et al. (2018) of each Neural ODE
block to obtain $\\{H^{t+\Delta t},H^{t+2\Delta t},H^{t+3\Delta t}\\}$, where
$t$ is starting time point of the historical time segment and $\Delta t$ is
the $\frac{1}{3}$ time length of the whole period. This means that
$t_{w}+3\Delta t_{w}=t_{d}+3\Delta t_{d}=t_{r}+3\Delta t_{r}=t_{p}$. We
consider $H^{t_{p,w}},H^{t_{p,d}},H^{t_{p,r}}$ as the hidden feature of
predicted time segment of three independent ODE blocks.

We use a fully connected neural network as a fusion layer to aggregate the
hidden features $H^{t_{p,w}},H^{t_{p,d}},H^{t_{p,r}}$. The output of the
fusion layer $H^{t_{p}}$ is the same dimension as each hidden state of
$\\{H^{t_{p,w}},H^{t_{p,d}},H^{t_{p,r}}\\}$. As shown in Figure 2, we will
output a set of intermediate hidden states $\mathbb{H}=\\{\\{H^{t_{w}+\Delta
t_{w}},H^{t_{w}+2\Delta t_{w}},H^{t_{d}+\Delta t_{d}},H^{t_{d}+2\Delta
t_{d}},H^{t_{r}+\Delta t_{r}},$ $H^{t_{r}+2\Delta t_{r}}\\}$ and hidden state
of the predicted time segment $H^{t_{p}}$. We use a decoder, which is also a
fully-connected neural network, as a mapping from the hidden state to the
original data space. Combined with the decoder, We then can obtain the
prediction of the predicted time segment $X^{t_{p}}$.

We perform supervised training for our model by minimizing the mean square
error between predicted and ground-truth traffic data using stochastic
gradient descent algorithms. The training loss $L$ is defined as:

$L=MSE(\mathbb{X},\mathbb{X}_{GT})$

where $X^{t_{p}}_{GT}$ are the ground-truth time segment data. The details of
our model framework is shown in Figure 2.

Figure 2: The framework of the proposed attention-based spatial temporal graph
Neural ODE is shown in this figure. We apply different ODE blocks to different
input traffic segments, which are used to learn the traffic dynamics of
different time intervals. A fully-connected fusion layer is applied to the
aggregated features of three blocks and outputs the final aggregated features.

###  3.2 Neural ODE block

Each ODE block in Figure 2 consists of a spatial attention layer, a temporal
attention layer, and a chebyshev graph convolution layer.

The traffic conditions of one location can affect the traffic condition of
other neighbor locations. For this consideration, we adopt the chebyshev graph
convolution to aggregate the information of the neighbor vertices Kipf and
Welling (2016). In order to allow for capturing multi-hop neighbors’ traffic
conditions, we use up to third order of chebyshev polynomials.

However, the influence of neighbor vertices is highly dynamic. We may need to
pay different attention to the same neighbor vertices under different traffic
conditions when we aggregate the information of the neighbor vertices. Also,
in the temporal dimension, correlations between traffic conditions in
different time slices vary in different situations. Hence, we use the
attention mechanism Feng et al. (2017) to capture spatial and temporal
correlations by calculating the attention scores $A_{S}$ and $A_{T}$, using
the spatial and temporal attention layers. We couple all information in
chebyshev graph convolution layer to predict the change of traffic condition
$\frac{dH^{t}}{dt}$ over a specific time length $\Delta t$. The output of the
ODE block is the hidden state of future time segment $H^{t+\Delta
t}=H^{t}+\frac{dH^{t}}{dt}$. The details of the ODE block architecture is
shown in Figure 3.

Figure 3: The architecture of an individual neural ODE block is shown. We will
compute the spatial attention scores and temporal attention scores based on
the input traffic features. Combined with exact graph connectivity, we perform
diffusion graph convolution to derive the traffic features of the next time
step.

##  4 Numerical Experiments

###  4.1 Data set description

We use two real-world traffic data sets, PeMS-BAY Li et al. (2017) and PeMS04
Guo et al. (2019), to validate our model. These two data sets are collected by
the Caltrans Performance Measurement System (PeMS) every 30 seconds [STGODE].
Both the traffic data are aggregated into 5-minutes intervals. In PeMS-BAY
data set, we only have the data of traffic flow velocity. There are 325
sensors in Bay Area collecting six months of data from Jan $1^{st}$ 2017 to
May $31^{st}$ 2017\. As for PeMS04 data set, we have 307 sensors on the
highway of major metropolitan areas in California measuring three traffic
features, including traffic flow, average speed, and average occupancy, for
almost two months.

We use our model to predict traffic speed using PeMS-BAY data set and traffic
flow using PeMS04 data set. Both data sets are sorted by time from the past to
the present. We split the data into three parts for training (70%), validation
(10%), and testing (20%). Data normalization with the mean and standard
deviation of the training data is applied to these three parts.

###  4.2 Baseline models

We compare our model performance with several baseline models:

  * •

Historical average (HA) Ermagun and Levinson (2018): it estimates seansonal
traffic pattern and uses weighted average of traffic flow data as prediction.

  * •

Auto-regressive integrated moving average model (ARIMA) Ermagun and Levinson
(2018): it is a traditional parametric for analysis of time series data.

  * •

Fully connected LSTM (FC-LSTM) Sutskever et al. (2014): it is an RNN-based
sequence model whose encoder and decoder are both LSTM layers.

  * •

Diffusion convolutional recurrent neural network (DCRNN) Li et al. (2017): it
purposes a combination of diffusion convolution operator and GRU Fu et al.
(2016) to capture spatio-temporal correlations.

  * •

Spatial-temporal graph ODE (STGODE) Fang et al. (2021): aggregating the
spatial and temporal information of the traffic segments, it uses Neural ODE
to learn the system’s dynamic.

  * •

Attention based spatial temporal graph convolutional networks (ASTGCN) Guo et
al. (2019): it uses attention mechanism to capture spatial and temporal data
dependency for better future prediction.

  * •

Graph multi-attention network (GMAN) Zheng et al. (2020): it uses a multi-head
self-attention mechanism to capture the spatial and temporal correlation of
traffic series data and transform attention to capture dependency between
historical time segments and future time segments.

###  4.3 Comparison of traffic prediction performance

We use two different metrics to evaluate the model performance: the root mean
squared error (RMSE) and mean absolute error (MAE). We will exclude missing
data when evaluating the performance of both data sets. Using Adam optimizer,
we use a learning rate of 0.0001 and train our model for 50 epochs. The
tradeoff coefficient $\alpha$ is set to be 0.1, and the batch size is set to
be 32.

We report the traffic prediction error of 15-min ahead prediction, 30-min
ahead prediction, and 60-min ahead prediction in Table 1, using PeMS-BAY data
set for validation. We observe that our method outperforms all the baseline
models in 15-min ahead prediction and 30-min ahead prediction. For 60-min
ahead prediction, our model achieves significant improvement in decreasing
RMSE of the prediction and slight improvement in decreasing MAE. These results
demonstrate the effectiveness of learning traffic dynamics to predict future
traffic conditions.

Table 1: Traffic speed prediction performance of different models on PeMS-BAY data set model | 15 min | 30 min | 60 min
---|---|---|---
| RMSE | MAE | RMSE | MAE | RMSE | MAE
HA | 5.60 | 2.88 | 5.60 | 2.88 | 5.60 | 2.88
ARIMA | 3.30 | 1.61 | 4.76 | 2.33 | 6.50 | 3.39
FC-LSTM | 4.19 | 2.05 | 4.55 | 2.2 | 4.96 | 2.36
DCRNN | 2.95 | 1.38 | 3.97 | 1.74 | 4.74 | 2.05
ASTGCN | 2.80 | 1.35 | 3.82 | 1.66 | 4.56 | 2.06
GMAN | 2.82 | 1.34 | 3.72 | 1.62 | 4.32 | 1.86
Our method | 2.68 | 1.30 | 3.36 | 1.62 | 4.10 | 2.00

We also compare our model’s performance on the average results of traffic flow
prediction using PeMS04 data set over the next one hour in Table 2. The
results show that our model also outperforms other baseline models in traffic
flow prediction. Neural ODE-related model (STGODE and our model) has higher
performance than other baseline models. Moreover, it also shows that our model
has good performance in decreasing RMSE of the prediction.

Table 2: Traffic flow prediction performance of different models on PeMS04 data set model | RMSE | MAE
---|---|---
HA | 54.11 | 36.76
ARIMA | 68.13 | 32.15
FC-LSTM | 45.76 | 29.50
DCRNN | 37.61 | 24.63
ASTGCN | 35.22 | 22.93
STGODE | 32.82 | 20.84
Our method | 30.65 | 20.23

###  4.4 Effect of the fusion layer

We consider every component in our model has some specific physical meaning
except for the fusion layer. The encoder and decoder are two opposite mappings
between observed traffic data and the hidden state of the dynamic traffic
system. Neural ODE imitates the evolution of the dynamic system. However, it
is hard to explain the output of the fusion layer. So we will investigate the
effect of the fusion layer in our model.

We remove the fusion layer and use the prediction of three independent ODE
blocks for comparison. We consider that dynamic of different periods will have
different effects on our prediction. As shown in Figure 4, recent traffic data
has a significant effect on short-term traffic prediction (5min-20min), while
weekly period data is more important in long-term traffic prediction
(45min-60min). The ODE block of daily period data does not have ideal short-
term or long-term prediction performance. We consider that daily traffic
pattern is more difficult to capture due to abrupt change in rush hour traffic
conditions. However, by coupling features of different traffic patterns, we
can leverage different historical data effects on the prediction of different
time points, which helps us increase the model accuracy of future traffic
prediction.

Figure 4: Prediction error comparison of different neural ODE blocks and our
purposed model. We used RMSE to compare the performance of using different
traffic features. We observed that weekly period segment was useful for long-
term future prediction and recent traffic segment was beneficial for short-
term future prediction. The performance of our model showed the effectiveness
of the fusion layer to achieve higher performance in our model.

###  4.5 Effect of adjoint training method

The adjoint training method highly reduces the memory cost of Neural ODE
training Chen et al. (2018). However, Some researchers claimed that the
adjoint training method of Neural ODE introduces error into gradient
information, which will negatively influence the final well-trained model
Daulbaev et al. (2020). Here we also implement two training methods, adjoint
training and regular training. Then we compare the performance of models
trained by these two methods.

In Table 3, we conclude that the adjoint training method decreases prediction
accuracy, and GMAN has higher accuracy in MAE if we use the adjoint method to
train our model. However, the model of adjoint training still outperforms GMAN
in the metric of RMSE.

Table 3: Effect of adjoint training on prediction error. model | 15 min | 30 min | 60 min
---|---|---|---
| RMSE | MAE | RMSE | MAE | RMSE | MAE
GMAN | 2.82 | 1.34 | 3.72 | 1.62 | 4.32 | 1.86
adjoint training | 2.72 | 1.33 | 3.52 | 1.72 | 4.22 | 2.05
non-adjoint training | 2.68 | 1.30 | 3.36 | 1.62 | 4.10 | 2.00

We also plot the average prediction error of our models in each training epoch
in Figure 5. We can observe that the prediction error of the adjoint training
model has more significant vibration, which indicates that adjoint training is
less stable than non-adjoint training.

Figure 5: The effect of the adjoint method on training stability is shown. We
compared the error evolution over the training process under different
training methods. We observed that using adjoint training slightly reduced the
stability of the model training. However, it will not affect the final model
performance if we train the model with enough epochs.

Considering the factors of model performance and training stability, it is
suggestive that not use adjoint training unless dealing with an extremely
large neural ODE model.

##  5 Conclusion

In this paper, we propose a novel deep learning architecture to capture the
dynamics of the traffic network system. We help our Neural ODE better imitate
the evolution of the traffic system by introducing the attention mechanism to
capture the spatial and temporal correlation between traffic data. We also
purpose a fusion layer to aggregate features of different periodic dynamics to
perform more accurate predictions for the future traffic data. The results
show that our model can outperform most of the existing models in the root
mean square error metric.

However, the attention mechanism is not the exact dynamics of the traffic
system. To make the model more explainable, we can introduce more physics-
related components in the neural ODE block, such as diffusion processes. Also,
adjoint training is an efficient computation method of training our model. If
we can increase the stability of the adjoint training, we can construct a more
complicated deep learning model to capture the temporal dependency better.

## References

  * Li et al. [2017] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _arXiv preprint arXiv:1707.01926_ , 2017.
  * Yu et al. [2017] Bing Yu, Haoteng Yin, and Zhanxing Zhu.  Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting.  _arXiv preprint arXiv:1709.04875_ , 2017.
  * Guo et al. [2019] Shengnan Guo, Youfang Lin, Ning Feng, Chao Song, and Huaiyu Wan.  Attention based spatial-temporal graph convolutional networks for traffic flow forecasting.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 33, pages 922–929, 2019.
  * Sun et al. [2018] Yang Sun, Zuoqin Shi, Hanhan Deng, Wenquan Wu, Ling Zhao, Sumin Li, Zhenshi Zhang, and Ronghua Du.  Identifying important spatial structures in cities based on graphwave–taking shanghai as an example.  In _2018 11th International Conference on Intelligent Computation Technology and Automation (ICICTA)_ , pages 209–213. IEEE, 2018.
  * Song et al. [2020] Chao Song, Youfang Lin, Shengnan Guo, and Huaiyu Wan.  Spatial-temporal synchronous graph convolutional networks: A new framework for spatial-temporal network data forecasting.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 914–921, 2020.
  * Zhong et al. [2022] Weiheng Zhong, Tanwi Mallick, Hadi Meidani, Jane Macfarlane, and Prasanna Balaprakash.  Explainable graph pyramid autoformer for long-term traffic forecasting, 2022.
  * Chen et al. [2018] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018.
  * Zhou et al. [2021] Fan Zhou, Liang Li, Kunpeng Zhang, and Goce Trajcevski.  Urban flow prediction with spatial–temporal neural odes.  _Transportation Research Part C: Emerging Technologies_ , 124:102912, 2021.
  * Fang et al. [2021] Zheng Fang, Qingqing Long, Guojie Song, and Kunqing Xie.  Spatial-temporal graph ode networks for traffic flow forecasting.  In _Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining_, pages 364–373, 2021.
  * Feng et al. [2017] Xiaocheng Feng, Jiang Guo, Bing Qin, Ting Liu, and Yongjie Liu.  Effective deep memory networks for distant supervised relation extraction.  In _IJCAI_ , volume 17, 2017.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Ermagun and Levinson [2018] Alireza Ermagun and David Levinson.  Spatiotemporal traffic forecasting: review and proposed directions.  _Transport Reviews_ , 38(6):786–814, 2018.
  * Sutskever et al. [2014] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  _Advances in neural information processing systems_ , 27, 2014.
  * Fu et al. [2016] Rui Fu, Zuo Zhang, and Li Li.  Using lstm and gru neural network methods for traffic flow prediction.  In _2016 31st Youth Academic Annual Conference of Chinese Association of Automation (YAC)_ , pages 324–328. IEEE, 2016.
  * Zheng et al. [2020] Chuanpan Zheng, Xiaoliang Fan, Cheng Wang, and Jianzhong Qi.  Gman: A graph multi-attention network for traffic prediction.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 1234–1241, 2020.
  * Daulbaev et al. [2020] Talgat Daulbaev, Alexandr Katrutsa, Larisa Markeeva, Julia Gusak, Andrzej Cichocki, and Ivan Oseledets.  Interpolated adjoint method for neural odes.  _arXiv preprint arXiv:2003.05271_ , 2020.