Role of shift constant in Energy Shifted SAV for Hamiltonian Systems

In recent years, considerable work has been devoted to the design of energy-stable numerical methods, a class of geometrical integration technique in which the discrete energy or pseudo-energy remains conserved. This property finds practical applications in systems for which the energy is bounded from below since the growth of the solutions is then itself bounded, yielding a form of stability. Such a property may be further exploited to transform the energy overall into a quadratic form, representing the core of recent numerical techniques such as the invariant energy quadratization (IEQ) and the scalar auxiliary variable (SAV) approaches. These methods have been applied to a large class of problems due to their remarkable efficiency. Yet, several aspects of such techniques have seen little investigation. In this work, the role of the “shift” constant in the expression of the potential energy in Hamiltonian systems is investigated. This is a positive gauge constant that may be used to augment the expression of the energy. Since Hamilton’s equations are given in terms of gradients of the Hamiltonian, such a constant has no influence on the resulting dynamics in the continuous system. However, empirical evidence has suggested that the convergence properties of the associated SAV approaches are affected by the magnitude of the shift factor. In this work, the behaviour of the relative global error of SAV schemes in the cases of the harmonic and Duffing oscillators is numerically investigated. The results reveal an optimal shift factor that increases the convergence rate. Using the optimal shift factor, the proposed method displays variable order accuracy ranging from second to twelfth order.


Introduction
Hamiltonian systems are widespread in physics, describing a large class of dynamical systems.Under time-invariant or autonomous conditions, the Hamiltonian corresponds to the physical energy of a given system, and is conserved.In the simulation setting, mimicking this property through a conserved numerical energy (or pseudo-energy) conservation is, thus, a desirable feature.In many numerical designs, the enforcing of such a conservation property is achieved by means of fully implicit schemes requiring the iterative solution of nonlinear algebraic equations at each time step, or via the exact evaluation of continuous-time integrals [1,2].In the former case, the iterative nature of the time-stepping routine often represents a computational bottleneck, while in the latter, exact numerical energy conservation is achieved in the few cases where the integrals can be computed exactly.Recently, fully explicit, energy-stable numerical designs for a class of Hamiltonian systems have been proposed [3], overcoming all such limitations.These designs are part of a larger class of techniques rooted in energy quadratisation, and appearing under various guises such as the Invariant Energy Quadratisation (IEQ) [4,5], and the Scalar Auxiliary Variable (SAV) methods [6].These methods have seen a growing body of applications, particularly in gradient flow systems (see e.g.[7,8]), and including diagonally-implicit Runge-Kutta methods [9].
Hamilton's dynamical equations are derived from the gradients of the Hamiltonian, and thus the equations are gauge-invariant: that is, Hamiltonians differing by a constant yield the same equations.Numerically, however, gauge constants have been shown to influence the behaviour of the numerical error, as in the Navier-Stokes simulation routines presented in [10].The influence of the gauge constant in SAV-like methods for Hamiltonian systems will be the subject of this article, building on previous work by the authors [3].A rule for selecting the shift constant is tested here on two scalar problems: the simple harmonic oscillator and the Duffing oscillator (a particularly simple nonlinear oscillator).The experimental analysis reveals that an optimal shift factor exists for a given potential, which can improve the local convergence rate up to several orders, and which may be obtained as a fraction of the system's energy.
The manuscript is organized as follows: Section 2 describes the proposed method, applied to the test cases described in Section 3. Section 4 includes the numerical tests and the estimates of the optimal gauge constant's values.

Energy-shifted Hamiltonian systems
Hamilton's equations are a set of ordinary differential equations (ODEs) giving the rate of change of the displacements and momenta of a system of particles.In this work, a single particle with position q(t) and momentum p(t) is considered, for which: Here, H = H(q, p) : R 2 → R is the Hamiltonian, restricted here to the form: where V (q) is the potential energy of the particle.Here, and in the following, the Hamiltonian is scaled by the particle's mass for the sake of conciseness.It is furthermore assumed that: implying H ≥ 0 ∀ (p, q).The non-negativity of the potential V is verified in many Hamiltonian systems, though not all -amongst others, the gravitational potential is an exception.When (3) holds, the potential energy may be written as [3]: where ψ ∈ R + 0 is referred to as the "auxiliary variable" [3,4].Consequently, Hamilton's equations take the form: q = p ṗ = −gψ . ( with g := ∂ψ ∂q .The equations in (5) can be combined through differentiation, and the rate of change of the auxiliary variable may itself be obtained by a simple application of the chain rule.Together, these yield the following system, forming the basis for SAV-like approaches: Initialisation is provided by: q(0) = q 0 q(0) = p 0 ψ(0) = 2V (q 0 ), (7) where q 0 , p 0 are constants.Note that the same equations can be obtained by replacing V with V ϵ in (2), a shifted potential defined as: where the non-negative shift ϵ has the interpretation of a gauge constant.

Time stepping routine
System ( 6) is now integrated in time using the method proposed in [3].This is a finite difference scheme discretizing the time domain into uniformly-spaced grid intervals at the times t n = nk, where n ∈ N is the time index, and k is time step.This grid is used to compute the discrete-time position q n ≈ q(t n ).Moreover, an interleaved grid centered at the midpoints t n + k 2 is employed to approximate the auxiliary variable ).An approximation to ( 6) is obtained as: where: All finite difference approximations in the scheme are centered, making the scheme reversible and (at least) second-order accurate.Furthermore, it possesses a conserved, non-negative energy of the form: from which bounds on |q n − q n−1 | and |ψ n− 1 2 | in terms of the initial energy h 1 2 are easily derived, leading to unconditional stability.
The numerical initialisation of scheme ( 8) is given in terms of the exact solution q(t), as: Whilst q is here a scalar, scheme (8) may be generalised to the vector case; an explicit update exists, making the scheme particularly attractive from the standpoint of real-time simulation [3,[11][12][13].

Test Cases
Two test problems are considered, with closed-form solution.For both problems, in order to simplify the analysis, null velocity initial conditions are considered, such that p 0 = 0 in (7).

Simple Harmonic Oscillator
The first test case is the classical simple harmonic oscillator, for which: where ν > 0 is the oscillator's angular frequency in radians per second.An exact solution to (6) exists as: q(t) = q 0 cos(νt) .
Substituting the definition of V (q) as per ( 12) in ( 9), one has: Duffing Oscillator The second test case is the Duffing oscillator, a nonlinear oscillator with a mixed linear/cubic restoring force.It models various oscillatory systems, such as the simple pendulum under a moderately large vibration amplitude [14] and, in the distributed case, it governs the dynamics of plates, shells and strings [12,15].In the scalar case considered here, one has: The sign of γ may be positive or negative, leading to a "hardening" or a "softening" nonlinearity, respectively.Because ( 3) is required to hold, only the case γ ≥ 0 is considered here.An analytic solution exists, generalising (13), as: Here, cn(z; m) is the Jacobi elliptic function of argument z and parameter m.Using ( 15) in ( 9), the nonlinear gradient is given in this case by:

Numerical Results
To investigate the impact of the shift constant on the convergence properties of scheme (8), the influence of ϵ on the relative error is first assessed.The relative error is here defined as: that is, the relative difference between the output of scheme (8) at the final time step N and the corresponding exact solution q(t N ).Note that q n , from (8), depends on ϵ via g n ϵ , and, hence, so does δ N ϵ .Then, define: provided this minimum exists.Here, ϵ ⋆ will be sought as a fraction of h 0 the total numerical energy in the absence of shift, such that: with The total energy h 0 provides in this case a useful scaling factor.Hence, the global minimum will be defined in terms of the scaled shift, as:

Simple Harmonic Oscillator
The behaviour of the relative error δ N ϵ is first checked as a function of η, for a set of frequencies Ω = ν | ν ∈ (1, 2, 5, 10, 20, 50) • 10 2 , and under various choices of the time step k.It is convenient to perform a check on a coarse grid along the η axis, and then to restrict the search to a useful range.The coarse search reveals the presence of minima in the range 0.01 ≤ η ⋆ ≤ 0.03, as seen in Figure 1.An analysis of this figure reveals that a single (and, hence, global) minimum Figure 1.Coarse behaviour of the relative error δ N ϵ as a function of the scaled shift η, under various choices of the time step k, given on top of each panel, and linear frequency ν, given by different colours.Colour scheme: ν = 100 (blue), 200 (red), 500 (yellow), 1000 (purple), 2000 (green), 5000 (cyan).Here, N = floor(0.3/k),and the η axis is sampled using 500 linearly spaced grid points between 0.01 and 0.03.When 10 −6 ≲ k ≲ 10 −5 , a single, global minimum is detected ∀ ν.For larger values of k, a single minimum is detected only for frequencies such that νk ≲ 0.1 whereas for smaller values of k, the effects of round-off are visible as a noise-like distortion.For all tests, the system is initialised with q 0 = 1, p 0 = 0. exists for δ N ϵ when the time step k becomes sufficiently small compared to the radian frequency ν, such that νk ≲ 0.1.When 10 −6 ≲ k ≲ 10 −5 , a clear global minimum appears for all frequencies.The definition of the scaled shift η as per (20) results appropriate from the inspection of this figure, where the minimum points η ⋆ appear to be more or less aligned for all frequencies, though a small dependency of η ⋆ on ν is nonetheless observed.Interestingly, for smaller values of the time step k, the behaviour of δ N ϵ is dominated by round-off errors appearing as a noise-like distortion.
that, while Figure 1 displays the case q 0 = 1, the minimum points η ⋆ are independent of the magnitude of the initial condition q 0 for fixed ν.Furthermore, such values for η ⋆ are independent of k for fixed ν, provided that k is small enough to yield sufficient resolution, but large enough to keep away from round-off errors.In practice, the optimal scaled shift value appears to be a function of the radian frequency only, so that η ⋆ = η ⋆ (ν), and is thus independent of k and q 0 .Let W = {η | η = η ⋆ (ν), ν ∈ Ω} denote the set of the computed minimum points from Table 1.The convergence rate of scheme ( 8) as a function of the time step k is now checked, as per Figure 2.Each panel considers one frequency from the set Ω, and the nine coloured lines correspond to the eight values of the shift constant obtained from W∪{0, 10 12 } (that is, from the set E = {ϵ | ϵ ∈ (W ∪ {0, 10 12 }) h 0 }), plus the second-order leapfrog scheme, used as a reference.For all panels, as expected, the value of δ N ϵ is considerably smaller in correspondence of ϵ ⋆ .When k ≈ 10 −5 (a typical value of the time step in acoustics simulation), the relative error is several orders of magnitude smaller for the energy-shifted SAV scheme compared to the leapfrog scheme.These results suggest computing ϵ ⋆ with a precision at least 4 significant digits, since a smaller precision yields rather large amplifications of the relative error δ N ϵ .This conclusion is supported by the behaviour of the curves in Figure 1, all presenting a very steep descent in the neighbourhood of η ⋆ .
The dependence of the convergence rate on ϵ is another evident feature of the curves in Figure 2, and one of the most interesting findings of this work.Whilst scheme (8) remains, generally, second-order convergent, the rate of convergence seems to be largely affected by the magnitude of the shift factor.Locally, slopes range from 3 to 12 when ϵ = ϵ ⋆ .When ϵ is chosen in some neighbourhood around ϵ ⋆ , some improvements are still observed for larger time steps.This suggests that increasing the precision on the estimate of ϵ ⋆ may improve the convergence rates further, though this behaviour can only be understood by a formal analysis of the structure of the error, left as future work.
The limiting cases ϵ = {0, ∞} are also worth commenting.It is remarked that the somewhat "natural" choice ϵ = 0 yields a very slow convergence rate, if any: for larger values of the oscillator frequency ν, the convergence of the scheme seems to be completely impaired.On the other hand, a very large value of ϵ produces a scheme virtually indistinguishable from the leapfrog algorithm, as one can show immediately from (8) (in this case, g n ϵ ≈ ( √ ϵ) −1 ∂V /∂q, ψ n+ 1 2 ≈ √ ϵ, from which the leapfrog is recovered).Before proceeding, it is worth remarking that the case ϵ = 0 may be treated differently.First, note that, according to definition (14), g n ϵ = ν sign(q n ) when ϵ = 0. Whilst ψ, from (4) was restricted to be a positive constant, it may well change sign according to ψ = ± √ 2V + ϵ.Selecting the sign accordingly, one may then get g n ϵ = ν, yielding a different three-step scheme.The assessment of this scheme is left as future work.

Duffing oscillator
The analysis of the Duffing oscillator may be done analogously.It is convenient in this case to keep the linear frequency ν fixed, and to study the behaviour of the scheme for varying (γ, q 0 ): both these parameters influence the amplitude-dependent nonlinearity of the system.In the following, γ ∈ {0.1, 1, 10, 100} and q 0 ∈ {0.1, 1, 10}.
First, a coarse check of the η axis is performed, as per Figure 3.In the figure, plots of δ N ϵ under two different choices of the initial condition q 0 and time step k are given.Whilst a single, global minimum is recovered for fixed ν, γ, k, the value of η ⋆ now clearly depends on the intial condition q 0 .Thus, η ⋆ = η ⋆ (q 0 , ν, γ).Table 2 reports values of the minimum points η ⋆ under various choices for γ, q 0 .Note that, as expected, for small values of the nonlinear parameter γ and q 0 , the same value as from Table 1 is recovered (η ⋆ = 0.02171 for ν = 500).
Let now V = {η | η = η ⋆ (ν, γ, q 0 ), ν = 500, γ ∈ {0.1, 1, 10, 100}, q 0 ∈ {0.1, 1, 10}} be the set of the optimal scaled shift values from Table 2.In Figure 4, the behaviour of the relative error δ N ϵ is assessed against the time step k, using the values of the shift constant from V ∪ {0}, plus the leapfrog scheme.Much like the simple harmonic oscillator case, it can be appreciated that the convergence rate improves drastically whenever ϵ = ϵ ⋆ .Values in the neighbourhood of ϵ ⋆ also lead to a considerably faster convergence for larger k compared to the leapfrog.Again, the no-shift case (ϵ = 0) leads to poor or no convergence.Nonlinear parameter γ and intial condition q 0 as indicated.For the numerical search, k = 3 • 10 −6 was used, and the η axis was scanned between 0.015 and 0.03 using 5000 grid points.
Figure 3. Coarse behaviour of the relative error δ N ϵ for the Duffing oscillator as a function of the scaled shift η, under various choices of the time step k, given on top of each panel, and γ, given by different colours.Colour scheme: γ = 0.1 (blue), 1.0 (red), 10 (yellow), 100 (purple).Here, N = floor(0.3/k),and the η axis is sampled using 500 linearly spaced grid points between 0.015 and 0.03.For the tests, the system is initialised with q 0 as indicated, showing the dependence of η ⋆ on the initial condition.

Conclusions
In this work, an assessment of the gauge constant in energy-shifted SAV methods for Hamiltonian system was given.The scalar cases of the simple harmonic oscillator and the Duffing oscillator were treated in detail.In both cases, the value of the shift affects the convergence rate, and numerically computed optimal values were given as a fraction of the system's energy.For the harmonic oscillator, a linear problem, such optimal value depends exclusively on the linear oscillator frequency, whilst for the Duffing oscillator, a nonlinear problem, such optimal value depends on the amplitude of the initial condition as well as the linear frequency and nonlinear parameter.When the shift is chosen in the neighbourhood of the numerically computed optimal values, local changes in the convergence rate up to twelfth order are observed, as well as a large reduction of the magnitude of the relative error compared to standard algorithms such as the leapfrog.On the other hand, when the schemes are run with no shift, the error does not decrease with decreasing time step, and hence convergence is halted.These findings suggest studying the behaviour of relative error as a function of the shift analytically, which is left as future work.

Table 1 .
Numerically-computed global minimum points from Figure