# Explicit Exactly Energy-conserving Methods for Hamiltonian Systems

Stefan Bilbao  Acoustics and Audio Group, University of Edinburgh, Edinburgh,
UK Michele Ducceschi  Department of Industrial Engineering, University of
Bologna, Bologna, Italy  Fabiana Zama  Department of Mathematics, University
of Bologna, Bologna, Italy

###### Abstract

For Hamiltonian systems, simulation algorithms that exactly conserve numerical
energy or pseudo-energy have seen extensive investigation. Most available
methods either require the iterative solution of nonlinear algebraic equations
at each time step, or are explicit, but where the exact conservation property
depends on the exact evaluation of an integral in continuous time. Under
further restrictions, namely that the potential energy contribution to the
Hamiltonian is non-negative, newer techniques based on invariant energy
quadratisation allow for exact numerical energy conservation and yield
linearly implicit updates, requiring only the solution of a linear system at
each time step. In this article, it is shown that, for a general class of
Hamiltonian systems, and under the non-negativity condition on potential
energy, it is possible to arrive at a fully explicit method that exactly
conserves numerical energy. Furthermore, such methods are unconditionally
stable, and are of comparable computational cost to the very simplest
integration methods (such as Störmer-Verlet). A variant of this scheme leading
to a conditionally-stable method is also presented, and follows from a
splitting of the potential energy. Various numerical results are presented, in
the case of the classic test problem of Fermi, Pasta and Ulam, as well as for
nonlinear systems of partial differential equations, including those
describing high amplitude vibration of strings and plates.

##  1 Introduction

The design of exact energy conserving numerical methods for nonlinear
Hamiltonian systems goes back as least as far as the work of LaBudde and
Greenspan, for the single particle subject to a central potential [1],
Marciniak [2] for the $N$-body problem, Strauss and Vazquez [3] for the
nonlinear Klein-Gordon equation, and Greenspan on the nonlinear harmonic
oscillator [4], and was greatly generalized by Simo et al. [5] in their work
on energy- and momentum-conserving algorithms. The concept of energy
conservation in time stepping schemes goes back much further in the case of
linear systems—see the discussion and references in [5], particularly with
regard to the energy conservation of the Crank-Nicholson method for linear
systems. An exact numerical conservation of a “pseudo-energy” (henceforth
simply “energy” or “numerical energy” in this article) of a numerical scheme
for a conservative system may be contrasted with schemes that conserve energy
approximately—see, e.g., early work by Hughes et. al [6]. The same point of
view is taken with other geometrical numerical integration techniques, such
as, e.g., symplectic methods, where energy is conserved approximately,
following the result from Zhong and Marsden [7] that, under some restrictions,
exact energy conservation and symplecticity for Hamiltonian systems cannot be
obtained simultaneously [8]. Various studies examine the degree to which
energy is conserved for symplectic and non-symplectic methods [9, 10]. In this
paper, we focus on exact conservation of numerical energy rather than
approximate conservation, and we will not consider momentum conservation. Part
of the reason for the focus here on exact numerical energy conservation is the
possibility of determining conditions for numerical stability as has been
pointed out earlier—see [5] and Remark 8 in [11].

Most exact numerical energy-conserving algorithms presented to date are
implicit [3, 2, 5, 12], requiring the solution of a nonlinear system of
algebraic equations at each time step, normally through an iterative method
(such as, e.g., Newton-Raphson). This is obviously a computational bottleneck,
and a fully explicit formulation is preferable. Note, however, that for linear
systems, explicit and exact energy-conserving methods are available—see, e.g.,
[13]. Exact explicit and conservative intergrators have also been designed on
a case-by-case basis—normally new ad hoc stability considerations appear in a
problem-dependent way regarding the choice of the time step. See [14] for a
treatment of the three-body problem, and [15] for the Kepler problem.

In a recent article by Marazzato et al. [11], an explicit Hamiltonian
integrator is presented that preserves a numerical energy exactly, for general
nonlinear systems. Two features are worth noting: a) the conserved quantity is
defined through a continuous time integration, and thus must itself be
approximated, leading effectively to approximate (non-exact) energy
conservation in a discrete time implementation—with fine enough approximation,
numerical energy conservation to machine precision may be achieved; and b)
when the energy of the model system is non-negative, as it commonly is in
practice, the numerical energy does not inherit this property, and thus
conditions for numerical stability do not immediately follow. Other exactly
conservative methods, also reliant on the exact evaluation of a continuous
integration have been proposed previously [16, 17].

In the context of gradient flows for diffusive systems, recent progress has
been made in developing energy-stable methods through a variable
transformation representing the energy, generally a nonlinear function of the
state, as a quadratic form—such methods are referred to as invariant energy
quadratisation (IEQ) approaches [18, 19]. The main result is that it is
possible to arrive at updating equations for a time-stepping method that are
linearly implicit—so that the update depends only on the solution of a linear
system, rather than the solution of a system of nonlinear algebraic equations.
This allows for the sidestepping of the many difficulties associated with
iterative solvers, and reduced computational cost. Alongside the more recently
introduced scalar auxiliary variable (SAV) approaches [20], such techniques
have been applied to a wide variety of problems [21]. More recently, IEQ/SAV
approaches have been applied to Hamiltonian systems, as in the case of
diagonally-implicit Runge Kutta methods [22]. For an interesting overview of
the relationship between quadratisation techniques and linearly implicit
schemes, see the recent article by Sato et al. [23].

A linearly implicit method will require the solution of a linear system at
each time step, and will most likely constitute the computational bottleneck
in a simulation as a whole. As will be shown here, using IEQ and SAV
approaches, it is possible to arrive at energy-stable (indeed exactly lossless
in machine arithmetic) numerical designs that are fully explicit. The explicit
character of the update derives from the availability of a closed form inverse
for the linearly implicit system that arises, and furthermore the ability to
solve an $N\times N$ linear system in $O(N)$ operations, through the Sherman-
Morrison inversion theorem [24]. This property is fully general, and not
dependent on the particular form of the Hamiltonian, except through the
additional non-negativity requirement on the potential energy. As a result,
computer execution time for exactly energy conserving methods for Hamiltonian
systems is very nearly on par with that of the simplest explicit numerical
schemes, such as, e.g., Störmer-Verlet integration. This design addresses the
two points a) and b) above with reference to the scheme presented in [11].

In Section 2, a restricted class of Hamiltonian systems is introduced. Under
the constraint that the potential energy is non-negative, quadratisation
techniques can be used to arrive at a simplified system of three equations in
momentum, position, and a single scalar auxiliary variable derived from the
potential energy. A further generalised form is also discussed, where a
quadratic term is extracted or split from the expression for the potential
energy before quadratisation, leading to a distinct starting point for
numerical designs. Time stepping methods are introduced in Section 3,
including first the rudimentary Störmer-Verlet method, the energy-conserving
method presented in [11], and then an exact energy-conserving and
unconditionally stable method obtained using IES/SAV approaches, and allowing
an explicit update. A variation of this last scheme to the case of a split
potential energy is also exactly energy-conserving, and leads to a distinct
conditionally-stable method. In Section 4, three examples are presented: the
Fermi-Pasta-Ulam ODE system, and two PDE systems then reduced, by semi-
discretisation, to Hamiltonian ODE systems: the nonlinear vibration of a
string, and the high-amplitude vibration of a thin plate. Various numerical
results are presented, illustrating convergence rates, exact numerical energy
conservation to machine precision, and relative computation times. Some
concluding remarks appear in Section 5.

Preliminary results have been presented recently at the 2022 European
Nonlinear Dynamics Conference [25].

##  2 Hamiltonian Systems

Consider a system of particles, with $N$ generalised coordinates ${\bf
q}=[q_{1},\ldots,q_{N}]^{\intercal}$ and momenta ${\bf
p}=[p_{1},\ldots,p_{N}]^{\intercal}$. Here, ⊺ indicates a transposition
operation, so ${\bf p}$ and ${\bf q}$ are $N\times 1$ column vectors. Both are
functions of time $t\geq 0$, so ${\bf p}={\bf p}(t)$ and ${\bf q}={\bf q}(t)$.
Suppose also that an associated Hamiltonian $H({\bf p},{\bf q})$ is defined by

$H\left({\bf p},{\bf q}\right)=\tfrac{1}{2}{\bf p}^{\intercal}{\bf M}^{-1}{\bf
p}+{V}\left({\bf q}\right)$ (1)

for some constant symmetric positive definite $N\times N$ matrix ${\bf M}$
referred to as the mass matrix. $V({\bf q})$ is the potential energy for the
system of particles. $H$ and $V$ are scalar functions of time $t$, through
their dependence on ${\bf p}$ and ${\bf q}$.

The system dynamics follow from Hamilton’s equations:

$\dot{{\bf q}}=\nabla_{{\bf p}}H\qquad\qquad\dot{{\bf p}}=-\nabla_{{\bf
q}}H\,,$ (2)

where dots indicate differentiation with respect to time $t$, and
$\nabla_{{\bf p}}$ and $\nabla_{{\bf q}}$ are gradients with respect to ${\bf
p}$ and ${\bf q}$, respectively. For the particular form of Hamiltonian given
in (1), Hamilton’s equations become

$\dot{{\bf q}}={\bf M}^{-1}{\bf p}\qquad\qquad\dot{{\bf p}}=-\nabla_{{\bf
q}}V\,.$ (3)

Through time differentiation of the first of (3), and substitution of the
second, a second order system of ordinary differential equations in ${\bf q}$
results:

${\bf M}\ddot{{\bf q}}=-\nabla_{{\bf q}}V\,.$ (4)

Such a second order form serves as the starting point for many numerical
methods, including the classic Störmer-Verlet method. See Section 3.1. In the
first order representation (3), initial conditions are required for ${\bf q}$
and ${\bf p}$, so ${\bf q}(0)={\bf q}^{(0)}$ and ${\bf p}(0)={\bf p}^{(0)}$,
for given $N$ vectors ${\bf q}^{(0)}$ and ${\bf p}^{(0)}$. For the second
order form (4), one may set $\dot{{\bf q}}(0)={\bf M}^{-1}{\bf p}^{(0)}$.

The defining feature of a Hamiltonian system is that, given initial conditions
${\bf p}^{(0)}$ and ${\bf q}^{(0)}$, the energy $H(t)$ is constant for all
$t\geq 0$:

$\dot{H}=0\quad\rightarrow\qquad H(t)=H(0)={\rm constant}\quad\forall t\geq
0\,,$ (5)

where $H(0)=H({\bf p}^{(0)},{\bf q}^{(0)})$, evaluated from (1) using initial
conditions ${\bf q}^{(0)}$ and ${\bf p}^{(0)}$.

###  2.1 Comments

The system and Hamiltonian in (1) are not the most general possible. A few
remarks are offered here:

  * •

The mass matrix ${\bf M}$ is assumed constant—a common assumption [5]. In many
cases of interest, it is furthermore diagonal, but we will consider the more
general form here, and indicate cases in which a diagonal form of ${\bf M}$
will have an impact on computational performance.

  * •

The Hamiltonian in (1) is separable, so that $H({\bf p},{\bf q})=T({\bf
p})+V({\bf q})$, for kinetic energy $T$ and potential energy $V$. In some
cases of interest [8], the mass matrix ${\bf M}$ may be of the form ${\bf
M}={\bf M}({\bf q})$, disturbing this splitting, but these cases will not be
considered here.

  * •

Following from the points above, the Hamiltonian is quadratic in ${\bf p}$,
meaning that the associated dynamical system is linear in ${\bf p}$.

  * •

The individual displacements $q_{i}$ and momenta $p_{i}$, $i=1,\ldots,N$, are
assumed scalar here, and represent either displacements/momenta in a signal
coordinate direction, or individual components of more general vector
displacements/momenta. Both cases will be seen in the numerical examples in
Section 4.

###  2.2 Non-negativity of the Potential Energy

As a further constraint, assume that

$V({\bf q})\geq 0\qquad\forall{\bf q}\in{\mathbb{R}}^{N}\,.$ (6)

This constraint is a natural one in many applications, but not all—for
example, the gravitational potential used in the calculations of planetary
orbits is not of this form [26]. This further implies, from (1), that

$H\geq 0\qquad\forall{\bf p},{\bf q}\in{\mathbb{R}}^{N}\,.$ (7)

Furthermore, as ${\bf M}>0$, from (5), one has the following bound on ${\bf
p}(t)$ in terms of the initial conditions:

$\|{\bf p}(t)\|\leq\sqrt{2\lambda_{{\rm max}}({\bf M})H(0)}\qquad\forall t\geq
0\,,$ (8)

where $\lambda_{{\rm max}}({\bf M})$ is the maximal eigenvalue of ${\bf M}$,
and $\|\cdot\|$ indicates a Euclidean norm. If $V({\bf q})$ is radially
unbounded [27], so that $V({\bf q})\rightarrow+\infty$ as $\|{\bf
q}\|\rightarrow+\infty$, then a further bound follows for $\|{\bf q}\|$.

The non-negativity property of $V$ is essential to invariant energy
quadratisation methods—indeed, it can be generalized to the case of $V$
bounded from below, so that $V({\bf q})\geq c$, for any constant $c$ [20], but
the non-negativity condition (6) above will be used here for simplicity.

###  2.3 Potential Energy Quadratisation

Under the non-negativity condition on $V$, from (6), one may define

$V=\tfrac{1}{2}\psi^{2}\,.$ (9)

Hamilton’s equations become

$\dot{{\bf q}}={\bf M}^{-1}{\bf p}\qquad\qquad\dot{{\bf p}}=-\psi\nabla_{{\bf
q}}\psi\,.$ (10)

Furthermore, using the chain rule,

$\dot{\psi}=(\nabla_{{\bf q}}\psi)^{\intercal}\dot{{\bf q}}\,.$ (11)

Finally, introducing

${\bf g}\triangleq\nabla_{{\bf q}}\psi=\frac{1}{\sqrt{2V}}\nabla_{{\bf
q}}V\,,$ (12)

a system of three equations results:

$\dot{{\bf q}}={\bf M}^{-1}{\bf p}\qquad\qquad\dot{{\bf p}}=-\psi{\bf
g}\qquad\qquad\dot{\psi}={\bf g}^{\intercal}\dot{{\bf q}}\,.$ (13)

An auxiliary initial condition $\psi(0)$ may be set as $\psi(0)=\sqrt{2V({\bf
q}^{(0)})}$ in terms of the given initial condition ${\bf q}(0)={\bf
q}^{(0)}$.

The key feature of system (13) is that, if ${\bf g}$ is assumed known at at
any given time instant, the resulting equations are linear in ${\bf p}$, ${\bf
q}$ and $\psi$. Though in the present case, ${\bf g}$ does indeed have a
dependence on ${\bf q}$, in the numerical setting, ${\bf g}$ can be evaluated
directly using previously computed values of the solution, and thus the update
becomes a matter of solving a linear system. This is the essence of the
“linearly implicit” property of schemes arising from energy quadratisation. As
will be shown subsequently here, under the appropriate numerical design, the
linear system is of a particularly simple form with a known easily-computed
inverse leading, effectively, to an explicit update. Such a property is
independent of the particular form of the potential energy $V({\bf q})$,
provided the non-negativity constraint (6) is satisfied.

###  2.4 Potential Energy Splitting

In some cases, the expression for the potential energy naturally takes the
form

$V({\bf q})=\frac{1}{2}{\bf q}^{\intercal}{\bf K}{\bf q}+V^{\prime}({\bf
q})\,,$ (14)

where ${\bf K}$ is a symmetric positive semi-definite $N\times N$ matrix, and
where $V^{\prime}\geq 0$. The first term could represent the stored energy of
the system due to linear mechanisms, and $V^{\prime}$ that due to additional
nonlinear effects. Now, using

$V^{\prime}=\tfrac{1}{2}\psi^{2}\,,$ (15)

one arrives at a system of equations generalizing (13):

$\dot{{\bf q}}={\bf M}^{-1}{\bf p}\qquad\qquad\dot{{\bf p}}=-{\bf K}{\bf
q}-\psi{\bf g}\qquad\qquad\dot{\psi}={\bf g}^{\intercal}\dot{{\bf q}}\,.$ (16)

Though equivalent to (13), in a numerical setting, such a splitting allows for
a larger family of numerical designs, treating the linear and nonlinear parts
of the problem separately. See Section 3.4.

###  2.5 Potential Energy Gauge and Regularisation

A well-known technique in IEQ/SAV approaches is the introduction of a constant
shift of the global energy (see e.g. [20, 28])—Hamilton’s equations (3) are
unchanged (i.e. gauge invariant), but some numerical schemes, such as IEQ/SAV-
based approaches are affected [29]. This amounts to the replacement

$V\rightarrow V+\epsilon$ (17)

for a suitably chosen shift constant $\epsilon\geq 0$. Under the non-
negativity constraint (6), $V+\epsilon$ is thus bounded away from zero,
regularizing the calculation of ${\bf g}$ as defined in (12). The
regularization approach applies equally to the case of a split potential,
under the replacement $V^{\prime}\rightarrow V^{\prime}+\epsilon$.

##  3 Numerical Methods

In this section, we assume time discretization using a constant time step $k$,
such that solutions are computed at times $t^{n}=nk$, $n=0,1,\ldots$. A
discrete-time vector ${\bf u}^{n}$ represents an approximation to a continuous
time vector ${\bf u}(t)$ at times $t=t^{n}$.

###  3.1 Störmer-Verlet Integration

The most basic approach to the integration of the Hamiltonian system as
defined in Section 2 is through direct approximation of the second order form
in (4), using centered differences, as:

${\bf q}^{n+1}=2{\bf q}^{n}-{\bf q}^{n-1}-k^{2}{\bf M}^{-1}\nabla_{{\bf
q}}V|_{{\bf q}={\bf q}^{n}}\,.$ (18)

This discretisation is referred to as Störmer-Verlet, and was known to
Newton—see [30] and the references therein. It is a fully explicit two-step
scheme—in order to advance the solution to time step $n+1$, a direct
evaluation of the gradient of the potential energy $V$ at time step $n$ is
required. The values of the state ${\bf q}^{n}$ at $n=0$ and $n=1$ may be
initialised, using the initial conditions ${\bf q}^{(0)}$ and ${\bf p}^{(0)}$,
as

${\bf q}^{0}={\bf q}^{(0)}\qquad{\bf q}^{1}={\bf q}^{(0)}+k{\bf M}^{-1}{\bf
p}^{(0)}$ (19)

to second order in $k$, using a Taylor series approximation.

Störmer-Verlet is second-order accurate in the time step $k$, time-reversible,
and symplectic, but not energy-conserving except in an approximate sense [30].
It can exhibit numerical instability depending on both the choice of time step
and the size of the initial conditions. Such instabilities will be mentioned
in Section 4.

###  3.2 Explicit Energy-conserving Method (Marazzato et al.)

Consider the following scheme approximating (3), proposed by Marazzato et al.
in [11]:

${\bf q}^{n+1}={\bf q}^{n}+k{\bf M}^{-1}{\bf
p}^{n+\tfrac{1}{2}}\qquad\qquad{\bf p}^{n+\tfrac{3}{2}}={\bf
p}^{n-\tfrac{1}{2}}-2\int_{nk}^{(n+1)k}\nabla_{{\bf q}}V(\tilde{{\bf
q}}^{n}(t))dt\,.$ (20)

This is a time interleaved scheme, with ${\bf p}^{n}$ and ${\bf
q}^{n+\tfrac{1}{2}}$ defined at alternating multiples of $k/2$. The first of
(20) is a standard interleaved approximation to the first of (3). The second
of (20) relies on a continuous integration of $\nabla_{{\bf q}}V$ over the
time interval $t\in[nk,(n+1)k]$, and over the known free-flight trajectory
$\tilde{{\bf q}}^{n}(t)$ defined by

$\tilde{{\bf q}}^{n}(t)={\bf q}^{n}+(t-nk){\bf M}^{-1}{\bf
p}^{n+\tfrac{1}{2}}\,.$ (21)

Other conservative schemes proposed also rely on such a continuous integration
[16]. Except for particular functional forms of the potential energy $V({\bf
q})$, this integral is not available in closed form and must be approximated.
However, once the integral is evaluated, the scheme (20) is fully explicit.

The scheme (20) has an associated numerical energy that is exactly conserved
[11]:

${\bf H}^{n}=\tfrac{1}{2}{\bf p}^{n+\tfrac{1}{2}}{\bf M}^{-1}{\bf
p}^{n-\tfrac{1}{2}}+V({\bf q}^{n})={\rm constant}\,.$ (22)

There are two important points to mention here:

  * •

Scheme (20) is exactly conservative, but depends on a continuous integration
over a free-flight trajectory as given in (21) in order to achieve this
property. Thus the exact conservation property is approached in the limit of
increasing accuracy in the approximation of the continuous integration. It is
also important to point out here that a fine-grained approximation of the
continuous integration will require multiple evaluations of the gradient
$\nabla_{{\bf q}}V$ which, depending on the functional form of $V$, represents
an additional cost that grows with the desired accuracy of the approximation.

  * •

When $V\geq 0$, the numerical energy (22) does not inherit the non-negativity
property of the model system, and thus cannot be used in order to bound
solution growth, as pointed out in Remark 8 of [11].

Under the additional non-negativity constraint (6) on $V$, it is possible to
demonstrate a method that conserves a pseuodenergy that does not depend on a
continuous integration, and that furthermore inherits the non-negativity
property of the Hamiltonian of the model system, allowing for useful bounds on
solution size. Furthermore, through the exploitation of matrix structure, it
will be shown that such a method is fully explicit. Additionally, only one
function evaluation $\nabla_{{\bf q}}V$ is required per time step.

###  3.3 Explicit Exactly Energy-conserving Method

Returning now to the form (13) of Hamilton’s equations obtained under
quadratisation of the potential energy, and the introduction of the new
variable $\psi$ as in (9), consider the following scheme, written in terms of
the interleaved time series ${\bf p}^{n+\tfrac{1}{2}}$ and ${\bf q}^{n}$, and
the scalar time series $\psi^{n+\tfrac{1}{2}}$:

$\displaystyle{\bf q}^{n+1}$ $\displaystyle=$ $\displaystyle{\bf q}^{n}+k{\bf
M}^{-1}{\bf p}^{n+\tfrac{1}{2}}$ (23a) $\displaystyle{\bf p}^{n+\tfrac{1}{2}}$
$\displaystyle=$ $\displaystyle{\bf p}^{n-\tfrac{1}{2}}-\tfrac{k}{2}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)$ (23b)
$\displaystyle\psi^{n+\tfrac{1}{2}}$ $\displaystyle=$
$\displaystyle\psi^{n-\tfrac{1}{2}}+\tfrac{1}{2}\left({\bf
g}^{n}\right)^{\intercal}\left({\bf q}^{n+1}-{\bf q}^{n-1}\right)\,.$ (23c)

Here, ${\bf g}^{n}$ is defined as

${\bf g}^{n}=\nabla_{{\bf q}}\psi|_{{\bf q}={\bf
q}^{n}}=\frac{1}{\sqrt{2V({\bf q}^{n})}}\nabla_{{\bf q}}V|_{{\bf q}={\bf
q}^{n}}\,.$ (24)

The system (23) may be manipulated into an explicit update form in the
following way. Beginning from (23a), one may write:

$\displaystyle{\bf q}^{n+1}$
$\displaystyle\stackrel{{\scriptstyle\eqref{fd1_1}}}{{=}}$ $\displaystyle
2{\bf q}^{n}-{\bf q}^{n-1}+k{\bf M}^{-1}\left({\bf p}^{n+\tfrac{1}{2}}-{\bf
p}^{n-\tfrac{1}{2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{fd1_2}}}{{=}}$ $\displaystyle
2{\bf q}^{n}-{\bf q}^{n-1}-\tfrac{k^{2}}{2}{\bf M}^{-1}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)$
$\displaystyle\stackrel{{\scriptstyle\eqref{fd1_3}}}{{=}}$ $\displaystyle
2{\bf q}^{n}-{\bf q}^{n-1}-k^{2}{\bf M}^{-1}{\bf
g}^{n}\psi^{n-\tfrac{1}{2}}-\tfrac{k^{2}}{4}{\bf M}^{-1}{\bf g}^{n}({\bf
g}^{n})^{\intercal}\left({\bf q}^{n+1}-{\bf q}^{n-1}\right)\,.$

Finally, the update has the form

${\bf A}^{n}{\bf q}^{n+1}={\bf b}^{n}\,,$ (26)

where

${\bf A}^{n}={\bf
I}+{\bm{\alpha}}^{n}({\bm{\beta}}^{n})^{\intercal}\qquad\qquad{\bf
b}^{n}=2{\bf q}^{n}-2k{\bm{\alpha}}^{n}\psi^{n-\tfrac{1}{2}}-\left({\bf
I}-{\bm{\alpha}}^{n}({\bm{\beta}}^{n})^{\intercal}\right){\bf q}^{n-1}\,.$
(27)

Here, ${\bf I}$ is the $N\times N$ identity matrix, and the vectors
${\bm{\alpha}}^{n}$ and ${\bm{\beta}}^{n}$ are defined in terms of ${\bf
g}^{n}$ by ${\bm{\alpha}}^{n}=\tfrac{k}{2}{\bf M}^{-1}{\bf g}^{n}$ and
${\bm{\beta}}^{n}=\tfrac{k}{2}{\bf g}^{n}$. Thus, given ${\bf q}^{n}$, ${\bf
q}^{n-1}$ and $\psi^{n-\tfrac{1}{2}}$, both ${\bf A}^{n}$ and ${\bf b}^{n}$
may be explicitly constructed. Notice that the matrix ${\bf A}^{n}$ is
positive definite by construction, due to the positive definiteness of ${\bf
M}$, and thus the update (26) always has a unique solution. Once the update in
(26) has been performed, $\psi^{n+\tfrac{1}{2}}$ may be computed from (23c)
explicitly, using ${\bf q}^{n-1}$, ${\bf q}^{n+1}$ and
$\psi^{n-\tfrac{1}{2}}$. Scheme (26) requires initial values for ${\bf q}^{0}$
and ${\bf q}^{1}$, which may be set in the same way as for Störmer-Verlet, as
in (19), and also $\psi^{\tfrac{1}{2}}$ which may be set to second order in
$k$ in terms of the initial displacement ${\bf q}^{(0)}$ and momentum ${\bf
p}^{(0)}$ as

$\displaystyle\psi^{\tfrac{1}{2}}$ $\displaystyle=$
$\displaystyle\sqrt{2V({\bf q}^{(0)})}+\frac{k}{2\sqrt{2V({\bf
q}^{(0)})}}\left(\nabla_{{\bf q}}V|_{{\bf q}={\bf
q}^{(0)}}\right)^{\intercal}{\bf M}^{-1}{\bf p}^{(0)}\,.$ (28)

Though the solution of (26) apparently requires an $N\times N$ linear system
solution, in fact, the inverse is available in closed form, and allows for an
explicit solution in $O(N)$ operations. ${\bf A}^{n}$ is a rank-1 perturbation
of the identity, and thus Sherman Morrison inversion [24] yields:

$\left({\bf A}^{n}\right)^{-1}={\bf
I}-\frac{{\bm{\alpha}}^{n}({\bm{\beta}}^{n})^{\intercal}}{1+({\bm{\beta}}^{n})^{\intercal}{\bm{\alpha}}^{n}}\,.$
(29)

Computational cost is thus on par with the other methods presented here.
Notice in particular, though, that the solution of a linear system involving
${\bf M}^{-1}$ is common to all the methods in this section; the computational
cost of performing this operation, which may indeed be the bottleneck, is not
considered here. In many cases, though, ${\bf M}$ is diagonal, or even a
simple multiple of the identity, and thus its inversion is trivial.

Exact conservation of a numerical energy follows directly from scheme (23).
Left-multiplying (23b) by $\tfrac{1}{2}\left({\bf p}^{n+\tfrac{1}{2}}+{\bf
p}^{n-\tfrac{1}{2}}\right)^{\intercal}{\bf M}^{-1}$ gives:

$\displaystyle\tfrac{1}{2}\left({\bf p}^{n+\tfrac{1}{2}}+{\bf
p}^{n-\tfrac{1}{2}}\right)^{\intercal}{\bf M}^{-1}\left({\bf
p}^{n+\tfrac{1}{2}}-{\bf p}^{n-\tfrac{1}{2}}\right)+\tfrac{k}{4}\left({\bf
p}^{n+\tfrac{1}{2}}+{\bf p}^{n-\tfrac{1}{2}}\right)^{\intercal}{\bf
M}^{-1}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)\stackrel{{\scriptstyle\eqref{fd1_2}}}{{=}}$
$\displaystyle 0$ (30) $\displaystyle\tfrac{1}{2}\left({\bf
p}^{n+\tfrac{1}{2}}+{\bf p}^{n-\tfrac{1}{2}}\right)^{\intercal}{\bf
M}^{-1}\left({\bf p}^{n+\tfrac{1}{2}}-{\bf
p}^{n-\tfrac{1}{2}}\right)+\tfrac{1}{4}\left({\bf q}^{n+1}-{\bf
q}^{n-1}\right)^{\intercal}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)\stackrel{{\scriptstyle\eqref{fd1_1}}}{{=}}$
$\displaystyle 0$ $\displaystyle\tfrac{1}{2}\left({\bf
p}^{n+\tfrac{1}{2}}+{\bf p}^{n-\tfrac{1}{2}}\right)^{\intercal}{\bf
M}^{-1}\left({\bf p}^{n+\tfrac{1}{2}}-{\bf
p}^{n-\tfrac{1}{2}}\right)+\tfrac{1}{2}\left(\psi^{n+\tfrac{1}{2}}-\psi^{n-\tfrac{1}{2}}\right)\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)\stackrel{{\scriptstyle\eqref{fd1_3}}}{{=}}$
$\displaystyle 0\,.$

This may be identified as

$H^{n+\tfrac{1}{2}}=H^{n-\tfrac{1}{2}}={\rm constant}\,,$ (31)

where the numerical energy $H^{n+\tfrac{1}{2}}$ is

$H^{n+\tfrac{1}{2}}=\tfrac{1}{2}\left({\bf
p}^{n+\tfrac{1}{2}}\right)^{\intercal}{\bf M}^{-1}{\bf
p}^{n+\tfrac{1}{2}}+\tfrac{1}{2}\left(\psi^{n+\tfrac{1}{2}}\right)^{2}\geq
0\,.$ (32)

It is worth comparing this expression, for which non-negativity is ensured,
with the expression (22) for the scheme in Section 3.2. The momentum ${\bf
p}^{n+\tfrac{1}{2}}$ may be bounded in terms of the energy $H$, for all time
steps $n$, as

$\|{\bf p}^{n+\tfrac{1}{2}}\|\leq\sqrt{2\lambda_{{\rm max}}({\bf M})H}\,,$
(33)

which is identical to the bound (8) for the model system. The scheme (23) is
thus unconditionally numerically stable.

###  3.4 Split Potential Form

Consider now the split form of the potential energy $V$, as described in
Section 2.4, where a positive semi-definite quadratic form has been separated
from $V$ to leave a residual energy contribution $V^{\prime}\geq 0$, from
which an auxiliary variable $\psi$ may be defined as in (15). Consider the
following scheme, now modified with respect to (23) through the addition of a
linear term in (34b) below:

$\displaystyle{\bf q}^{n+1}$ $\displaystyle=$ $\displaystyle{\bf q}^{n}+k{\bf
M}^{-1}{\bf p}^{n+\tfrac{1}{2}}$ (34a) $\displaystyle{\bf p}^{n+\tfrac{1}{2}}$
$\displaystyle=$ $\displaystyle{\bf p}^{n-\tfrac{1}{2}}-k{\bf K}{\bf
q}^{n}-\tfrac{k}{2}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)$ (34b)
$\displaystyle\psi^{n+\tfrac{1}{2}}$ $\displaystyle=$
$\displaystyle\psi^{n-\tfrac{1}{2}}+\tfrac{1}{2}\left({\bf
g}^{n}\right)^{\intercal}\left({\bf q}^{n+1}-{\bf q}^{n-1}\right)\,.$ (34c)

where now,

${\bf g}^{n}=\nabla_{{\bf q}}\psi|_{{\bf q}={\bf
q}^{n}}=\frac{1}{\sqrt{2V^{\prime}({\bf q}^{n})}}\nabla_{{\bf
q}}V^{\prime}|_{{\bf q}={\bf q}^{n}}\,.$ (35)

An explicit update follows as in the case of the non-split form in (26), with
${\bf A}^{n}$ as given in (27), but with ${\bf b}^{n}$ now defined as

${\bf b}^{n}=\left(2{\bf I}-k^{2}{\bf M}^{-1}{\bf K}\right){\bf
q}^{n}-2k{\bm{\alpha}}^{n}\psi^{n-\tfrac{1}{2}}-\left({\bf
I}-{\bm{\alpha}}^{n}({\bm{\beta}}^{n})^{\intercal}\right){\bf q}^{n-1}\,.$
(36)

An expression for a conserved numerical energy follows as before, now taking
the form:

$H^{n+\tfrac{1}{2}}=\tfrac{1}{2}\left({\bf
p}^{n+\tfrac{1}{2}}\right)^{\intercal}{\bf M}^{-1}{\bf
p}^{n+\tfrac{1}{2}}+\tfrac{1}{2}({\bf q}^{n+1})^{\intercal}{\bf K}{\bf
q}^{n}+\tfrac{1}{2}\left(\psi^{n+\tfrac{1}{2}}\right)^{2}={\rm constant}\,.$
(37)

Because the second term is of indefinite sign, the numerical energy is not
necessarily non-negative. Because it is a quadratic form, however, it may be
bounded, using (34a), as:

$\tfrac{1}{2}({\bf q}^{n+1})^{\intercal}{\bf K}{\bf
q}^{n}\geq-\tfrac{1}{8}({\bf q}^{n+1}-{\bf q}^{n})^{\intercal}{\bf K}({\bf
q}^{n+1}-{\bf q}^{n})=-\tfrac{k^{2}}{8}\left({\bf
p}^{n+\tfrac{1}{2}}\right)^{\intercal}{\bf M}^{-1}{\bf K}{\bf M}^{-1}{\bf
p}^{n+\tfrac{1}{2}}\,.$ (38)

It then follows that

$H^{n+\tfrac{1}{2}}\geq\tfrac{1}{2}\left({\bf
p}^{n+\tfrac{1}{2}}\right)^{\intercal}\left({\bf M}^{-1}-\frac{k^{2}}{4}{\bf
M}^{-1}{\bf K}{\bf M}^{-1}\right){\bf p}^{n+\tfrac{1}{2}}\,.$ (39)

Given that ${\bf M}>0$, ${\bf M}^{-1}$ may be factored into unique upper and
lower triangular factors as ${\bf M}^{-1}={\bf M}^{-\tfrac{T}{2}}{\bf
M}^{-\tfrac{1}{2}}$, a condition for non-negativity of $H$ is

$k\leq\frac{2}{\lambda_{{\rm max}}\left({\bf M}^{-\tfrac{1}{2}}{\bf K}{\bf
M}^{-\tfrac{T}{2}}\right)}\,.$ (40)

The scheme is now conditionally stable, with the stability condition (40)
corresponding to that for a linear system when $V^{\prime}=0$.

###  3.5 Remarks

  * •

Linear Conditions: One of the interesting features of the non-split form (23)
is that, even if the potential $V$ is quadratic in ${\bf q}$, implying a
linear system, the scheme is not linear. If, however, the split form (34) is
used, then under linear conditions, $V^{\prime}=0$, and scheme (34) reduces
exactly to Störmer-Verlet (18). Note that under linear conditions, Störmer-
Verlet does indeed possess an exactly conserved numerical energy, as given in
(37) with $\psi^{n+\tfrac{1}{2}}=0$.

  * •

Generalized Splitting: The non-split form (23) and the split form (34) may be
combined to yield a large family of conservative schemes in an obvious way. If

$V({\bf q})=\tfrac{1}{2}{\bf q}^{\intercal}{\bf K}{\bf q}+V_{{\rm
nonlinear}}({\bf q})\,,$ (41)

where $V_{{\rm nonlinear}}\geq 0$ and consists of higher order terms in ${\bf
q}$, then any splitting of the form

$V({\bf q})=\tfrac{1}{2}{\bf q}^{\intercal}{\bf K}_{0}{\bf q}+V^{\prime}({\bf
q})\qquad V^{\prime}({\bf q})=V_{{\rm nonlinear}}+\tfrac{1}{2}{\bf
q}^{\intercal}{\bf K}^{\prime}{\bf q}\,,$ (42)

with ${\bf K}={\bf K}_{0}+{\bf K}^{\prime}$, ${\bf K}_{0}\geq 0$, ${\bf
K}^{\prime}\geq 0$ will yield a conservative form. The resulting numerical
scheme, with $V^{\prime}=\tfrac{1}{2}\psi^{2}$ will inherit conservation of
energy, which will be non-negative under a condition analogous to (40), but
depending on ${\bf K}_{0}$.

  * •

Generalized Update for ${\bf q}^{n}$: Consider again the non-split scheme
(23), which may be rewritten as

$\displaystyle{\bf q}^{n+1}$ $\displaystyle=$ $\displaystyle{\bf q}^{n}+k{\bf
M}^{-1}{\bf p}^{n+\tfrac{1}{2}}$ (43a) $\displaystyle{\bf p}^{n+\tfrac{1}{2}}$
$\displaystyle=$ $\displaystyle{\bf p}^{n-\tfrac{1}{2}}-\tfrac{k}{2}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)$ (43b)
$\displaystyle\psi^{n+\tfrac{1}{2}}$ $\displaystyle=$
$\displaystyle\psi^{n-\tfrac{1}{2}}+\tfrac{k}{2}\left({\bf
g}^{n}\right)^{\intercal}{\bf M}^{-1}\left({\bf p}^{n+\tfrac{1}{2}}+{\bf
p}^{n-\tfrac{1}{2}}\right)\,.$ (43c)

In this form, exact energy conservation follows from (43b) and (43c) only—it
is independent of the values of ${\bf g}^{n}$, which are derived solely from
${\bf q}^{n}$. Further opportunities for generalization are thus
available—(43a) could be replaced by any consistent update for ${\bf
q}^{n+1}$, and the exact numerical energy conservation property remains
undisturbed.

  * •

Variable Time Steps: Though the case of variable time steps will not be
discussed here in detail (and was indeed investigated in [11]), an exact
energy-conserving scheme follows immediately from the form given in (43)
above. Consider time instants $t^{n}$, and $t^{n+1/2}$, for integer $n$, where
$t^{n}<t^{n+1/2}<t^{n+1}$. From these, one may define two sequences of time
steps: $k^{n}=t^{n+1/2}-t^{n-1/2}$ and $k^{n+1/2}=t^{n+1}-t^{n}$. Supposing
that ${\bf q}^{n}$ is a time series defined for $t=t^{n}$, and similarly ${\bf
p}^{n+\tfrac{1}{2}}$ and ${\bf\psi}^{n+\tfrac{1}{2}}$ are defined for
$t=t^{n+\tfrac{1}{2}}$, then a scheme follows as:

$\displaystyle{\bf q}^{n+1}$ $\displaystyle=$ $\displaystyle{\bf
q}^{n}+k^{n+\tfrac{1}{2}}{\bf M}^{-1}{\bf p}^{n+\tfrac{1}{2}}$ (44a)
$\displaystyle{\bf p}^{n+\tfrac{1}{2}}$ $\displaystyle=$ $\displaystyle{\bf
p}^{n-\tfrac{1}{2}}-\tfrac{k^{n}}{2}{\bf
g}^{n}\left(\psi^{n+\tfrac{1}{2}}+\psi^{n-\tfrac{1}{2}}\right)$ (44b)
$\displaystyle\psi^{n+\tfrac{1}{2}}$ $\displaystyle=$
$\displaystyle\psi^{n-\tfrac{1}{2}}+\tfrac{k^{n}}{2}\left({\bf
g}^{n}\right)^{\intercal}{\bf M}^{-1}\left({\bf p}^{n+\tfrac{1}{2}}+{\bf
p}^{n-\tfrac{1}{2}}\right)\,.$ (44c)

It is direct to show that, regardless of the choices of $t^{n}$ and
$t^{n+\tfrac{1}{2}}$, the scheme (44) conserves the energy (32) exactly. Thus
such a generalisation is unconditionally stable, under the same reasoning as
in the case of constant time steps.

  * •

Regularisation: Regularisation through a potential energy shift, as described
in Section 2.5, impacts on the calculation of ${\bf g}$ in (24) and (35),
through a replacement $V\rightarrow V+\epsilon$ or $V^{\prime}\rightarrow
V^{\prime}+\epsilon$, respectively.

##  4 Examples

In this section, various Hamiltonian systems are simulated using schemes (23)
and (34), beginning with the Fermi-Pasta-Ulam ODE problem in Section 4.1, and
then progressing to more complex ODE systems derived as semi-discretisations
to PDE systems. These include the coupled transverse-longitudinal vibration of
a string at high amplitudes, in Section 4.2, and the high amplitude vibration
of a thin plate in Section 4.3.

In all cases, a useful measure of the exact energy conservation property is
the relative energy deviation error,

${\rm
Relative\,\,error}=\frac{H^{n+\tfrac{1}{2}}-H^{\tfrac{1}{2}}}{H^{\tfrac{1}{2}}}\,.$
(45)

In double precision floating point arithmetic, it is normally on the order of
machine precision, or approximately 10-16. Depending on the state size of the
system in question, however, larger deviations are possible.

###  4.1 Fermi-Pasta-Ulam Problem

As a simple example, consider the classic system of a linear arrangement
masses connected by linear and nonlinear springs, as proposed originally by
Fermi, Pasta and Ulam [31], and later adapted as a test problem by various
authors [32, 11], the form of which is followed here.

Consider a system of $N=2M$ masses in a linear arrangement, which longitudinal
displacements $q_{i}$ and momenta $p_{i}$, $i=1,\ldots,2M$. The system
Hamiltonian is defined by

$H=\frac{1}{2}\sum_{i=1}^{2M}p_{i}^{2}+V\qquad{\rm where}\qquad
V=\frac{\omega^{2}}{4}\sum_{i=1}^{M}\left(q_{2i}-q_{2i-1}\right)^{2}+\sum_{i=0}^{M}\left(q_{2i+1}-q_{2i}\right)^{4}\,,$
(46)

where, in the expression above, $q_{0}=q_{2M+1}=0$. Thus each mass is
connected, in an alternating arrangement, to a linear spring, and a cubic
nonlinear spring.

$V$ is clearly non-negative here, and thus, after consolidation of the
displacements $q_{i}$ and $p_{i}$ into vectors ${\bf p}$ and ${\bf q}$ of size
$N\times 1$, the scheme as presented in (23) follows immediately, with ${\bf
M}={\bf I}_{N}$, the $N\times N$ identity matrix, and with $V$ as given in
(46) above. This scheme is explicit, exactly conservative, and unconditionally
stable.

A split form (14) of the Hamiltonian also follows, using

${\bf K}=\frac{\omega^{2}}{2}{\bf
I}_{M}\otimes\begin{bmatrix}\,\,\,1&\\!\\!\\!-1\\\ -1&1\\\
\end{bmatrix}\qquad\qquad
V^{\prime}=\sum_{i=0}^{M}\left(q_{2i+1}-q_{2i}\right)^{4}\geq 0\,,$ (47)

where here, $\otimes$ indicates a Kronecker product. An exactly energy
conserving method follows, as in (34), and is stable under the condition (40),
which reduces in this case to

$k\leq\frac{2}{\omega}\,.$ (48)

####  4.1.1 Numerical Results

We use the same settings as in [11], and choose $\omega=50$, and $M=3$.
Simulations are run here with a time step of $k=10^{-3}$ s; the reference
solution is computed using Störmer-Verlet, with a time step of
$k=2^{-20}\approxeq 10^{-6}$ s. For initial conditions, we set ${\bf
p}^{(0)}={\bf 0}$, and ${\bf q}^{(0)}=[0,0,0,\alpha,0,0]^{\intercal}$, for
different values of $\alpha$. See Figure 1, illustrating c comparison between
time histories using the exactly conservative schemes (23) and (34) with
Störmer-Verlet (18). In all cases, the onset of errors is slightly faster for
(23) than for the other schemes, indicating a higher error (see also Figure
3).

Figure 1: Time histories of $q_{1}$, under the initial condition
$q_{4}=\alpha$, and under different values of $\alpha$ as indicated. Results
for schemes (23), (34) and (18) are plotted against a reference solution.

For the schemes (23) and (34), numerical energy, as defined by (32) and (37)
respectively, is conserved to machine precision. See Figure 2, illustrating
the relative deviation in energy, as defined in (45), for scheme (23) for the
Fermi-Pasta-Ulam system, under the conditions as described above, and for
$\alpha=100$. It is easily seen that the relative error is on the order of
machine roundoff error in double-precision floating point, or approximately
$10^{-16}$. Furthermore, it is possible to directly observe the quantisation
of the relative energy to the lowest bits in the machine number
representation.

Figure 2: Relative variation in numerical energy, as defined by (45), for
scheme (23) for the Fermi-Pasta-Ulam system.

As a basic test of convergence, consider the $L^{2}$ error defined, over a
simulation duration $n=0,\ldots,N_{f}$, by

${\rm Error}=\sqrt{\sum_{n=0}^{N_{f}}k\|{\bf q}^{n}-{\bf q}_{{\rm
ref}}(t=nk)\|^{2}}\,,$ (49)

where ${\bf q}_{{\rm ref}}$ are computed using Störmer-Verlet with a time step
of $k=2^{-20}\approxeq 10^{-6}$ s. The error appears in Figure 3, for the
exact energy conserving method in (23), the split potential method in (34),
and using Störmer-Verlet, over a range of time steps. In this case, the total
simulation duration is 1 s, and the initial condition is chosen to be large,
with $\alpha=100$. No regularisation was employed in this case (see Section
2.5). In general, the errors for Störmer-Verlet and the split scheme (34)
track each other quite closely, with the non-split method (23) performing
somewhat worse. Second order accuracy is easily observed in all cases.
Stability for Störmer-Verlet is dependent on the size of the initial
condition; for $\alpha=100$, the range of time steps over which Störmer-Verlet
is unstable is indicated in the figure.

Figure 3: Error, as defined in (49), as a function of time step $k$ for the
explicit Hamiltonian scheme (23), the split potential form (34), and Störmer-
Verlet (18), for the Fermi-Pasta-Ulam problem, where ${\bf
q}^{(0)}=[0,0,0,100,0,0]$.

###  4.2 Nonlinear String Vibration

As a second example, consider the motion of a taut string. For sufficiently
large displacements, the linear wave equation must be augmented by appropriate
nonlinear terms to account for the amplitude-dependent physical phenomena
observed during motion. One may obtain a geometrically exact model by
considering large strains, and applying Hooke’s law [33]. This model is
written compactly as

$\rho
A\partial_{t}^{2}{u}=\partial_{x}\left(\frac{\partial{\mathcal{V}}}{\partial\zeta}\right),\quad\rho
A\partial_{t}^{2}{v}=\partial_{x}\left(\frac{\partial{\mathcal{V}}}{\partial\eta}\right).$
(50)

Here, $u=u(x,t):\mathcal{D}\times\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}$
represents the transverse displacement, for a spatial domain
$\mathcal{D}=[0,L]$, where $L$ is the length of the unstretched string.
Similarly, $v=v(x,t)$ is the longitudinal displacement. In (50),
$\partial_{t}$ and $\partial_{x}$ represent partial differentiation with
respect to time $t$ and spatial coordinate $x$. Furthermore, we define
$\zeta=\partial_{x}u$, $\eta=\partial_{x}v$. The function
${\mathcal{V}}={\mathcal{V}}(\zeta,\eta):\mathbb{R}^{2}\rightarrow\mathbb{R}_{0}^{+}$
is a non-negative potential density, which may be given in two equivalent
forms as

$\displaystyle{\mathcal{V}}(\zeta,\eta)=$
$\displaystyle\frac{EA}{2}\left(\zeta^{2}+\eta^{2}\right)-(EA-
T_{0})\left(\sqrt{(1+\eta)^{2}+\zeta^{2}}-1-\eta\right),$ (51a)
$\displaystyle{\mathcal{V}}(\zeta,\eta)=$
$\displaystyle\frac{T_{0}}{2}\left(\zeta^{2}+\eta^{2}\right)+\frac{EA-
T_{0}}{2}\left(\sqrt{(1+\eta)^{2}+\zeta^{2}}-1\right)^{2}.$ (51b)

Here, the various constants that appear are: $\rho$, the volume density in
kg$\cdot$m-3; $E$, Young’s modulus, in kg$\cdot$ s-2 $\cdot$ m-1; $A$, the
area of the string cross section in m2; and $T_{0}$, the applied tension in
kg$\cdot$ m $\cdot$ s-2. It is easy to show that (51a) is the same as (51b),
and that ${\mathcal{V}}\geq 0\,\,\forall\,\,(\zeta,\,\eta)$. In (51b), one may
easily split $\mathcal{V}$ into a quadratic part, plus a non-negative
nonlinear term, provided that $EA>T_{0}$ (a condition that is generally
satisfied, for instance by all strings of interest in musical acoustics):
ultimately, this allows for numerical solution using a split potential form as
described in Section 3.4. Both forms have been used in previous works: (51a)
in e.g. [34, 11]; (51b) in [35]. System (50) is Hamiltonian, with the total
energy defined by

$H=\int_{\mathcal{D}}\frac{\rho
A}{2}\,\left((\partial_{t}u)^{2}+(\partial_{t}v)^{2}\right)+\mathcal{V}\,\mathrm{d}x.$
(52)

Energy conservation holds under a suitable set of boundary conditions. Here,
conditions of fixed type are considered, such that $u=v=0$ at $x=\\{0,L\\}$.

####  4.2.1 Semi-discrete Form

The domain $\mathcal{D}$ may be divided into segments of length $h$, the grid
spacing. Let $M$ be the total number of segments, yielding $M-1$ grid points,
not including the end points, where the solution is fixed to zero. The
continuous functions $u(x,t)$, $v(x,t)$ may then be approximated by grid
functions $u_{l}(t)$, $v_{l}(t)$, at the grid point $x_{l}=lh$, $l=1,...,M-1$.
Approximations to $\partial_{x}$ may be given as difference operators,
expressed here in terms of their action on a grid function $u_{l}(t)$:

$D_{+}u_{l}=\frac{1}{h}\left(u_{l+1}-u_{l}\right),\,\,\,D_{-}u_{l}=\frac{1}{h}\left(u_{l}-u_{l-1}\right).$
(53)

From these, one may also define the second difference operator as
$D_{2}=D_{+}D_{-}$. Furthermore, let $\zeta_{l}=D_{-}u_{l}$,
$\eta_{l}=D_{-}v_{l}$. A semi-discrete form for (50) is then obtained as

$\rho
A\ddot{u}_{l}=D_{+}\left(\frac{\partial{\mathcal{V}}_{l}}{\partial\zeta_{l}}\right),\qquad\rho
A\ddot{v}_{l}=D_{+}\left(\frac{\partial{\mathcal{V}}_{l}}{\partial\eta_{l}}\right),$
(54)

where ${\mathcal{V}}_{l}\triangleq{\mathcal{V}}(\zeta_{l},\eta_{l})$. It is
convenient to write system (54) compactly using the consolidated state vector
${\bf q}=[{\bf u}^{\intercal},{\bf v}^{\intercal}]^{\intercal}$. In this form,
the difference operators are represented by matrices, such that $D_{-}$
becomes the $M\times(M-1)$ matrix ${\bf D}_{-}$ given in terms of its action
on e.g. $\bf u$ as

${\bf D}_{-}{\bf u}=\frac{1}{h}([{\bf u}^{\intercal},0]-[0,{\bf
u}^{\intercal}]).$ (55)

Then, define ${\bf D}_{+}=-{\bf D}_{-}^{\intercal}$, and ${\bf D}_{2}={\bf
D}_{+}{\bf D}_{-}$. In vector form, (54) becomes

$\rho A\ddot{\bf q}=\frac{1}{h}\begin{bmatrix}{\bf
D}_{+}\nabla_{\bm{\zeta}}{V}\\\ {\bf
D}_{+}\nabla_{\bm{\eta}}{V}\end{bmatrix},$ (56)

where ${\bm{\zeta}}={\bf D}_{-}{\bf u}$, ${\bm{\eta}}={\bf D}_{-}{\bf v}$.
This system conserves the semi-discrete energy of the form (1), with

${\bf M}=\rho Ah{\bf I}_{2M-2},\quad V=h\sum_{l=1}^{M}{\mathcal{V}}_{l}.$ (57)

Here ${\bf I}_{2M-2}$ is the $(2M-2)\times(2M-2)$ identity matrix. A proof is
obtained immediately by noting that

$\dot{V}=\dot{\bm{\zeta}}^{\intercal}\nabla_{\bm{\zeta}}V+\dot{\bm{\eta}}^{\intercal}\nabla_{\bm{\eta}}V=-\dot{\bf
u}^{\intercal}{\bf D}_{+}\nabla_{\bm{\zeta}}V-\dot{\bf v}^{\intercal}{\bf
D}_{+}\nabla_{\bm{\eta}}V=-\dot{\bf q}^{\intercal}\begin{bmatrix}{\bf
D}_{+}\nabla_{\bm{\zeta}}{V}\\\ {\bf
D}_{+}\nabla_{\bm{\eta}}{V}\end{bmatrix}.$ (58)

Thus, left-multiplying (56) by $h\dot{\bf q}^{\intercal}$ yields (57). A split
potential form is also available, via (51b). In this case,

${\bf K}=-T_{0}h\begin{bmatrix}{\bf D}_{2}&{\bf 0}\\\ {\bf 0}&{\bf
D}_{2}\end{bmatrix},\quad V^{\prime}=\frac{h(EA-
T_{0})}{2}\sum_{l=1}^{M}\left(\sqrt{(1+\eta_{l})^{2}+\zeta_{l}^{2}}-1\right)^{2}$
(59)

####  4.2.2 Numerical Methods

System (56) may be integrated in time in a number of ways. First, introduce
the discrete time vector ${\bf q}^{n}$, approximating ${\bf q}(t)$ at the time
$t_{n}=kn$. Furthermore, let $V^{n}=V({\bm{\zeta}}^{n},{\bm{\eta}}^{n})$. An
explicit time stepping scheme is obtained by application of the Störmer-Verlet
algorithm:

${\bf q}^{n+1}=2{\bf q}^{n+1}-{\bf q}^{n-1}+\frac{k^{2}}{\rho
Ah}\begin{bmatrix}{\bf D}_{+}\nabla_{{\bm{\zeta}}^{n}}{V}^{n}\\\ {\bf
D}_{+}\nabla_{{\bm{\eta}}^{n}}{V}^{n}\end{bmatrix}.$ (60)

While simple, this scheme does not conserve a positive discrete energy, and
instabilities may occur at large displacements. A stable scheme may be
obtained by an energy-conserving discretisation of the gradient [36, 37]:

${\bf q}^{n+1}=2{\bf q}^{n+1}-{\bf q}^{n-1}+\frac{k^{2}}{\rho
A}\begin{bmatrix}{\bf D}_{+}\,\mathfrak{g}^{n}_{{\bm{\zeta}}}\\\ {\bf
D}_{+}\,\mathfrak{g}^{n}_{{\bm{\eta}}}\end{bmatrix}$ (61)

where the discrete gradients are defined as

$(\mathfrak{g}^{n}_{{\bm{\zeta}}})_{l}=\frac{\mathcal{V}({\zeta}^{n+1}_{l},{\eta}^{n}_{l})-\mathcal{V}({\zeta}_{l}^{n-1},{\eta}_{l}^{n})}{{\zeta}_{l}^{n+1}-{\zeta}^{n-1}_{l}},\quad(\mathfrak{g}^{n}_{{\bm{\eta}}})_{l}=\frac{\mathcal{V}({\zeta}^{n}_{l},{\eta}^{n+1}_{l})-\mathcal{V}({\zeta}_{l}^{n},{\eta}_{l}^{n-1})}{{\eta}_{l}^{n+1}-{\eta}^{n-1}_{l}},$
(62)

with $l=1,..,M$. This scheme leads to discrete energy conservation, and
unconditional stability, but it is fully implicit, and will generally require
the use of iterative root finding routines such as Newton-Raphson [37].

Finally, scheme (23) results from Hamiltonian (1) with definitions (57). The
split-potential form (34) is also available, via (59). In the latter case, a
stability condition arises as per (40), such that

$k\leq\sqrt{\rho A/T_{0}}\,\,h.$ (63)

This is the standard Courant-Friedrichs-Lewy stability condition for the one-
dimensional linear wave equation [38].

####  4.2.3 Numerical Results

As a first example, let the string be initialised in its first linear mode of
vibration in the transverse direction, that is

$u(x,0)=\alpha\sqrt{A}\sin(\pi x/L),\quad v(x,0)=0,$ (64)

and let the initial velocity be zero for both transverse and longitudinal
motion. The amplitude parameter $\alpha$ is nondimensional. The string
parameters are taken from [39], for the C3 piano string, and are: $\rho=7850$
kg m-3; $A=8.87\cdot 10^{-7}$ m2; $L=1.259$ m; $E=2.02\cdot 10^{11}$ kg$\cdot$
s-2 $\cdot$ m-1; $T_{0}=759$ kg$\cdot$ m $\cdot$ s-2. Figure 4 shows snapshots
of the computed solution using scheme (34), and under various initial
amplitudes $\alpha$. In this case, a regularisation of the potential energy
$V^{\prime}$ is employed (see Section 2.5), with a shift of $\epsilon=10^{8}$.
Typical amplitude-dependent phenomena are visible: the frequency of vibration
increases with the initial amplitude, and the initial shape deforms
progressively.

Figure 4: Snapshots of the geometrically exact nonlinear string, computed
using scheme (34), with split potential form as per (59). Here, $k=2.6\cdot
10^{-7}$, and the grid spacing is chosen as $h=1.05\sqrt{E/\rho}\,k$. The
string’s initial normalised amplitude $\alpha$ is as indicated.

In Figure 5, the output waveform of scheme (34) are checked against a
reference solution obtained using the Störmer-Verlet algorithm (60), for the
large input amplitude $\alpha=300$, indicating that the output of the two
schemes converges to a common solution.

Figure 5: Transverse displacement $u$ (in green) for a string initialised with
$\alpha=300$, at $x=0.5L$. In black is the reference solution, computed using
scheme (60) with a time step $k=1.04\cdot 10^{-7}$. In green are the waveforms
computed using (34), and with time steps as indicated. The grid spacing is
chosen as $h=1.05\sqrt{E/\rho}\,k$.

The relative energy error, as defined in (45), is shown in Figure 6, showing
conservation to near machine accuracy. The larger range of variation here,
compared with the case of the Fermi-Pasta-Ulam system is a result of the much
larger state size and the resulting accumulation of errors. Notice in
particular that here, in contrast to the case of the Fermi-Pasta-Ulam system,
even though the energy variation is extremely small, there is now a clear
correlation with the numerical solution. Such an effect is highly dependent on
finite wordlength effects in double precision floating point, including the
precise order of operations in the final update, and within the expression
used to calculate the relative error (45), and is not well understood.

Figure 6: Energy variation of scheme (34) for the geometrically exact
nonlinear string, using the split potential form as per (59). Here, the energy
error is as per (45). Here, the time step is $k=2.4\cdot 10^{-7}$, and the
grid spacing is chosen as $h=1.05\sqrt{E/\rho}\,k$. The string is initialised
with using $\alpha=300$.

Figure 7 presents the compute times for schemes (60), (61) and (34). The fully
implicit, conservative scheme (61) is the slowest, requiring a few iterations
of the Newton-Raphson routine per time step. The proposed scheme (34) has
compute times of the same order of magnitude as the fully explicit Störmer-
Verlet algorithm, and a few orders of magnitude smaller than the fully
implicit scheme. All simulations were run in Matlab, using a 2016 MacBook Pro
with a 2.9 GHz Quad-Core Intel Core i7 chip.

Figure 7: Computation times, in s, for the geometrically exact nonlinear
string, using the fully implicit scheme (61) (black), Störmer-Verlet (60)
(blue) and the proposed scheme (34) (green). The grid spacing for a given time
step is chosen as $h=1.05\sqrt{E/\rho}\,k$. For the fully implicit scheme,
Newton-Raphson is run with a tolerance of $10^{-13}$, and the maximum number
of iterations is capped at 20.

###  4.3 Nonlinear Plate Vibration: The Föppl-von Kármán System

As a final example, consider the problem of the transverse vibration of a thin
plate at high amplitudes. A commonly used model is the so-called dynamic
analogue of the Föppl-von Kármán equations (see, e.g., [40, 41, 42]), that
have been used extensively recently in studies of wave turbulence [43, 44,
45]. Time-stepping methods have been developed [46], including a linearly-
implicit energy-conserving method [47].

Though the dynamic Föppl-von Kármán equations can be written directly in
Hamiltonian form, they are most commonly presented as the following pair of
coupled partial differential equations:

$\rho\xi\partial_{t}^{2}q=-Q\Delta\Delta
q+{\mathcal{L}}(q,F)\qquad\qquad\frac{2}{E\xi}\Delta\Delta
F=-{\mathcal{L}}(q,q)\,.$ (65)

Here, $q(x,y,t)$ and $F(x,y,t)$ are the transverse displacement of the plate
and Airy stress function respectively; both are functions of spatial
coordinates $(x,y)\in{\mathcal{D}}\subset{\mathbb{R}}^{2}$, and time $t\geq
0$. In this simple example, the spatial domain ${\mathcal{D}}$ will be taken
to be the square of side length $L$ m, so that ${\mathcal{D}}=[0,L]^{2}$.
$\partial_{t}$ and $\Delta$ represent partial differentiation with respect to
time $t$ and the 2D Laplacian operator, respectively. $\Delta\Delta$ is thus
the biharmonic operator. For simplicity, boundary conditions are assumed to be
of simply supported type over the boundary $\partial{\mathcal{D}}$ of
${\mathcal{D}}$, so that

$q=\Delta q=0\qquad\qquad F=\Delta F=0\qquad{\rm
over}\quad\partial{\mathcal{D}}\,.$ (66)

The various constants that appear in (65) are: $\rho$, the material density,
in kg$\cdot$ m-3; $\xi$, the plate thickness, in m; $E$, Young’s modulus, in
kg$\cdot$ s-2 $\cdot$ m-1; and the flexural rigidity
$Q=E\xi^{3}/12(1-\nu^{2})$, where $\nu$ is Poisson’s ratio for the plate
material. ${\mathcal{L}}$ is a bilinear operator, defined in terms of its
action on two functions $f(x,y)$ and $g(x,y)$ as

${\mathcal{L}}(f,g)=\partial_{x}^{2}f\partial_{y}^{2}g+\partial_{y}^{2}f\partial_{x}^{2}g-2\partial_{x}\partial_{y}f\partial_{x}\partial_{y}g\,,$
(67)

where $\partial_{x}$ and $\partial_{y}$ represent partial differentiation with
respect to $x$ and $y$, respectively. Notice that only the first of (65) is
dynamic; the pair of equations (65) could be rewritten as a single equation in
displacement $q$ alone, and would constitute a second order in time cubic
nonlinear PDE.

System (65) is Hamiltonian, with the total energy defined by

$H=\iint_{{\mathcal{D}}}\frac{\rho\xi}{2}(\partial_{t}q)^{2}+\frac{Q}{2}(\Delta
q)^{2}+\frac{1}{2E\xi}(\Delta F)^{2}d\sigma\,.$ (68)

This particular form of the energy holds under fixed edge boundary conditions,
(such as simply supported (66)), and must be modified under other types of
conditions, such as free-edge. Notice that the final two terms under the
integral above, which correspond to the potential energy, separate into a
quadratic form in $q$, representing stored energy due to linear effects, and a
quadratic form in $F$, representing additional nonlinear effects; both terms
are individually non-negative, signalling that in numerical design, a
splitting of the form described in Section 2.4 is available.

####  4.3.1 Semi-discrete Form

For the square region ${\mathcal{D}}$, of side length $L$, one may start by
defining grid locations $x_{l}=lh$, $y_{m}=mh$, where $l,m=1,\ldots,M-1$, for
some integer $M$ such that $M=L/h$, where $h$ is a grid spacing. The semi-
discrete grid functions $q_{l,m}(t)$ and $F_{l,m}(t)$ thus represent
approximations to $q(x,y,t)$ and $F(x,y,t)$ at $x=x_{l}$ and $y=y_{m}$,
respectively.

Approximations to the spatial derivative operators $\partial_{x}$ and
$\partial_{y}$ may be written, in terms of their action on a grid function
$u_{l,m}(t)$ (such as $q_{l,m}$ or $F_{l,m}$ as defined above), as

$D_{x\pm}u_{l,m}=\mp\frac{1}{h}\left(u_{l,m}-u_{l\pm 1,m}\right)\qquad\qquad
D_{y\pm}u_{l,m}=\mp\frac{1}{h}\left(u_{l,m}-u_{l,m\pm 1}\right)\,.$ (69)

Under simply supported conditions, when grid points outside the range
$l,m=1,\ldots,M-1$ are referred to in the definitions above, such values are
assumed to be zero. These are the most basic forward and backward
approximations to derivatives, but are sufficient for the present
purposes—more elaborate approximations (such as those of spectral type [48])
are available.

From the basic operations defined in (69), centered approximations to the
Laplacian $\Delta$ and biharmonic operator $\Delta\Delta$ follow as

$D_{\Delta}=D_{x+}D_{x-}+D_{y+}D_{y-}\qquad\qquad
D_{\Delta\Delta}=D_{\Delta}D_{\Delta}\,.$ (70)

It is important to note that this particular construction of the biharmonic
approximation $D_{\Delta\Delta}$, through a product of Laplacian
approximations $D_{\Delta}$ under fixed conditions ensures that the simply
supported conditions (66) are satisfied.

A semi-discrete approximation to (65) then follows as

$\rho\xi\ddot{q}_{l,m}=-QD_{\Delta\Delta}q_{l,m}+\ell(q,F)=0\qquad\qquad\frac{2}{E\xi}D_{\Delta\Delta}F_{l,m}=-\ell(q,q)\,.$
(71)

The operator $\ell(\cdot,\cdot)$ approximates ${\mathcal{L}}(\cdot,\cdot)$, as
defined in (67). One useful centered approximation, operating on two grid
functions $f_{l,m}$ and $g_{l,m}$ is [47]:

$\displaystyle\ell(f,g)$ $\displaystyle=$ $\displaystyle
D_{x+}D_{x-}fD_{y+}D_{y-}g+D_{y+}D_{y-}fD_{x+}D_{x-}g$
$\displaystyle-\tfrac{1}{2}D_{x+}D_{y+}fD_{x+}D_{y+}g-\tfrac{1}{2}D_{x+}D_{y-}fD_{x+}D_{y-}g-$
$\displaystyle\tfrac{1}{2}D_{x-}D_{y+}fD_{x-}D_{y+}g-\tfrac{1}{2}D_{x-}D_{y-}fD_{x-}D_{y-}g\,.$

It is useful to represent the semi-discrete ODE system (71) in vector form,
using $(M-1)^{2}\times 1$ vectors ${\bf q}$ and ${\bf F}$. In this
representation, the operators $D_{x+}$ and $D_{y+}$ become $(M-1)^{2}\times
M(M-1)$ matrices ${\bf D}_{x+}$ and ${\bf D}_{y+}$, and ${\bf D}_{x-}=-{\bf
D}_{x+}^{\intercal}$ and ${\bf D}_{y-}=-{\bf D}_{y+}^{\intercal}$. The
operators $D_{\Delta}$ and $D_{\Delta\Delta}$ become
$(M-1)^{2}\times(M-1)^{2}$ matrices ${\bf D}_{\Delta}$ and ${\bf
D}_{\Delta\Delta}$ respectively. One arrives at the form

$\rho\xi\ddot{{\bf q}}=-Q{\bf D}_{\Delta\Delta}{\bf q}+\ell({\bf q},{\bf
F})=0\qquad\qquad\frac{2}{E\xi}{\bf D}_{\Delta\Delta}{\bf F}=-\ell({\bf
q},{\bf q})\,.$ (72)

This second order in time ODE system serves as the starting point for methods
such as Störmer-Verlet and other linearly-implicit energy-conserving methods,
as described below.

By introducing the momentum variable ${\bf p}={\bf M}\dot{{\bf q}}$, the first
order system (3) equivalent to (72) results from a Hamiltonian of the form of
(1), with $N=(M-1)^{2}$, and where

${\bf M}=\rho\xi h^{2}{\bf I}_{(M-1)^{2}}\qquad V=\frac{Qh^{2}}{2}\|{\bf
D}_{\Delta}{\bf q}\|^{2}+\frac{h^{2}}{2E\xi}\|{\bf D}_{\Delta}{\bf
F}\|^{2}\,.$ (73)

Here, ${\bf I}_{(M-1)^{2}}$ is the $(M-1)^{2}\times(M-1)^{2}$ identity matrix.
A natural splitting of the potential energy $V$ as in (14) follows, with

${\bf K}=Qh^{2}{\bf D}_{\Delta\Delta}\qquad\qquad
V^{\prime}=\frac{h^{2}}{2E\xi}\|{\bf D}_{\Delta}{\bf F}\|^{2}\,.$ (74)

####  4.3.2 Numerical Methods

Beginning from the second order system (72), one may introduce the discrete
time vectors ${\bf q}^{n}$ and ${\bf F}^{n}$, representing approximations to
$q$ and $F$ at $t=nk$, for integer $n$, and where $k$ is the time step in s.
Störmer-Verlet integration results immediately in:

${\bf q}^{n+1}=\left(2{\bf I}_{(M-1)^{2}}-\frac{Qk^{2}}{\rho\xi}{\bf
D}_{\Delta\Delta}\right){\bf q}^{n}-{\bf
q}^{n-1}+\frac{k^{2}}{\rho\xi}\ell({\bf q}^{n},{\bf
F}^{n})\qquad\qquad\frac{2}{E\xi}{\bf D}_{\Delta\Delta}{\bf F}^{n}=-\ell({\bf
q}^{n},{\bf q}^{n})\,.$ (75)

The first of these updates is explicit, but relies on ${\bf F}^{n}$, which
must be obtained from the second equation through the solution of a linear
system involving the biharmonic operator ${\bf D}_{\Delta\Delta}$.

The Störmer-Verlet scheme above is not conservative, and is prone to numerical
instability. A slight variant, however, leads to a scheme with exact energy
conservation:

$\displaystyle{\bf q}^{n+1}=\left(2{\bf
I}_{(M-1)^{2}}-\frac{Qk^{2}}{\rho\xi}{\bf D}_{\Delta\Delta}\right){\bf
q}^{n}-{\bf q}^{n-1}+\frac{k^{2}}{2\rho\xi}\ell({\bf q}^{n},{\bf F}^{n+1}+{\bf
F}^{n-1})$ (76a) $\displaystyle\frac{1}{E\xi}{\bf D}_{\Delta\Delta}\left({\bf
F}^{n+1}+{\bf F}^{n}\right)=-\ell({\bf q}^{n+1},{\bf q}^{n})\,.$ (76b)

Due to the bilinearity of the operator $\ell$, this scheme is linearly
implicit—at each time step, ${\bf q}^{n+1}$ and ${\bf F}^{n+1}$ must be solved
simultaneously, using a linear system constructed anew at each time step;
Störmer-Verlet, in contrast, requires only the solution of a linear system in
${\bf D}_{\Delta\Delta}$, which is of a known form. Thus it may be expected
that scheme (76), while energy-conserving and provably numerically stable
[47], is significantly slower to execute than the Störmer-Verlet scheme (75).

Finally, from the Hamiltonian form given in (3), with $V$ and ${\bf M}$ as
given in (73), as well as the splitting of the potential energy as in (74), a
time-interleaved scheme of the form (34) results. This scheme possesses a non-
negative exactly conserved numerical energy under the condition (40) which, in
this case, reduces to a lower bound on the grid spacing $h$ in terms of the
time step $k$:

$h\geq h_{{\rm min}}=2\sqrt{k}\left(D/\rho\xi\right)^{\tfrac{1}{4}}\,.$ (77)

The scheme, like Störmer-Verlet, relies on a linear system solution involving
the biharmonic operator ${\bf D}_{\Delta\Delta}$, but is otherwise explicit.

####  4.3.3 Numerical Examples

As an example, consider initialisation of the Föppl-von Kármán system (65)
through its lowest linear mode shape:

$q(x,y,0)=\alpha\xi\sin(\pi x/L)\sin(\pi
y/L)\qquad\partial_{t}q|_{x,y,t=0}=0\,,$ (78)

where the dimensionless parameter $\alpha$ controls the maximum amplitude of
the initial condition relative to the plate thickness $\xi$. Furthermore, the
plate is assumed to be made of steel with $E=2\times 10^{11}$ Pa, $\rho=7850$
kg$\cdot$ m-3, and $\nu=0.3$, and to be of thickness $\xi=2$ mm and side
length $L=0.5$ m.

Using the scheme (34), with the grid spacing and time step chosen according to
(77), typical amplitude-dependent behaviour is observed. See Figure 8. At a
low initial condition amplitude of $\alpha=0.01$, behaviour is essentially
linear. At amplitudes near the plate thickness at $\alpha=2$, the period of
oscillation decreases, and spontaneous mode generation is observed, and for
large amplitudes, such as $\alpha=10$, turbulent behaviour is observed. Under
such stringent high-amplitude conditions ($\alpha=10$), scheme (34) remains
stable, and results converge to those of the reference solution, computed
using Störmer-Verlet with a time step of $k=2.5\times 10^{-6}$ s. No
regularisation (see Section 2.5) is used in this case. See Figure 9. Under
these conditions, Störmer-Verlet is unstable for $k>2\times 10^{-5}$ s.

Figure 8: Snapshots of the time evolution of a Föppl-von Kármán plate, at
times as indicated, and for different initial condition amplitudes
$\alpha=0.01$, $\alpha=2$ and $\alpha=10$. Results are computed using scheme
(34), with a time step of $k=5\times 10^{-6}$.  Figure 9: The output waveform
$q_{{\rm o}}(t)$ drawn from a plate vibrating at high amplitude, with
$\alpha=10$. The reference solution, generated using Störmer-Verlet with
$k=2.5\times 10^{-6}$ s is shown in black; results are shown for scheme (34),
in green, at different time steps $k$, as indicated.

The simulation of the Föppl-von Kármán system is computationally intensive.
Most interesting in this case is a comparison of computation times between
Störmer-Verlet (75), the linearly-implicit energy-conserving method (76), and
the explicit energy-conserving scheme (34). Unavoidable in all cases is some
form of linear system solution; for (75) and (34), this will involve the
biharmonic operator, which is constant over the course of a simulation, and
thus amenable to factorisation techniques (e.g. Cholesky) to decrease solution
times. For the linearly-implicit exact energy conserving method (76), however,
this is not the case—the linear system to be solved must be constructed anew
at each time step. This is reflected in timings, as shown in Figure 10, for
simulations of 1 second for different choices of time step $k$. Computation
was performed in Matlab on a Lenovo P50 with an Intel Xeon E3 v5. As can be
seen, computation time for the explicit scheme (34) is far lower than for
scheme (76), and very nearly on par with Störmer-Verlet.

Figure 10: Computation times, in s, for the Föppl-von Kármán plate, for a 1 s
simulation duration, using the linearly implicit exact energy conserving
method (76), the explicit exact energy conserving method (34), and Störmer-
Verlet (75).

Finally, see Figure 11, illustrating the relative energy variation for scheme
(34) for the Föppl-von Kármán plate; as in the case of nonlinear string
vibration, the energy variation is of the order of $10^{-15}$, with some
correlation with the numerical solution visible—see Figure 6 for comparison.

Figure 11: Relative numerical energy variation, as defined in (45), for scheme
(34) for the Föppl-von Kármán plate, with a high initial condition of the form
of (78), with $\alpha=10$. The time step is $k=10^{-5}$ s.

##  5 Concluding Remarks

The design of exactly energy-conserving methods for Hamiltonian systems has
progressed from fully implicit designs through, more recently, to explicit
methods for which exact energy conservation can be attained through an
approximation to a continuous integral of the potential energy, or, more
importantly, to linearly implicit designs based on invariant energy
quadratisation. The main novelty in this paper is to call attention to
structure within such linearly implicit designs—structure that can be
exploited in order arrive at fully explicit methods. These exactly
conservative methods are of roughly the same computational cost as the most
efficient non-conservative explicit methods—with the additional feature of a
clear means of ensuring numerical stability, either unconditionally, or, if a
splitting of the potential energy is employed, under well-defined conditions
on the time step that are independent of the initial conditions. Accuracy is
of second order, and is plainly evident in the centered (but interleaved)
discretisation approach, and borne out by simulation results (see Section
4.1). It is not clear whether it is possible to extend this framework to
obtain higher order accuracy.

These methods are not completely general, and require, additionally, a
condition of non-negativity on the potential energy, as per invariant energy
quadratisation approaches. More generally, given that the dynamics of a system
are independent of shifts in the potential energy by a constant (i.e., a
gauge), a more general condition is that the potential energy is bounded from
below. The useful technique of splitting of the potential energy introduced
here is a further restriction. But the restriction (6) mentioned above to non-
negative expressions for potential energy $V$ (or bounded from below) is
slightly more strict than necessary. More general is a restriction to
expressions $V({\bf q})$ that are single-signed for all ${\bf q}$—and even
more generally, bounded either from above or below. An important example here
is the $N$-body problem, under a gravitational potential, for which $V({\bf
q})\leq 0$. In this case, one may set, instead of (9),
$V=-\tfrac{1}{2}\psi^{2}$, and the main development follows as above, with
this sign change, and an explicit exactly energy conserving method follows as
before. In this case, however, global bounds on solution size are not
available, as the total energy itself is no longer necessarily non-negative.
See the comments in [26] regarding general dynamic stability for Hamiltonian
systems.

Only briefly alluded to here, in Section 3.5, is the extension to the case of
variable time steps—useful in modeling systems with slow/fast dynamics, as
discussed in [11]. The behaviour of such schemes remains unexplored. Another
open question is the need for regularization, as introduced in Section 2.5. In
two of the three cases presented here, the Fermi-Pasta-Ulam problem, and the
Föppl-von Kármán system, good results were obtained without the use of such
regularisation. In the case of string vibration, however, such regularisation
was necessary. It would be of great interest to know the origin of this
distinction. In all cases, however, regularisation has no impact on the
explicit nature of the schemes proposed here, and only negligible impact on
computational cost.

## Acknowledgments

M. Ducceschi was supported by the European Research Council (ERC), under grant
2020-StG-950084-NEMUS. For the purpose of open access, the first author has
applied a creative commons attribution (CC BY) licence to any author accepted
manuscript version arising.

## References

  * [1] R. LaBudde and D. Greenspan.  Energy and momentum conserving methods of arbitrary order for the numerical integration of equations of motion I. Motion of a single particle.  Numerical Mathematics, 25:323–346, 1976.
  * [2] A. Marciniak.  Energy conserving, arbitrary order numerical solutions of the N-body problem.  Numerische Mathematik, 45:207–218, 1984.
  * [3] W. Strauss and L. Vazquez.  Numerical solution of a nonlinear Klein-Gordon equation.  Journal of Computational Physics, 28:271–278, 1978.
  * [4] D. Greenspan.  Conservative numerical methods for $\ddot{x}=f(x)$.  Journal of Computational Physics, 56:28–41, 1984.
  * [5] J. Simo, N. Tarnow, and K. Wong.  Exact energy-momentum conserving algorithms for symplectic schemes for nonlinear dynamics.  Computer Methods in Applied Mechanics and Engineering, 100:63–116, 1992.
  * [6] T. Hughes, T. Caughey, and W. Liu.  Finite-Element methods for nonlinear elastodynamics which conserve energy.  Journal of Applied Mechanics, 45:366–370, 1978.
  * [7] G. Zhong and J. Marsden.  Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators.  Physics Letters A, 133(3):134–139, 1988.
  * [8] E. Hairer.  Long-time energy conservation of numerical integrators.  In L. Pardo, A. Pinkus, E. Suli, and M. Todd, editors, Foundations of Computational Mathematics: Santander 2005, pages 162–180. Cambridge University Press, Cambridge, UK, 2006.
  * [9] R. McLachlan and M. Perlmutter.  Energy drift in reversible time integration.  Journal of Physics A: Mathematical and General, 37:L593–L598, 2004\.
  * [10] E. Faou, E. Hairer, and T. Pham.  Energy conservation with non-symplectic methods: Examples and counter-examples.  BIT, 44:699–709, 2004.
  * [11] F. Marazzato, A. Ern, C. Mariotti, and L. Monasse.  An explicit pseudo-energy conserving time-integration scheme for Hamiltonian dynamics.  Computer Methods in Applied Mechanics and Engineering, 347:906–927, 2019.
  * [12] L. Vu-Quoc and S. Li.  Invariant-conserving finite difference algorithms for the nonlinear Klein-Gordon equation.  Computer Methods in Applied Mechanics and Engineering, 107:341–391, 1993.
  * [13] J. Diaz and M. Grote.  Energy conserving explicit local time stepping for second-order wave equations.  SIAM Journal of Scientific Computing, 31(3):1985–2014, 2009.
  * [14] Y. Chin and C. Qin.  Explicit energy-conserving schemes for the three-body problem,.  Journal of Computational Physics, 83:485–493, 1989.
  * [15] B. Shadwick, J. Bowman, and P. Morrison.  Exactly conservative integrators.  SIAM Journal of Applied Mathematics, 59(3):1112–1133, 1998.
  * [16] G. Quispel and D. McLaren.  A new class of energy-preserving numerical integration methods.  Journal of Physics A: Mathematical and Theoretical, 41(4):1–7, 2008.
  * [17] L. Brugnano, F. Iavernaro, and D. Trigiante.  A two-step, fourth-order method with energy preserving properties.  Computer Physics Communications, 183:1860–1868, 2012.
  * [18] X. Yang and D. Han.  Linearly first- and second-order, unconditionally energy stable schemes for the phase field crystal model.  Journal of Computational Physics, 330:1116–1134, 2016.
  * [19] X. Yang, J. Zhao, and Q. Wang.  Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method.  Journal of Computational Physics, 333:104–127, 2017.
  * [20] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  Journal of Computational Physics, 353:407–416, 2018.
  * [21] Y. Gong and J. Zhao.  Energy-stable Runge–Kutta schemes for gradient flow models using the energy quadratization approach.  Applied Mathematics Letters, 94:224–231, 2019.
  * [22] H. Zhang, X. Qian, and S. Song.  Novel high-order energy-preserving diagonally implicit Runge–Kutta schemes for nonlinear Hamiltonian ODEs.  Applied Mathematics Letters, 102:1–9, 2020.
  * [23] S. Sato, Y. Miyatake, and J. Butcher.  High-order linearly implicit schemes conserving quadratic invariants, 2021\.  Available at: arXiv:2203.00944v1.
  * [24] J. Sherman and W. Morrison.  Adjustment of an inverse matrix corresponding to a change in one element of a given matrix.  Annals of Mathematical Statistics, 21(1):124–127, 1950.
  * [25] S. Bilbao and M. Ducceschi.  Fast explicit algorithms for Hamiltonian numerical integration.  In Proceedings of the European Nonlinear Dynamics Conference, Lyon, France, July 2022.
  * [26] O. Gonzalez and J. Simo.  On the stability of symplectic and energy-momentum algorithms for non-linear Hamiltonian systems with symmetry.  Computer Methods in Applied Mechanics and Engineering, 134(4):197–222, 1996.
  * [27] W. Terrell.  Stability and Stabilization.  Princeton University Press, Princeton, New Jersey, 2009.
  * [28] Z. Liu and X. Li.  The exponential scalar auxiliary variable (e-sav) approach for phase field models and its explicit computing.  SIAM Journal on Scientific Computing, 42(3):B630–B655, 2020.
  * [29] L. Lin, Z. Yang, and S. Dong.  Numerical approximation of incompressible navier-stokes equations based on an auxiliary energy variable.  Journal of Computational Physics, 388:1–22, 2019.
  * [30] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration illustrated by the Störmer–Verlet method.  Acta Numerica, 12:399–450, 2003.
  * [31] E. Fermi, J. Pasta, and S. Ulam.  Studies of nonlinear problems.  Technical Report Los Alamos LA-1940, 1955.
  * [32] E. Hairer, C. Lubich, and G. Wanner.  Geometric Numerical Integration.  Springer Series in Computational Mathematics, 2006.
  * [33] P. Morse and U. Ingard.  Theoretical Acoustics.  Princeton University Press, Princeton, NJ, USA, 1968.
  * [34] J. Chabassier and P. Joly.  Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string.  Computer Methods in Applied Mechanics and Engineering, 199(45):2779–2795, 2010.
  * [35] M. Ducceschi and S. Bilbao.  Simulation of the geometrically exact nonlinear string via energy quadratisation.  Journal of Sound and Vibration, 534:117021, 2022.
  * [36] T. Itoh and K. Abe.  Hamiltonian-conserving discrete canonical equations based on variational difference quotients.  Journal of Computational Physics, 76(1):85–102, 1988.
  * [37] M. Ducceschi and S. Bilbao.  Non-iterative, conservative schemes for geometrically exact nonlinear string vibration.  In Proceedings of the International Conference on Acoustics (ICA 2019), Aachen, Germany, September 2019.
  * [38] R. Courant, K. Friedrichs, and H. Lewy.  On the partial differential equations of mathematical physics.  Mathematische Annalen, 100:32–74, 1928.  In German.
  * [39] J. Chabassier and M. Duruflé.  Physical parameters for piano modeling.  Technical report, 2012.  Available at https://hal.inria.fr/hal-00688679v1/document.
  * [40] T. von Kármán.  Festigkeitsprobleme im maschinenbau.  Encyklopädie der Mathematischen Wissenschaften, 4(4):311–385, 1910.
  * [41] A. Föppl.  Vorlesungen über technische Mechanik.  Druck und Verlag von B.G. Teubner, Leipzig, 1907.
  * [42] A. Nayfeh and D. Mook.  Nonlinear Oscillations.  John Wiley and Sons, New York, New York, 1979.
  * [43] M. Ducceschi, O. Cadot, C. Touzé, and S. Bilbao.  Dynamics of the wave turbulence spectrum in vibrating plates: A numerical investigation using a conservative finite difference scheme.  Physica D: Nonlinear Phenomena, 280-281:73–85, 2014.
  * [44] N. Yokoyama and M. Takaoka.  Weak and strong wave turbulence spectra for elastic thin plate.  Physical Review Letters, 110(10):105501, 2013.
  * [45] G. Düring, C. Josserand, and S. Rica.  Wave turbulence theory of elastic plates.  Physica D: Nonlinear Phenomena, 347:42–73, 2017.
  * [46] R. Kirby and Z. Yosibash.  Solution of von Kármán dynamic non-linear plate equations using a pseudo-spectral method.  Computer Methods in Applied Mechanics and Engineering, 193(6–8):575–599, 2004.
  * [47] S. Bilbao.  A family of conservative finite difference schemes for the dynamical von Kármán plate equations.  Numerical Methods for Partial Differential Equations, 24(1):193–216, 2008.
  * [48] R. Vichnevetsky and J. Bowles.  Fourier Analysis of Numerical Approximations of Hyperbolic Equations.  SIAM, Philadelphia, Pennsylvania, 1982.


$W_{p}(\mathcal{G},\Omega)=W_{p}(\mathtt{prune}(\mathcal{G},\Omega,d,p))\cap\mathcal{S}.$

###### Proof

Let
$\mathtt{prune}(\mathcal{G},\Omega,d,p)=(\mathcal{G}^{\prime},\Omega^{\prime})$
with semantics
$\llbracket\mathcal{G}^{\prime}\rrbracket=(\mathcal{S}^{\prime},\mathcal{S}_{\mathit{Env}}^{\prime},\mathcal{S}_{\mathit{Sys}}^{\prime},\rho^{\prime})$.
Furthermore, let $W=W_{p}(\mathcal{G},\Omega)$ and
$W^{\prime}=W_{p}(\mathcal{G}^{\prime},\Omega^{\prime})$. As
$\llbracket\mathcal{G}\rrbracket$ (resp.
$\llbracket\mathcal{G}^{\prime}\rrbracket$) have no dead-end, a play $\pi$ is
winning for Player $\mathit{Sys}$ iff $\pi\in\Omega$ (resp.
$\pi\in\Omega^{\prime}$). Hence, let us denote the objectives for
$\mathit{Sys}$ and $\mathit{Env}$ as $\Omega_{\mathit{Sys}}=\Omega$ and
$\Omega_{\mathit{Env}}=\mathcal{S}^{\omega}\setminus\Omega$, respectively in
$\mathcal{G}$. $\Omega_{\mathit{Sys}}^{\prime}$ and
$\Omega_{\mathit{Env}}^{\prime}$ are defined analogously. Now, we will show
both direction of $W=W^{\prime}\cap\mathcal{S}$ separately.

$(\subseteq)$ Let $s_{0}\in W$, then there exists a strategy $\sigma_{p}$ for
Player $p$ in $\llbracket\mathcal{G}\rrbracket$ such that
$\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})\subseteq\Omega_{p}$.
To show that $s_{0}\in W^{\prime}$, we will define a strategy
$\sigma_{p}^{\prime}$ for Player $p$ in
$\llbracket\mathcal{G}^{\prime}\rrbracket$ such that
$\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})\subseteq\Omega_{p}^{\prime}$.

First observe that for every $s=((l,\mathbf{x}),\mathbf{i})\in
W\cap\mathcal{S}_{\mathit{Sys}}$, we have $s\not\in\llbracket d\rrbracket$ and
hence, $\mathbf{x}\uplus\mathbf{i}\models_{T}\neg d(l)$. So, it holds that

$\displaystyle(((l,\mathbf{x}),\mathbf{i}),(l^{\prime},\mathbf{x}^{\prime}))\in\rho$
$\displaystyle\iff\exists(g,u)\in\mathit{Labels}(l,l^{\prime}).\leavevmode\nobreak\
\mathbf{x}\uplus\mathbf{i}\models_{T}g,\mathbf{x}^{\prime}(x)=\chi_{\mathbf{x}\uplus\mathbf{i}}(u(x))$
$\displaystyle\iff\exists(g,u)\in\mathit{Labels}(l,l^{\prime}).\leavevmode\nobreak\
\mathbf{x}\uplus\mathbf{i}\models_{T}g\wedge\neg
d(l),\mathbf{x}^{\prime}(x)=\chi_{\mathbf{x}\uplus\mathbf{i}}(u(x))$
$\displaystyle\iff(((l,\mathbf{x}),\mathbf{i}),(l^{\prime},\mathbf{x}^{\prime}))\in\rho^{\prime}.$

Furthermore, by definition, for every $s\in W\cap\mathcal{S}_{\mathit{Env}}$,
$(s,(s,\mathbf{i}))\in\rho$ iff $(s,(s,\mathbf{i}))\in\rho^{\prime}$.
Therefore,
$\rho\cap(W\times\mathcal{S})=\rho^{\prime}\cap(W\times\mathcal{S})$.

With this, we define
$\sigma_{p}^{\prime}:(\mathcal{S}^{\prime})^{*}\cdot\mathcal{S}_{p}^{\prime}\to\mathcal{S}$
as follows. Let
$\tau\in(\mathcal{S}^{\prime})^{*}\cdot\mathcal{S}_{p}^{\prime}$. If $\tau\in
W^{*}$, then we define $\sigma_{p}^{\prime}(\tau)=\sigma_{p}(\tau)$, else we
set $\sigma_{p}^{\prime}(\tau)$ to an arbitrary successor of
$\mathit{last}(\tau)$.

Now, let
$\pi\in\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})$.
By definition, $\pi$ is an infinite play. If $\pi\in W^{\omega}$, then, by
definition,
$\pi\in\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})$ and
hence, $\pi\in\Omega_{p}\subseteq\Omega_{p}^{\prime}$.

Now, assume $\pi\not\in W^{\omega}$, then there exists $k\in\mathbb{N}$ s.t.
$\pi[k]\not\in W$ and $\pi[j]\in W$ for all $j<k$. As
$\rho\cap(W\times\mathcal{S})=\rho^{\prime}\cap(W\times\mathcal{S})$,
$\pi[0,k]$ is also a prefix of a play in
$\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})$. Since
$\pi[k]\not\in W$, Player $1-p$ has a strategy $\sigma_{1-p}$ such that
$\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(\pi[k],\sigma_{1-p})\subseteq\Omega_{1-p}$.
For safety objective or any prefix-independent objective $\Omega_{p}$, if
Player $p$ uses a strategy that is consistent with $\pi$ until $\pi[k]$ and
then switches to $\sigma_{1-p}$, then this gives us an infinite play
$\pi^{\prime}\in\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})\cap\Omega_{1-p}$.
This is a contradiction to the assumption that
$\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})\subseteq\Omega_{p}$.
Similarly, for reachability objective $\Omega_{p}=\mathit{Reach}(R)$, if
$\pi[j]\not\in R$ for all $j\leq k$, then the same argument as above gives us
a contradiction. If there exists $j\leq k$ s.t. $\pi[j]\in R$, then
$\pi\in\Omega_{p}^{\prime}$ by definition.

$(\supseteq)$ Let $s_{0}\in W^{\prime}\cap\mathcal{S}$, then there exists a
strategy $\sigma_{p}^{\prime}$ for Player $p$ in
$\llbracket\mathcal{G}^{\prime}\rrbracket$ such that
$\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})\subseteq\Omega_{p}^{\prime}$.
To show that $s_{0}\in W$, we will define a strategy $\sigma_{p}$ for Player
$p$ in $\llbracket\mathcal{G}\rrbracket$ such that
$\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})\subseteq\Omega_{p}$.

We define $\sigma_{p}:\mathcal{S}^{*}\cdot\mathcal{S}_{p}\to\mathcal{S}$ as
follows. Let $\tau\in\mathcal{S}^{*}\cdot\mathcal{S}_{p}$. If
$\tau\in\mathcal{S}^{*}$ and $\sigma_{p}^{\prime}(\tau)\in\mathcal{S}$, then
we define $\sigma_{p}(\tau)=\sigma_{p}^{\prime}(\tau)$, else we set
$\sigma_{p}(\tau)$ to an arbitrary successor of $\mathit{last}(\tau)$.

Now, let
$\pi\in\mathit{Plays}_{\llbracket\mathcal{G}\rrbracket}(s_{0},\sigma_{p})$. If
$\pi\in\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})$,
then $\pi\in\Omega_{p}^{\prime}\cap\mathcal{S}^{\omega}=\Omega_{p}$.
Otherwise, there exists prefix $\tau$ of $\pi$ s.t.
$\sigma_{p}^{\prime}(\tau)\in\\{\mathsf{sink}_{p},\mathsf{sink}_{1-p}\\}$ and
$\tau\sigma_{p}^{\prime}(\tau)$ is a prefix of a play in
$\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})$.
By construction, there is no transition from
$\mathcal{S}\cup\\{\mathsf{sink}_{p}\\}$ to $\mathsf{sink}_{1-p}$, hence,
$\sigma_{p}^{\prime}(\tau)=\mathsf{sink}_{p}$. Furthermore, by construction,
there is no play in $\Omega_{p}^{\prime}$ that visits state
$\mathsf{sink}_{p}$. This leads to a contradiction to the assumption that
$\mathit{Plays}_{\llbracket\mathcal{G}^{\prime}\rrbracket}(s_{0},\sigma_{p}^{\prime})\subseteq\Omega_{p}^{\prime}$.∎

The next statement follows directly from Lemma 5 and Lemma 4 and allows us to
soundly prune the states determined to be losing for the over-approximated
player $\mathit{OverapproxP}(\widehat{G}^{\circ})$ in an abstract game
$\widehat{G}^{\circ}$.

###### Corollary 2

Let $\widehat{G}^{\circ}$ for some $\circ\in\\{\uparrow,\downarrow\\}$ be a
$(\mathcal{P}_{\mathbb{X}},\mathcal{P}_{\mathbb{X}\cup\mathbb{I}})$-induced
abstraction of $\mathcal{G}$ and let
$p=\mathit{OverapproxP}(\widehat{G}^{\circ})$. If $d\in\mathcal{D}$ is such
that $\llbracket
d\rrbracket\subseteq\gamma(W_{1-p}(\widehat{G}^{\circ},\widehat{\Omega}))$,
then
$W_{p}(\mathcal{G},\Omega)=W_{p}(\mathtt{prune}(\mathcal{G},\Omega,d,p))\cap\mathcal{S}$.

Algorithm 5 shows the version RPGPruneCacheSolve of Algorithm 4 extended with
pruning. SolveAbstractWR works in the same way as SolveAbstract, but
additionally returns the winning region for the opponent. The correctness of
Algorithm 5 follows from Corollary 2 and Theorem 5.1.

######  Theorem 0.C.1 (Correctness of Algorithm 5)

Given a reactive program game structure $\mathcal{G}$ and a location-based
objective $\Omega$, for any $b\in\mathbb{N}$, if RPGPruneCacheSolve
terminates, then it returns
$W_{\mathit{Sys}}(\llbracket\mathcal{G}\rrbracket,\Omega).$

1 function  _RPGPruneCacheSolve(
$\mathcal{G}=(T,\mathbb{I},\mathbb{X},L,\mathit{Inv},\delta)$, $\Omega$,
$b\in\mathbb{N}$)_

2
$(\mathcal{P}_{\mathbb{X}},\mathcal{P}_{\mathbb{X}\cup\mathbb{I}}):=\mathsf{AbstractDomain}(\mathcal{G})$

3
$(\widehat{G}^{\uparrow},\widehat{G}^{\downarrow}):=\mathsf{AbstractRPG}(\mathcal{G},(\mathcal{P}_{\mathbb{X}},\mathcal{P}_{\mathbb{X}\cup\mathbb{I}}))$

4
$(U_{\mathit{Sys}},D_{\mathit{Sys}},\mathcal{H}_{\mathit{Sys}},\widehat{W}_{\mathit{Env}}):=\textsc{SolveAbstractWR}(\widehat{G}^{\uparrow},\Omega)$

5
$(U_{\mathit{Env}},D_{\mathit{Env}},\mathcal{H}_{\mathit{Env}},\widehat{W}_{\mathit{Sys}}):=\textsc{SolveAbstractWR}(\widehat{G}^{\downarrow},\Omega)$

6
$(\mathcal{G}^{\prime},\Omega^{\prime}):=\mathtt{prune}(\mathcal{G},\Omega,\widehat{W}_{\mathit{Env}},\mathit{Sys})$

7
$(\mathcal{G}^{\prime\prime},\Omega^{\prime\prime}):=\mathtt{prune}(\mathcal{G}^{\prime},\Omega^{\prime},\widehat{W}_{\mathit{Sys}},\mathit{Env})$

8
$C_{\mathit{Sys}}:=\textsc{GenerateCache}(\mathcal{G}^{\prime\prime},\widehat{G}^{\uparrow},\mathit{Sys},(U_{\mathit{Sys}},D_{\mathit{Sys}},\mathcal{H}_{\mathit{Sys}}),b)$

9
$C_{\mathit{Env}}:=\textsc{GenerateCache}(\mathcal{G}^{\prime\prime},\widehat{G}^{\downarrow},\mathit{Env},(U_{\mathit{Env}},D_{\mathit{Env}},\mathcal{H}_{\mathit{Env}}),b)$

10 return RPGSolveWithCache($\mathcal{G}^{\prime\prime},C_{\mathit{Sys}}\cup
C_{\mathit{Env}}$)

Algorithm 5 Procedure for solving reactive program games enhanced with
abstraction-based pruning and abstract template-based caching.

##  Appendix 0.D Benchmarks

All our benchmarks are modeled as reactive program games with Büchi objectives
for the system player. We describe all our benchmarks in detail below.

##### Scheduler.

This benchmark outlines two primary tasks for the system: a global task and
one repeated task. Initially, three program variables, namely $taskG$,
$taskR$, and $boundG$, are set to specific values determined by the
environment, where the system can choose to set $taskR$ as the negation of the
value specified by the environment. The system undergoes four distinct stages.
At each stage, if $taskR>0$, it proceeds to the next stage; otherwise, it has
the option to increment either of the variables $taskG$ or $taskR$. For the
repeated task, the system must attain the state $taskR>0$ in every stage,
whereas for the global task, the system must eventually ensure $taskG>boundG$.
Crucially, the system must strategically ensure that there is at least one
stage where $taskR\leq 0$. This specific condition allows the system to
increment the variable $taskG$, thereby progressing toward fulfilling the
global task. The entire process resets after the completion of the four
stages.

##### Item Processing.

In this benchmark, the system processes a tray of items, manipulating
variables such as $trayItems$ and $numItems\\_s$ through different locations.
Similar to the scheduler benchmark, these variables are initially set to the
values determined by the environment. At each location, the system is
presented with various operations it can perform on these variables. The
primary goal of the system is to guarantee the condition $numItems\\_s\geq
trayItems$ in the $done$ location. To meet this objective, the system must
strategically select the appropriate options at each step.

##### Chain Benchmarks.

We have two categories of chain benchmarks, both designed with the objective
$\mathit{Buchi}(\\{goal\\})$.

  * •

In chain $k$, we parameterize the number of locations and program variables
based on the value $k$. Specifically, the benchmark has $k+3$ locations
$\\{int,goal,sink,l_{1},l_{2},\ldots,l_{k}\\}$, $k+2$ program variables
$\\{y,c,x_{1},x_{2},\ldots,x_{k}\\}$, and one input variable $i$. Starting
from initial location $int$, it sets $c=0$ and goes to location $goal$. From
$goal$, if $c>0$, it goes to the $sink$ which only allows a self-loop. If
$c\leq 0$, it sets $x_{1}=i$, $y=i$ and goes to location $l_{1}$. From every
location $l_{j}$ with $j<k$, if $x_{j}=0$, then it can go to next location
$l_{j+1}$ by setting $x_{j+1}=y$. Alternatively, if $x_{j}\neq 0$, it can
either loop on $l_{j}$ while incrementing/decrementing $x_{j}$ by $1$, or go
to $goal$ by setting $c=i$. The transitions from $l_{k}$ are similar with the
exception that it directly goes to $goal$ if $x_{k}=0$.

  * •

In chain simple $k$, we parametrize only the number of locations based on the
value $k$. Similar to chain k, this benchmark has $k+3$ locations
$\\{int,goal,sink,l_{1},l_{2},\ldots,l_{k}\\}$ and one input varibale $i$.
However, it only has $3$ program variables $\\{y,c,x\\}$. The transition
relation is similar to that of chain k, but the program variable $x$ is used
instead of $x_{1},x_{2},\ldots,x_{k}$.

##### Robot Benchmarks.

These benchmarks describe a robot’s movement along a one-dimensional discrete
grid. The robot’s position is determined by a program variable, denoted as
$\mathit{trackPos}$ (for clarity, we use $\mathit{pos}$ in Figure 2). The
robot can move either one step forward, incrementing the variable
$\mathit{trackPos}$ by $1$, or one step backward, decrementing
$\mathit{trackPos}$ by $1$.

  * •

In robot analyze samples, the robot starts from a designated location
$\mathit{base}$. It then collects a number of samples (as specified by the
environment) from another location. Subsequently, it moves to another location
$\mathit{lab}$, conducts tests on all collected samples, and finally returns
to the $\mathit{base}$.

  * •

In robot repair, the robot starts from a designated location $\mathit{base}$.
If the device is deemed faulty (as determined by the environment) the robot
gets a number of spare parts (as specified by the environment) from the stock.
It proceeds to repair the device using all the spare parts and eventually
returns to the $\mathit{base}$.

  * •

In robot collect samples, the robot starts from a designated location
$\mathit{base}$, collects a number of samples (as specified by the
environment) from another location, and then returns to the $\mathit{base}$.

  * •

In robot deliver products $k$, the robot is tasked with purchasing $k$
distinct products for the $\mathit{office}$. It follows an alternating
pattern: choosing a product, storing the required quantity (as specified by
the environment), proceeding to the $\mathit{bank}$ to withdraw the
corresponding amount of money, purchasing the specified products from the
$\mathit{store}$, and finally returning to the $\mathit{office}$.

##### Smart Home Benchmarks.

These benchmarks describe the dynamic adjustments made by a smart home in
response to environmental inputs, involving temperature regulation, blinds
adjustment, and lighting control. The smart home turns on/off lights by
toggling a boolean program variable $\mathit{light}$ to true/false. The
blinds’ position is adjusted by incrementing or decrementing the program
variable $\mathit{blinds}$ by $0.1$. Additionally, the temperature level is
adjusted by incrementing or decrementing the program variable
$\mathit{temperature}$ by $1.0$, influenced by an input variable
$\mathit{disturbance}$.

  * •

In smart home day not empty, during the day and when the home is occupied (as
determined by the environment), the smart home activates $\mathit{daymode}$,
turns on the lights, raises the blinds, and raises the temperature to a range
specified by the program variable $\mathit{minimum}$.

  * •

In smart home day warm, the smart home also performs the tasks outlined in
smart home day not empty. Additionally, if the environment indicates the home
being too warm, it decrements the variable $\mathit{minimum}$ by $2.0$ and
continues its routine.

  * •

In smart home day cold, similar to smart home day warm, if the environment
indicates the home being too cold, it increments the variable
$\mathit{minimum}$ by $2.0$ and continues its routine.

  * •

In smart home day warm or cold, the smart home combines the tasks from both
smart home day warm and smart home day cold.

  * •

In smart home empty, during the day and when the home is empty (as determined
by the environment), the smart home activates $\mathit{daymode}$, turns off
the lights, lowers the blinds, and lowers the temperature to a range specified
by the program variable $\mathit{maximum}$.

  * •

In smart home night sleeping, during the night, the smart home deactivates
$\mathit{daymode}$ and lowers blinds. If the owner is sleeping (as determined
by the environment), it turns off lights; otherwise, it turns them on.
Subsequently, it raises the temperature to a range specified by the variable
$\mathit{minimum}$.

  * •

In smart home night empty, during the night, the smart home deactivates
$\mathit{daymode}$ and lowers blinds. If the home is empty (as determined by
the environment), it turns off the lights and lowers the temperature to a
range specified by the program variable $\mathit{maximum}$.

  * •

In smart home nightmode, the smart home combines the tasks from both smart
home night sleeping and smart home night empty.