Explicit Exactly Energy-conserving Methods for Hamiltonian Systems

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 Stormer-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.


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 "pseudoenergy" (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 × 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].

Hamiltonian Systems
Consider a system of particles, with N generalised coordinates q = [q 1 , . . . , q N ] and momenta p = [p 1 , . . . , p N ] . Here, indicates a transposition operation, so p and q are N × 1 column vectors. Both are functions of time t ≥ 0, so p = p(t) and q = q(t). Suppose also that an associated Hamiltonian H(p, q) is defined by for some constant symmetric positive definite N × N matrix M referred to as the mass matrix. V (q) is the potential energy for the system of particles. H and V are scalar functions of time t, through their dependence on p and q.
The system dynamics follow from Hamilton's equations: where dots indicate differentiation with respect to time t, and ∇ p and ∇ q are gradients with respect to p and q, respectively. For the particular form of Hamiltonian given in (1), Hamilton's equations becomeq Through time differentiation of the first of (3), and substitution of the second, a second order system of ordinary differential equations in q results: 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 q and p, so q(0) = q (0) and p(0) = p (0) , for given N vectors q (0) and p (0) . For the second order form (4), one may setq(0) = M −1 p (0) .

Comments
The system and Hamiltonian in (1) are not the most general possible. A few remarks are offered here: • The mass matrix 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 M will have an impact on computational performance.
• The Hamiltonian in (1) is separable, so that H(p, q) = T (p) + V (q), for kinetic energy T and potential energy V . In some cases of interest [8], the mass matrix M may be of the form M = M(q), disturbing this splitting, but these cases will not be considered here.
• Following from the points above, the Hamiltonian is quadratic in p, meaning that the associated dynamical system is linear in p.
• The individual displacements q i and momenta p i , i = 1, . . . , 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.

Non-negativity of the Potential Energy
As a further constraint, assume that 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 Furthermore, as M > 0, from (5), one has the following bound on p(t) in terms of the initial conditions: where λ max (M) is the maximal eigenvalue of M, and · indicates a Euclidean norm. If V (q) is radially unbounded [27], so that V (q) → +∞ as q → +∞, then a further bound follows for q . The non-negativity property of V is essential to invariant energy quadratisation methodsindeed, it can be generalized to the case of V bounded from below, so that V (q) ≥ c, for any constant c [20], but the non-negativity condition (6) above will be used here for simplicity.

Potential Energy Quadratisation
Under the non-negativity condition on V , from (6), one may define Hamilton's equations becomeq Furthermore, using the chain rule,ψ = (∇ q ψ) q .
The key feature of system (13) is that, if g is assumed known at at any given time instant, the resulting equations are linear in p, q and ψ. Though in the present case, g does indeed have a dependence on q, in the numerical setting, 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 (q), provided the non-negativity constraint (6) is satisfied.

Potential Energy Splitting
In some cases, the expression for the potential energy naturally takes the form where K is a symmetric positive semi-definite N ×N matrix, and where V ≥ 0. The first term could represent the stored energy of the system due to linear mechanisms, and V that due to additional nonlinear effects. Now, using one arrives at a system of equations generalizing (13): 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.

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 for a suitably chosen shift constant ≥ 0. Under the non-negativity constraint (6), V + is thus bounded away from zero, regularizing the calculation of g as defined in (12). The regularization approach applies equally to the case of a split potential, under the replacement V → V + .

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, . . .. A discrete-time vector u n represents an approximation to a continuous time vector u(t) at times t = t n .

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: 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 q n at n = 0 and n = 1 may be initialised, using the initial conditions q (0) and p (0) , as 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.

Explicit Energy-conserving Method (Marazzato et al.)
Consider the following scheme approximating (3), proposed by Marazzato et al. in [11]: This is a time interleaved scheme, with p n and q n+ 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 ∇ q V over the time interval t ∈ [nk, (n + 1)k], and over the known free-flight trajectoryq n (t) defined bỹ Other conservative schemes proposed also rely on such a continuous integration [16]. Except for particular functional forms of the potential energy V (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]: 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 ∇ 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 ≥ 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 ∇ q V is required per time step.

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 ψ as in (9), consider the following scheme, written in terms of the interleaved time series p n+ 1 2 and q n , and the scalar time series ψ n+ 1 2 : Here, g n is defined as The system (23) may be manipulated into an explicit update form in the following way. Beginning from (23a), one may write: where Here, I is the N × N identity matrix, and the vectors α n and β n are defined in terms of g n by α n = k 2 M −1 g n and β n = k 2 g n . Thus, given q n , q n−1 and ψ n− 1 2 , both A n and b n may be explicitly constructed. Notice that the matrix A n is positive definite by construction, due to the positive definiteness of M, and thus the update (26) always has a unique solution. Once the update in (26) has been performed, ψ n+ 1 2 may be computed from (23c) explicitly, using q n−1 , q n+1 and ψ n− 1 2 . Scheme (26) requires initial values for q 0 and q 1 , which may be set in the same way as for Störmer-Verlet, as in (19), and also ψ 1 2 which may be set to second order in k in terms of the initial displacement q (0) and momentum p (0) as Though the solution of (26) apparently requires an N × N linear system solution, in fact, the inverse is available in closed form, and allows for an explicit solution in O(N ) operations. A n is a rank-1 perturbation of the identity, and thus Sherman Morrison inversion [24] yields: Computational cost is thus on par with the other methods presented here. Notice in particular, though, that the solution of a linear system involving 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, 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 This may be identified as where the numerical energy H n+ 1 2 is 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 p n+ 1 2 may be bounded in terms of the energy H, for all time steps n, as which is identical to the bound (8) for the model system. The scheme (23) is thus unconditionally numerically stable.

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 ≥ 0, from which an auxiliary variable ψ 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: where now, An explicit update follows as in the case of the non-split form in (26), with A n as given in (27), but with b n now defined as An expression for a conserved numerical energy follows as before, now taking the form: 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: It then follows that Given that M > 0, M −1 may be factored into unique upper and lower triangular factors as The scheme is now conditionally stable, with the stability condition (40) corresponding to that for a linear system when V = 0.

Remarks
• Linear Conditions: One of the interesting features of the non-split form (23) is that, even if the potential V is quadratic in q, implying a linear system, the scheme is not linear. If, however, the split form (34) is used, then under linear conditions, V = 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 ψ n+ 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 where V nonlinear ≥ 0 and consists of higher order terms in q, then any splitting of the form with K = K 0 + K , K 0 ≥ 0, K ≥ 0 will yield a conservative form. The resulting numerical scheme, with V = 1 2 ψ 2 will inherit conservation of energy, which will be non-negative under a condition analogous to (40), but depending on K 0 .
• Generalized Update for q n : Consider again the non-split scheme (23), which may be rewritten as In this form, exact energy conservation follows from (43b) and (43c) only-it is independent of the values of g n , which are derived solely from q n . Further opportunities for generalization are thus available-(43a) could be replaced by any consistent update for 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 q n is a time series defined for t = t n , and similarly p n+ 1 2 and ψ n+ 1 2 are defined for t = t n+ 1 2 , then a scheme follows as: It is direct to show that, regardless of the choices of t n and t n+ 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 g in (24) and (35), through a replacement V → V + or V → V + , respectively.

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 transverselongitudinal 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, 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.

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, . . . , 2M . The system Hamiltonian is defined by 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 p and q of size N × 1, the scheme as presented in (23) follows immediately, with M = I N , the N × 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 where here, ⊗ 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

Numerical Results
We use the same settings as in [11], and choose ω = 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 10 −6 s. For initial conditions, we set p (0) = 0, and q (0) = [0, 0, 0, α, 0, 0] , for different values of α. 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). 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 α = 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.
As a basic test of convergence, consider the L 2 error defined, over a simulation duration n = 0, . . . , N f , by where q ref are computed using Störmer-Verlet with a time step of k = 2 −20 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 α = 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 α = 100, the range of time steps over which Störmer-Verlet is unstable is indicated in the figure.

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 Here, u = u(x, t) : D × R + 0 → R represents the transverse displacement, for a spatial domain D = [0, L], where L is the length of the unstretched string. Similarly, v = v(x, t) is the longitudinal displacement. In (50), ∂ t and ∂ x represent partial differentiation with respect to time t and spatial coordinate x. Furthermore, we define ζ = ∂ x u, η = ∂ x v. The function V = V(ζ, η) : R 2 → R + 0 is a non-negative potential density, which may be given in two equivalent forms as Here, the various constants that appear are: ρ, the volume density in kg·m −3 ; E, Young's modulus, in kg· s −2 · m −1 ; A, the area of the string cross section in m 2 ; and T 0 , the applied tension in kg·  m · s −2 . It is easy to show that (51a) is the same as (51b), and that V ≥ 0 ∀ (ζ, η). In (51b), one may easily split 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 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}.

Semi-discrete Form
The domain 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 ∂ x may be given as difference operators, expressed here in terms of their action on a grid function u l (t): From these, one may also define the second difference operator as (50) is then obtained as where V l V(ζ l , η l ). It is convenient to write system (54) compactly using the consolidated state vector q = [u , v ] . In this form, the difference operators are represented by matrices, such that D − becomes the M × (M − 1) matrix D − given in terms of its action on e.g. u as Then, define D + = −D − , and D 2 = D + D − . In vector form, (54) becomes where ζ = D − u, η = D − v. This system conserves the semi-discrete energy of the form (1), with Here I 2M −2 is the (2M − 2) × (2M − 2) identity matrix. A proof is obtained immediately by noting thatV Thus, left-multiplying (56) by hq yields (57). A split potential form is also available, via (51b). In this case,

Numerical Methods
System (56) may be integrated in time in a number of ways. First, introduce the discrete time vector q n , approximating q(t) at the time t n = kn. Furthermore, let V n = V (ζ n , η n ). An explicit time stepping scheme is obtained by application of the Störmer-Verlet algorithm: 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]: where the discrete gradients are defined as 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 This is the standard Courant-Friedrichs-Lewy stability condition for the one-dimensional linear wave equation [38].

Numerical Results
As a first example, let the string be initialised in its first linear mode of vibration in the transverse direction, that is and let the initial velocity be zero for both transverse and longitudinal motion. The amplitude parameter α is nondimensional. The string parameters are taken from [39], for the C3 piano string, and are: ρ = 7850 kg m −3 ; A = 8.87 · 10 −7 m 2 ; L = 1.259 m; E = 2.02 · 10 11 kg· s −2 · m −1 ; T 0 = 759 kg· m · s −2 . Figure 4 shows snapshots of the computed solution using scheme (34), and under various initial amplitudes α. In this case, a regularisation of the potential energy V is employed (see Section 2.5), with a shift of = 10 8 . Typical amplitude-dependent phenomena are visible: the frequency of vibration increases with the initial amplitude, and the initial shape deforms progressively. 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 α = 300, indicating that the output of the two schemes converges to a common solution.
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 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 4: Snapshots of the geometrically exact nonlinear string, computed using scheme (34), with split potential form as per (59). Here, k = 2.6 · 10 −7 , and the grid spacing is chosen as h = 1.05 E/ρ k. The string's initial normalised amplitude α is as indicated. In black is the reference solution, computed using scheme (60) with a time step k = 1.04 · 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 E/ρ k.

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: 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) ∈ D ⊂ R 2 , and time t ≥ 0. In this  (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 · 10 −7 , and the grid spacing is chosen as h = 1.05 E/ρ k. The string is initialised with using α = 300. simple example, the spatial domain D will be taken to be the square of side length L m, so that D = [0, L] 2 . ∂ t and ∆ represent partial differentiation with respect to time t and the 2D Laplacian operator, respectively. ∆∆ is thus the biharmonic operator. For simplicity, boundary conditions are assumed to be of simply supported type over the boundary ∂D of D, so that The various constants that appear in (65) are: ρ, the material density, in kg· m −3 ; ξ, the plate thickness, in m; E, Young's modulus, in kg· s −2 · m −1 ; and the flexural rigidity Q = Eξ 3 /12(1 − ν 2 ), where ν is Poisson's ratio for the plate material. L is a bilinear operator, defined in terms of its action on two functions f (x, y) and g(x, y) as where ∂ x and ∂ 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 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.

Semi-discrete Form
For the square region D, of side length L, one may start by defining grid locations x l = lh, y m = mh, where l, m = 1, . . . , 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 ∂ x and ∂ 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 Under simply supported conditions, when grid points outside the range l, m = 1, . . . , 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 purposesmore elaborate approximations (such as those of spectral type [48]) are available. From the basic operations defined in (69), centered approximations to the Laplacian ∆ and biharmonic operator ∆∆ follow as It is important to note that this particular construction of the biharmonic approximation D ∆∆ , through a product of Laplacian approximations D ∆ under fixed conditions ensures that the simply supported conditions (66) are satisfied. A semi-discrete approximation to (65) then follows as The operator (·, ·) approximates L(·, ·), as defined in (67). One useful centered approximation, operating on two grid functions f l,m and g l,m is [47]: It is useful to represent the semi-discrete ODE system (71) in vector form, using (M − 1) 2 × 1 vectors q and F. In this representation, the operators D x+ and D y+ become (M − 1) 2 × M (M − 1) matrices D x+ and D y+ , and D x− = −D x+ and D y− = −D y+ . The operators D ∆ and D ∆∆ become (M − 1) 2 × (M − 1) 2 matrices D ∆ and D ∆∆ respectively. One arrives at the form 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 p = Mq, the first order system (3) equivalent to (72) results from a Hamiltonian of the form of (1), with N = (M − 1) 2 , and where Here, A natural splitting of the potential energy V as in (14) follows, with

Numerical Methods
Beginning from the second order system (72), one may introduce the discrete time vectors q n and 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: The first of these updates is explicit, but relies on F n , which must be obtained from the second equation through the solution of a linear system involving the biharmonic operator D ∆∆ . 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: Due to the bilinearity of the operator , this scheme is linearly implicit-at each time step, q n+1 and 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 D ∆∆ , 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 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: The scheme, like Störmer-Verlet, relies on a linear system solution involving the biharmonic operator D ∆∆ , but is otherwise explicit.

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) = αξ sin(πx/L) sin(πy/L) ∂ t q| x,y,t=0 = 0 , where the dimensionless parameter α controls the maximum amplitude of the initial condition relative to the plate thickness ξ. Furthermore, the plate is assumed to be made of steel with E = 2×10 11 Pa, ρ = 7850 kg· m −3 , and ν = 0.3, and to be of thickness ξ = 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 α = 0.01, behaviour is essentially linear. At amplitudes near the plate thickness at α = 2, the period of oscillation decreases, and spontaneous mode generation is observed, and for large amplitudes, such as α = 10, turbulent behaviour is observed. Under such stringent high-amplitude conditions (α = 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 × 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 × 10 −5 s. 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 Figure 9: The output waveform q o (t) drawn from a plate vibrating at high amplitude, with α = 10. The reference solution, generated using Störmer-Verlet with k = 2.5 × 10 −6 s is shown in black; results are shown for scheme (34), in green, at different time steps k, as indicated.
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. 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.

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 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 α = 10. The time step is k = 10 −5 s. 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 (q) that are single-signed for all 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 (q) ≤ 0.
In this case, one may set, instead of (9), V = − 1 2 ψ 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 stepsuseful 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.