Geometric numerical integration of Li\'enard systems via a contact Hamiltonian approach

Starting from a contact Hamiltonian description of Li\'enard systems, we introduce a new family of explicit geometric integrators for these nonlinear dynamical systems. Focusing on the paradigmatic example of the van der Pol oscillator, we demonstrate that these integrators are particularly stable and preserve the qualitative features of the dynamics, even for relatively large values of the time step and in the stiff regime.


Introduction
Liénard systems are a class of 2-dimensional nonlinear dynamical systems that exhibit a stable limit cycle. Among them the most famous is the van der Pol oscillator [19,26].
Due to the existence of a stable limit cycle, such systems are of the utmost importance in modelling natural phenomena such as e.g. electrical circuits and neuronal dynamics, and therefore an accurate investigation of their dynamics is required. However, because of the nonlinear nature of such systems, analytical results are scarce and one has to recur to perturbative techniques and numerical integration.
An immediate and paramount problem for both the development of perturbative techniques and of stable numerical schemes is the lack of a geometric structure. Indeed, apart from very specific cases in which some integrability conditions are satisfied, and where one can use the Jacobi Last Multiplier to find a Lagrangian or Hamiltonian structure [10,22], in the general case such pursuit is hopeless. For instance, many Liénard systems present an attractor, a stable limit cycle, and thus they cannot be Hamiltonian in the symplectic sense. There have been several attempts in the literature in order to circumvent this problem. In [12] the authors suggested to enlarge the phase-space to a 4-dimensional manifold and define a particularly simple Hamiltonian system in this enlarged space so that the 2-dimensional projection onto the original space recovers the original dynamics, and then they showed that this approach allows for the use of perturbative methods. In [27], the classical Bateman trick for the harmonic oscillator has been extended to the van der Pol oscillator and then further generalised to all Liénard systems with a quadratic potential. Both these approaches involve a 4-dimensional phase-space and in both the authors have focused on the perturbation theory and have not explored the consequences of the Hamiltonisation for the numerical integration. From yet another perspective, in [11] the authors have presented various splitting schemes for "conditionally linear systems" (these include Liénard systems) which, although not being geometric, are based on the standard splitting schemes for symplectic Hamiltonian systems, and showed good qualitative and quantitative results.
In this work we contribute to the advancement of geometric integration for Liénard systems by using Hamiltonian flows on contact manifolds. Contact geometry was introduced in Sophus Lie's study of differential equations, and has been the subject of an intense research, especially related to low-dimensional topology [17]. In recent years, contact Hamiltonian systems have found many applications, first in the context of thermodynamics [7,21,29] and, more recently, in the context of the Hamiltonisation of several dissipative dynamical systems [6,8,9,13,15,16,30]. The large number of applications of contact systems that have appeared recently motivated research on geometric numerical integration [9,28,30]. Fortunately, contact flows possess geometric integrators (both variational and Hamiltonian) that precisely parallel their symplectic counterparts, and therefore they show remarkable numerical and analytical properties such as e.g. increased stability, near-preservation of invariant quantities, and modified Hamiltonians.
In this work, leveraging some of the ideas in [12], we start a treatment of Liénard systems from the point of view of contact Hamiltonian systems: we show that they can be given a particularly simple Hamiltonian formulation on a 3-dimensional contact manifold, and then we use this Hamiltonisation to construct splitting integrators for such systems and analyse their properties from an analytical point of view, exploiting the modified equations. Along the work we use the van der Pol oscillator as a paradigmatic example.
Our results show that the resulting geometric integrators are very stable, even when the system is stiff, and they preserve the qualitative features of the limit cycle even for large values of the time step, which permits to spare computational resources and is of primal importance in applications to e.g. neuronal dynamics [11]. Moreover, from the use of the modified equations, we can prove analytical results on the preservation and the period of the limit cycle that show a very good agreement with the numerical simulations.
The paper is organised as follows: in Section 2 we provide a Hamiltonian formulation of Liénard systems based on contact Hamiltonian dynamics, and then in Section 3 we introduce a new class of explicit geometric integrators for these systems that are naturally derived by splitting the Hamiltonian. Then in Sections 4 and 5 we thoroughly analyse the properties of these integrators both analytically and numerically by investigating the benchmark example of the van der Pol oscillator. We conclude in Section 6 with a discussion and a perspective on future work.
All the simulations are reproducible with the code provided in [33].

A brief review of Liénard Systems
Liénard systems are a family of planar coupled differential equations of the form [19] where F (x) is the antiderivative of an even function f (x) and g(x) is an odd function. Alternatively, (1) is equivalent to the second order scalar equation A third equivalent version of (1) is given by Example 1 (The van der Pol oscillator). Perhaps the most famous example of the family of Liénard systems is the van der Pol oscillator, which can be written using dimensionless variables as followsẍ and can be equivalently rewritten in the form (3) as from which we recognise that in this case f (x) = − (1 − x 2 ) and g(x) = x.
A crucial property of Liénard systems is encoded in the following theorem, guaranteeing the existence and uniqueness of a stable limit cycle for a large class of systems [24].
• F (x) has exactly one positive zero at x = a, is monotone increasing for x > a and lim x→+∞ F (x) = +∞; the dynamical system (1) presents a unique, stable limit cycle.
In particular, the theorem above implies that the van der Pol equation (4) with > 0 has a unique, stable limit cycle.
For additional information on the classical approach to the analysis of Liénard systems we refer to [24].

A brief review of contact Hamiltonian systems
Similarly to the fact that a symplectic manifold is a 2n-dimensional differentiable manifold endowed with a 2-form ω that is closed (dω = 0) and non-degenerate (ω n = 0), an exact contact manifold M is a (2n + 1)-dimensional manifold endowed with a 1-form η, called the contact form, that is non-degenerate, which means A contact version of Darboux's theorem [3] guarantees the local existence of coordinates (q i , p i , s) -called Darboux coordinates -which permit to express the contact form as η = ds − p i dq i , where Einstein's summation convention over repeated indices is being used here and in the following.
The contact form allows us to define in a natural way the concept of a Hamiltonian vector field on M . Let H be a real function on M , then the contact Hamiltonian vector field X H associated with H is defined by where ι X is the interior product and R is the Reeb vector field corresponding to η [17]. In Darboux coordinates X H takes the form Finally, contact manifolds carry a natural bracket structure, called the Jacobi bracket, which yields a Lie algebra on smooth functions on M and is defined as Again, in Darboux coordinates the Jacobi bracket reads We refer the reader to [3,7,8,16,18] for further details. For our scope, it will be important in the following to have an explicit expression for the Jacobi bracket of monomial functions, that is, where µ,μ ∈ R and i, j, r,ī,j,r ∈ N.

A contact Hamiltonian formulation of Liénard Systems
It is well known that any dynamical system on an n-dimensional manifold Q of the forṁ can be extended to a Hamiltonian system defined on the 2n-dimensional phase-space T * Q. This can be achieved with the introduction of the conjugate momentã p i in order to define the Hamiltonian A direct computation shows that when we consider only the dynamics on the original x-variables, then we recover the original n-dimensional system. For instance, in the case of Liénard systems (3), the Hamiltonian reads In [12], such approach has been used to derive a Hamiltonisation of Liénard systems in such extended phase-space that was then shown to be useful to perform perturbation theory. Moreover, in [25] a similar extension, but with a suitably defined new Hamiltonian that non-trivially couples the variables, has been used in order to develop geometric integrators in the extended phase-space and then used e.g. in the case of the van der Pol oscillator.
In principle one could use the Hamiltonian (13) and perform a splitting in order to obtain new geometric integrators that are symplectic in the extended phase-space. However, we see from the form of (13) that it is linear in the momenta, meaning that it is naturally associated with a contact Hamiltonian on the (2n − 1)-dimensional projectivised cotangent bundle P T * Q, endowed with the contact structure inherited from the canonical symplectic structure of T * Q [4,5]. The procedure to perform such reduction is quite simple in this case and it is reviewed e.g. in the recent work [29]. In order to avoid clutter of notation, from now on we focus on the case Q = R 2 , which is the relevant case for our study: we start with (13) and consider a connected component of the open set in whichp 2 = 0. On such set we can define the coordinates (q = x, s = y, p = −p 1 p 2 ), which serve as Darboux coordinates on P T * R 2 . Finally, we define the contact Hamiltonian A direct calculation then shows that the restriction of the resulting contact Hamiltonian system to the (q, s) plane recovers the original system. By means of the above prescription, we arrive at the following result for Liénard systems.
The associated contact Hamiltonian system iṡ From the first two equations we recover the original Liénard system in the (q, s)-space, while the third equation is decoupled.
Example 2 (The van der Pol oscillator revisited). As we have already seen in Section 2.1 the van der Pol equation is a particular case of a Liénard system, which is obtained by choosing f (x) and g(x) as Consequently the contact Hamiltonian in this case reads and the corresponding contact Hamiltonian systems is As expected, from the first two equations we recover the original van der Pol equation (4).
Remark 2.1. For s = 0 and setting the appropriate initial condition p 0 = −f (q 0 ) − g(q 0 )/s 0 , p(t) derived from (18) turns out to be the slope of the tangent ds dq to the orbit of the system at each point (q(t), s(t)) of its evolution. This stems from the fact that (16)- (18) are the characteristic equations of the Hamilton-Jacobi equation for (15). Details of this derivation are in preparation by [20].
Remark 2.2. The reduction procedure that led us to (14) is not unique. Indeed, we could have selected the connected component in whichp 1 = 0 and set (q = y, s = x, p = −p 2 p 1 ). The corresponding contact Hamiltonian for Liénard systems is: Beware that in this case X 1 (q, s) = q and X 2 (q, s) = −f (s)q − g(s), that is, the roles of q and s are switched, and the resulting system is which is equivalent to (16)- (18) for the (q, s) part, but not so much for p.
The choice of reduction, in the case at hand, was dictated by numerical convenience: the Hamiltonian H from (14) resulted in a simpler form of the algorithm providing better results.

Contact splitting integrators
Contact splitting integrators are a class of geometric integrators recently introduced in the context of celestial mechanics [9]. They are the contact analogues of the well-known symplectic splitting integrators.
Let H be a contact Hamiltonian which is separable into a sum of functions Then, the Hamiltonian vector field associated with H is separable as well If moreover, each of the X h j is exactly integrable, meaning that there exists a closed-form solution for its flow, then we can approximate the dynamics of X H to second order in τ with contact maps according to the following proposition. Proposition 3.1 (Contact splitting integrators). In the hypotheses above, let e tX h j denote the map given by time-t exact flow of each vector field X h j , for j = 1, . . . , N . Then is a second order contact numerical integrator, meaning that each map is a contactomorphism.
From knowledge of the second order contact integrator (27) and using Yoshida's standard formulation for the composition [32], we can construct two types of contact integrators of any even order; the difference between the two methods is that one involves exact coefficients for the calculation of the new time step, while the other uses approximated coefficients and involves a smaller number of map computations per iteration. The two methods are summarised in the following propositions.
with z 0 and z 1 given by is an integrator of order 2n + 2.
is an integrator of order 2n.
In Table 1 we list the values of the approximated coefficients {w j } m j=0 for three different 6th order integrators, labelled as A, B and C. Note that w 0 : Remark 3.1. The splitting integrator with approximate coefficients labeled as A is the better performer among the approximate splitting integrators of 6th order presented here. This can be related to the fact that its largest coefficient is the smallest among the approximate integrators.

Modified Hamiltonian and error analysis
One of the main advantages of using contact splitting integrators is the possibility to have a direct error control by using the modified equations obtained from the modified Hamiltonian that results from the Baker-Campbell-Hausdorff (BCH) formula (see [9] for further details on the derivation of the modified Hamiltonian in the contact case). Indeed, for an integrator of order 2n multiple applications of the BCH formula give [9] S 2n (τ ) where all the corrections X 2i+1 are Hamiltonian vector fields. Therefore S 2n (τ ) is the time-τ flow of a Hamiltonian vector field, and its associated Hamiltonian, called the modified Hamiltonian, can be written formally as the power series where the subscript 2n in H mod,2n denotes the fact that it is associated with an integrator of order 2n, and ∆H 2i are the Hamiltonian functions associated with the Hamiltonian vector fields X 2i+1 , that is, Plugging (32) into the contact Hamiltonian equations that stem from (8), we obtain the modified equations, which are the equations whose time-τ flow gives exactly the integrator S 2n (τ ). Therefore studying the modified equations and their relation with the original equations gives us important information on the modifications introduced by the integrator on the original system.

Geometric numerical integration of Liénard systems
The application of the contact splitting integrators introduced in Section 3.1 to Liénard systems starts with the splitting of the contact Hamiltonian (15) as and the consequent identification of the corresponding vector fields The structure of this splitting ensures the exact integrability condition for any choice of the functions f (q) and g(q). Indeed, the time-τ flow maps are explicitly given by and the corresponding time-τ flow maps are In the next section we present the numerical and analytical results of the application of various splitting integrators based on the maps (40) to the van der Pol oscillator. To fix the notation, when referring to a particular splitting, we will write e.g. S 2 (τ )(CBABC) to indicate that we are using the 2nd order integrator obtained using the splitting (27) of the maps (40) composed in the order indicated in parentheses.

Geometric numerical integration of the van der Pol
oscillator: numerical vs analytical results

Numerical results
We split the analysis into three different cases, labelled by the value of the nonlinear coupling parameter : for = 0 we recover the harmonic oscillator on the plane (q, s); for 1 and ∼ 1 we are in the non-stiff regime; for 1 we are in the stiff regime. It is well-known that to approximate the limit cycle with Euler-type methods, one cannot choose the time step τ independently of , even in the non-stiff case 1 [11]: for example the Euler method requires τ and the exponential midpoint method requires τ 3 . In the rest of this section we will focus on the performance of our algorithm in the preservation of the limit cycle. As we will see, our methods accurately preserve the limit cycle of the van der Pol oscillator when τ 1: this allows for much larger step sizes than Euler-type methods when integrating Liénard systems. Figure 1 shows the solutions in the (q, s)-plane for different time steps τ and with the same initial condition (q 0 , p 0 , s 0 ) = (2, 0, 0). We can observe that the integrator is stable at least until the surprisingly large value τ ∼ π/2 > 1. By increasing the time step the typical circular orbit of the harmonic oscillator becomes more elliptic, and the period changes. In Figure 2 we plot the relation between the time step and the period of the orbits obtained from numerical simulations. Even though the frequency changes, we can see that the variation remains well under control for all values of τ ∈ (0, 1]. In Figure 3 we show the persistence of the limit cycle for different values of ∈ {0.1, 0.5, 2, 4} and τ ∈ [π/256, π/2]. Clearly the limit cycle is preserved also for very large values of and τ in this range. Moreover, the very long integration time, with t ∈ [0, 10000], is an evidence of the stability of the integrator. Finally, the dependence of the period and the frequency of the limit cycle with respect to the time step shown in Fig. 4 is very similar to that of the harmonic oscillator.

1 (stiff regime)
To better understand what happens in the stiff case 1, it is convenient to perform, after the integration, the so-called Liénard transformation [11,19] This change of variables transforms the dynamics into and enables a nice geometric description of the limit cycle. Indeed, the q nullcline, which is the locus of points such thatq = 0, is given by the cubic s = q − q 3 3 . Since q evolves much faster than s, the solutions are quickly attracted by the cubic nullcline. Once there, they move slowly along the curve until they reach an extremum, at which point they quickly jump horizontally to the other branch of the nullcline. This periodic motion that jumps back and forth on the nullcline is the attractive limit cycle of the stiff van der Pol oscillator. Figure 5 shows the cubic nullcline and the numerically simulated attractor for ∈

Analytical results
In this section we provide an analytical study of the contact splitting integrators for the van der Pol oscillator based on the modified equations. We start by providing two general properties of the modified equations that are of special importance.
As we have seen in Example 2, in the contact formulation of the van der Pol oscillator the equations for q and s are independent of p, as it should be. Clearly, given that the maps for q and s in (40) are all independent of p, any splitting integrator will satisfy this property too. However, it is instructive to recover this result by using the modified Hamiltonian, since in the proof we will find out an important property of H mod , i.e. that it is linear in p, as it is the original Hamiltonian (20). This is the content of the next result. Proof. We prove first the second part: the claim is that if H mod is linear in p, then the corresponding modified equations for q and s do not depend on p. By a direct look at the general contact Hamiltonian equations (8), this is clearly true. Now let us prove that H mod is indeed linear in p: considering the splitting in (39), we have that A = − 1 − q 2 s, B = q, and C = ps, are all polynomials in q, p, s and that only C depends (linearly) on p. Therefore, we see from (11) that by commuting A, B and C we can only obtain terms that are at most linear p. Then again, by commuting two terms that are at most linear in p, we see from (11) that we always obtain terms that are at most linear in p. We conclude that the modified Hamiltonian is at most linear in p.
We conclude that H mod is indeed linear, because otherwise in the modified equations we would haveq = 0, which is clearly not the case.
Furthermore, we observe that when the time step τ = 0 any truncation of the modified equations is likely to possess new spurious equilibria. This is so since at any order the corresponding vector fields are polynomials in q, p, s of increasing order. Therefore it is important to actually prove that (q, s) = (0, 0) is the only fixed point (considering only the dynamics projected to the (q, s) plane) for the integrator and that it is unstable, as we show in the next result.
Proof. The proof is based on writing explicitly the action of the integrator on an initial condition, that is, we apply e τ /2X C e τ /2X B e τ X A e τ /2X B e τ /2X C to (q i , s i ), to obtain Now when we impose the condition for a fixed point using the second equation in (44) into the first equation in (43) we obtain which is true if and only if s i = 0.
Next, we substitute s i = 0 = s i+1 into the second equation in (43) and we obtain which is true if and only if q i = 0.
To prove that (0, 0) is unstable, we compute the Jacobian of the map (43) at (0, 0), and in particular we obtain that its determinant is e τ > 1, indicating that at least one eigenvalue has absolute value > 1, which proves the instability.
To conclude the proof, let > 0 and τ > 0. A direct computation shows that the eigenvalues of the Jacobian of the map (43) at (0, 0) are Depending on the sign of β we have two cases: the eigenvalues are both real or they are complex conjugates.
In what follows we split the analysis into three different cases depending on the value of , as we did in the Section 4.1.

= 0 (harmonic oscillator)
In this case we have a harmonic oscillator, for which each nontrivial trajectory has period T = 2π. Moreover, the maps (40) in this particular case are simplified (for instance, the map e τ X A becomes the identity) and the modified Hamiltonian takes the remarkably simple expression H mod,2 = psF (τ ) + q G(τ ) where The corresponding modified system is thus which is again exactly solvable (recall that τ is fixed), and the solution in q and s is a harmonic oscillator with frequency In Fig. 6 we compare (56) with the numerical results for the period and the frequency obtained in Section 4.1.1. We observe that there is a very good agreement between the analytical expression up to the 8th order in τ and the numerical results.

1 (non-stiff regime)
This regime can be studied using perturbation theory and therefore there are many results (see e.g. [1,2]). We study the persistence of the limit cycle for the contact splitting integrators in a way similar to [11], that means, we use the modified equations in order to provide some estimations on the amplitude and period of the limit cycle.

Proposition 4.3.
For any contact splitting integrator of order 2n based on the maps (40), the projection of the numerical solutions of the van der Pol system (21) onto the (q, s)plane have a limit cycle at the approximate radius r = 2 + O(τ 2n ). Moreover, the approximate radius of the S 2 (τ )(CBABC) integrator, up to order 4 in τ , is Proof. Let us consider a contact splitting integrator S 2n (τ ) of order 2n; using the BCH formula (see Section 3.2) we can argue that the modified Hamiltonian whose time-τ flow is given by S 2n (τ ) is of the form Thus the modified equations read We know form Proposition 4.1 that the equations forq andṡ are independent of p, and from Proposition 4.2 that the point (0, 0) in the (q, s)-plane is an unstable equilibrium of the system. If we rewrite the system in polar coordinates on the plane (q, s) with the change of variables q = r cos θ and s = r sin θ, then the equation forṙ readṡ r = r sin 2 (θ) 1 − r 2 cos 2 (θ) + τ 2n R 2n (r, θ) Since the modified Hamiltonian is by construction a polynomial in the variables q and s, the dependence on θ of R 2n is only through sums and products of trigonometric functions. In particular, this implies that the averaged dynamics ofṙ obtained by the integration along a period has the form 1 2π One now observes that, modulo high order terms in τ , the stationary points of the averaged dynamics are r = 0 and r = 2, which implies that the latter is the radius of the limit cycle, proving the first part of the theorem.
For any fixed order, it is possible to give a more refined estimate of the limit cycle radius by looking at the exact correction from the modified Hamiltonian.
In the non-stiff regime, we can also perform a perturbative analysis by applying the Poincaré-Lindstedt method to study the frequency (and hence the period) of the system (see e.g. [2]). The first step consists in the time reparametrisation t = ω t, which leads to the differential equation where the derivatives are now expressed in terms of t , instead of t, and, as usual, we omit the decoupled equation forṗ. Noticing that the modified Hamiltonian vector field depends on the two parameters and τ , we suppose, in analogy to the traditional approach [2], that all the terms appearing in the equations can be expanded in Taylor series with respect to such parameters as follows In particular, notice that we assume all the expressions to be of even order in τ , given that all the terms appearing in the modified equations are of even order. For convenience, and without loss of generality, we follow [2] and assume that This is equivalent to a convenient time shift that simplifies the initial conditions. The differential equation corresponding to the order 0 , τ 0 then reads ω 0,0 q 0,0 (t ) = s 0,0 (t ) ω 0,0 s 0,0 (t ) = −q 0,0 (t ) (70) whose solution is Since we want q 0,0 (t) and s 0,0 (t) to have period 2π, this fixes ω 0,0 = 1, while condition (69) implies A > 0 and B = 0.
To fix A, we need to consider the order 1 , τ 0 , which gives the differential equations Inserting the solution of the previous step we can solve (72). We find that in order to avoid secular behaviours, we need to fix ω 1,0 = 0 and A = 2.
By repeating this procedure for higher orders of and τ , we can compute the matrix ω i,j and the corresponding solutions. For instance, up to order 5 and τ 6 , we get The first important remark here is that the coefficients of the first row (corresponding to fixing i = 0 and taking j = 0, 1, 2, 3 in equation (66)) are exactly the same as for the approximation of the frequency obtained by using the modified Hamiltonian (cf. equation (56)), which shows a remarkable consistency between the two methods. Moreover, equation (66), with the coefficients ω i,j given in (73), allows us to extend the analytical analysis for the frequency and period of the limit cycle to the case = 0. In Figure 7 we compare the analytical results thus obtained with the numerical results from Section 4.1.2. Clearly the match is very accurate, as the curves are almost indistinguishable, even for very large values of the nonlinear coupling and of the time step τ .

1 (stiff regime)
This is allegedly the most difficult regime to study, because is large and therefore the nonlinear terms are important. Typically we must rely on the numerical results. However, we can give an argument for a reasonable measure of the distance between the simulated numerical dynamics and the original one: from a direct inspection of the modified Hamiltonian (see e.g. (58)), one can see directly that for any truncation up to order 2i in τ , we get a polynomial of the same order in . We formalise this observation in the following result. Proof. It is can be proved that (see [9]), each correction ∆H 2i in (32) is the result of taking 2i + 1 nested Jacobi brackets. Since the Jacobi bracket is anti-symmetric, we may have at most 2i equal terms inside the nested brackets. Considering that in the splitting (34) only the A term depends (linearly) on , and given the linearity of the Jacobi bracket, the greatest power in is given by the term {A, {A, {· · · , {A, P } η · · · } η } η } η , with P being either B or C. We conclude that the maximum degree in of ∆H 2i is just 2i.
From Proposition 4.4 it follows that the largest power in and τ in each correction ∆H 2i in the modified Hamiltonian is of the form ( τ ) 2i . Recalling that 1 in this case, one can expect that to keep the sum (32) under control, special attention should be given to the size of the product τ . This agrees with the results in Section 4.1.3, where we observed that the limit cycle presents a noticeable deformation for values of = 50, 100 and τ = 0.01, or = 100 and τ = 0.005, that is, when τ = 0.5, 1.

Geometric numerical integration of forced Liénard systems
To emphasise the applicability of contact integrators to general Liénard systems, we will now present a brief numerical application of contact integrators to Liénard systems with time-dependent forcing. As usual, we take the van der Pol oscillator as our benchmark example, and study this system under the influence of a forcing term that is known to give rise to chaotic behavior [23,25].
We stress that this section is meant as an example of possible further applications and the results presented here are by no means meant to be exhaustive analyses or comparisons with the previous literature. Moreover, we will focus on the numerical aspects and omit the analytical treatment of the modified Hamiltonians: since the computations are analogue to what we have already presented for the unforced van der Pol oscillator, we believe that adding them here would unnecessarily complicate the paper.
In the simulations that follow we proceed in analogy to [25]. We test the 2nd order contact integrator S 2 (τ )(CBABC) and two different 6th order integrators: S e 6 (τ )(CBABC), with exact coefficients, and S a 6 (τ )(CBABC), with approximate coefficients taken from family A in Table 1. (these are the integrators that have showed the best performance, cf. Remark 3.1). All the comparisons are made with respect to the LSODA solver provided by SciPy [31] with a relative accuracy parameter of 10 −13 and absolute accuracy parameter of 10 −15 .

The forced van der Pol oscillator
Following [23,25], we consider a forced van der Pol oscillator of the following form here A is the amplitude of the forcing and ω its frequency. Extending (20) to a time-dependent contact Hamiltonian, we observe that the equation above can be recovered from H(q, p, s, t) = ps − (1 − q 2 )s + q − A cos(ωt).
Indeed, on the (q, s) plane, the corresponding contact Hamiltonian system reduces to The nontrivial behaviour of this example is well known [23]: e.g. for the couplings A = = 5 one can show that the system undergoes a bifurcation cascade from a regular attractor (ω = 2.457) to a chaotic one (for ω = 2.463).
In the numerical experiments, we propagate the system until t = 500 and, unless differently specified, the time step is τ = 0.02.
As one can see in Figure 8, even though we are dealing with a stiff problem, the method is capable of capturing the attractor even for large value of the time step and  . The green dots correspond to the 2nd order integrator and the orange dots to a 6th order approximate integrator (CBABC) with the coefficients taken from family A in Table 1. Left: regular attractor. Right: strange attractor. From top to bottom the time step is decreasing. The inset plots contain the corresponding trajectory computed with LSODA. it is plotted separately because, besides the first row, it is virtually indistinguishable from the one obtained with the 6th order integrator.
long integration intervals, rapidly converging to the correct solution as the time-step decreases.
This system in the chaotic regime, ω = 2.463, was also the example used to analyze the performances of the modified leapfrog methods introduced in [25]. Even though the numerical test in [25] uses a 6th order integrator, we will still include a test for our second order integrator.
Even though both integrators are geometric in nature, explicit and with fixed timestep, the ones introduced in this paper present two main differences from those in [25]: they are based on contact geometry instead of symplectic one and they require the integration of only three variables (one of which is the time) instead of six.
In Figure 9 we show the trajectories computed by the aforementioned integrators. As one can see by comparing Figure 10 and [25, Figure 4], despite the simplicity of the contact methods, their performance is comparable to the ones presented in [25], with the approximate integrator performing better than the exact one: they give results comparable to an established differential equation solver, LSODA, with less computational work: for these simulations the amounts of vector field evaluations of LSODA is 1.4 · 10 6 , while our second order integrator requires 0.1 · 10 6 evaluations and the sixth order one 0.3 · 10 6 . These results can also be contrasted with the amount of evaluations for the corresponding algorithms in [25], which are 0.7 · 10 6 (Method 1) and 1.3 · 10 6 (Method 2). By avoiding the phase space extension we obtained two concrete advantages: we have reduced both the computational cost of the integrator and the number of possible combinations of the splitting maps.

Conclusions
In this work we have proposed a novel approach to the geometric numerical integration of an important class of nonlinear dynamical systems, that is, Liénard systems. Such systems are planar systems having a limit cycle, and therefore they cannot be Hamiltonian in the symplectic sense in their original variables. As a minimal extension, we have considered Liénard systems as 2-dimensional projections of contact Hamiltonian systems in three dimensions. This Hamiltonisation enables us to use the contact splitting integrators recently introduced in [9] and therefore to derive a new class of geometric numerical integrators for Liénard systems. We have used the paradigmatic example of the van der Pol oscillator to show that such formulation can be beneficial both for obtaining accurate numerical integrations of the dynamics at relatively small computational cost, and for deriving complementary analytical results, based on the use of the modified Hamiltonian and modified equations.
Although we have shown here some important results, several questions still remain to be addressed. For instance, we have not fully exploited the modified Hamiltonian and  Figure 10: Maximum absolute errors in x andẋ up to a given time for the reference contact integrators compared to the LSODA method along the orbit in Figure 9.
modified equations in the stiff case; we have not considered further theoretical properties related to the existence of a Hamiltonian structure, such as e.g. the preservation of volumes in the 3-dimensional manifold, or the associated Lagrangian structure. In this context, we remark that the approach investigated here is based on the simplest possible Hamiltonisation of Liénard systems by means of contact Hamiltonian systems, which is obtained by a Hamiltonian that is linear (hence singular) in the momenta. Therefore to derive an associated Lagrangian structure one would have to use the algorithm for singular contact Hamiltonian systems developed in [14]. From the numerical perspective, this could open the door to the use of contact variational integrators [28,30]. Moreover, other (contact) Hamiltonisations of Liénard and spiking systems might be possible, perhaps using non-standard contact structures, and therefore future work should also focus on alternative constructions.