Next Article in Journal
Analytical Solution of Laterally Loaded Free-Head Long Piles in Elasto-Plastic Cohesive Soils
Next Article in Special Issue
Warped Product Submanifolds in Locally Golden Riemannian Manifolds with a Slant Factor
Previous Article in Journal
DMO-QPSO: A Multi-Objective Quantum-Behaved Particle Swarm Optimization Algorithm Based on Decomposition with Diversity Control
Previous Article in Special Issue
Some Conditions on Trans-Sasakian Manifolds to Be Homothetic to Sasakian Manifolds
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach

1
Bernoulli Institute for Mathematics, Computer Science and Artificial Intelligence, University of Groningen, 9747 AG Groningen, The Netherlands
2
Instituto de Investigaciones en Matemáticas Aplicadas y en Sistemas (IIMAS–UNAM), Mexico City 04510, Mexico
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(16), 1960; https://doi.org/10.3390/math9161960
Submission received: 20 July 2021 / Accepted: 10 August 2021 / Published: 16 August 2021
(This article belongs to the Special Issue Differential Geometry: Structures on Manifolds and Their Applications)

Abstract

:
Starting from a contact Hamiltonian description of Liénard 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.
MSC:
65D30; 34K28; 34A26; 34C15

1. Introduction

Liénard systems are a class of two-dimensional nonlinear dynamical systems that exhibit a stable limit cycle. Among them the most famous is the van der Pol oscillator [1,2]. 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 [3,4], in the general case such a 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 to circumvent this problem. In [5] the authors suggested enlarging the phase–space to a four-dimensional manifold and define a particularly simple Hamiltonian system in this enlarged space so that the two-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 [6], the classical Bateman trick for the harmonic oscillator was 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 four-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 [7] the authors have presented various splitting schemes for “conditionally linear systems” (these include Liénard systems) which, although not 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 intense research, especially related to low-dimensional topology [8]. In recent years, contact Hamiltonian systems have found many applications, first in the context of thermodynamics [9,10,11] and, more recently, in the context of the Hamiltonisation of several dissipative dynamical systems [12,13,14,15,16,17,18,19]. The large number of applications of contact systems that have appeared recently motivated research on geometric numerical integration [15,16,20,21]. 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 ideas in [5], 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 three-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 sparing computational resources and is of primal importance in applications to, e.g., neuronal dynamics [7]. 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 Section 4 and Section 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 [22].

2. A Contact Hamiltonian Formulation of Liénard Systems

2.1. A Brief Review of Liénard Systems

Liénard systems are a family of planar coupled differential equations of the form [1]
x ˙ = y F ( x ) y ˙ = g ( x ) ,
where F ( x ) = f ( x ) d x for an even function f ( x ) and g ( x ) is an odd function. Alternatively, (1) is equivalent to the second order scalar equation
x ¨ = f ( x ) x ˙ g ( x ) .
A third equivalent version of (1) is given by
x ˙ = y , y ˙ = g ( x ) f ( x ) y .
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
x ¨ = ϵ ( 1 x 2 ) x ˙ x ,
and can be equivalently rewritten in the form (3) as
x ˙ = y , y ˙ = x + ϵ ( 1 x 2 ) y ,
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 [23].
Theorem 1
(Liénard’s Theorem). Under the conditions
  • F , g C 1 ( R ) ,
  • x g ( x ) > 0 if x 0 ,
  • F ( 0 ) = 0 and f ( 0 ) < 0 ,
  • 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 [23].

2.2. A Brief Review of Contact Hamiltonian Systems

Similarly to the fact that a symplectic manifold is a 2 n -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 ( 2 n + 1 ) -dimensional manifold endowed with a 1-form η , called the contact form, that is non-degenerate, which means
η ( d η ) n 0 .
A contact version of Darboux’s theorem [24] guarantees the local existence of coordinates ( q i , p i , s ) —called Darboux coordinates—which permit to express the contact form as η = d s p i d q 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
ι X H η = H ι X H d η = d H ι R H η ,
where ι X is the interior product and R is the Reeb vector field corresponding to η [8].
In Darboux coordinates X H takes the form
X H = H p i q ˙ q i + p i H s H q i p ˙ p i + p i H p i H s ˙ s ,
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
f , g η : = ι [ X f , X g ] η .
Again, in Darboux coordinates the Jacobi bracket reads
{ g , f } η = g f s g s f + p μ g s f p μ g p μ f s + g q μ f p μ g p μ f q μ .
We refer the reader to [11,13,17,24,25] 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,
μ q i p j s r , μ ¯ q i ¯ p j ¯ s r ¯ η = μ μ ¯ ( 1 j ) r ¯ + ( j ¯ 1 ) r q i + i ¯ p j + j ¯ s r + r ¯ 1 + i j ¯ i ¯ j q i + i ¯ 1 p j + j ¯ 1 s r + r ¯ ,
where μ , μ ¯ R and i , j , r , i ¯ , j ¯ , r ¯ N .

2.3. A Contact Hamiltonian Formulation of Liénard Systems

It is well known that any dynamical system on a n-dimensional manifold Q of the form x ˙ i = X i ( x ) can be extended to a Hamiltonian system defined on the 2 n -dimensional phase–space T * Q . This can be achieved with the introduction of the conjugate momenta p ˜ i in order to define the Hamiltonian
H ( x , p ˜ ) = p ˜ i X i ( x ) .
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
H ( x , y , p ˜ 1 , p ˜ 2 ) = p ˜ 1 y p ˜ 2 ( g ( x ) + f ( x ) y ) , ( p ˜ 1 , p ˜ 2 ) = ( p ˜ x , p ˜ y ) .
In [5], such an approach was used to derive a Hamiltonisation of Liénard systems in such an extended phase–space that was then shown to be useful to perform perturbation theory. Moreover, in [26] a similar extension, but with a suitably defined new Hamiltonian that non-trivially couples the variables, was 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 ( 2 n 1 ) -dimensional projectivised cotangent bundle P T * Q , endowed with the contact structure inherited from the canonical symplectic structure of T * Q [24,27]. The procedure to perform such reduction is quite simple in this case and it is reviewed, e.g., in the recent work [10]. 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 which p ˜ 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
H ( q , p , s ) = 1 p ˜ 2 H ( x , y , p ˜ 1 , p ˜ 2 ) = p X 1 ( q , s ) X 2 ( q , s ) .
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.
Theorem 2
(Hamiltonisation of Liénard systems). Liénard systems are contact Hamiltonian systems, with a Hamiltonian of the form
H = p s + f ( q ) s + g ( q ) .
The associated contact Hamiltonian system is
q ˙ = s
s ˙ = f ( q ) s g ( q ) ,
p ˙ = p 2 f ( q ) p f ( q ) s g ( q ) .
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
f ( x ) = ϵ ( 1 x 2 ) , g ( x ) = x .
Consequently the contact Hamiltonian in this case reads
H = p s ϵ ( 1 q 2 ) s + q ,
and the corresponding contact Hamiltonian systems is
q ˙ = s s ˙ = ϵ ( 1 q 2 ) s q p ˙ = 1 p 2 + ϵ ( 1 q 2 ) p 2 q s .
As expected, from the first two equations we recover the original van der Pol Equation (4).
Remark 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 d s d q 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 [28].
Remark 2.
The reduction procedure that led us to (14) is not unique. Indeed, we could have selected the connected component in which p ˜ 1 0 and set ( q = y , s = x , p = p ˜ 2 p ˜ 1 ) . The corresponding contact Hamiltonian for Liénard systems is:
K ( q , p , s ) = 1 p ˜ 1 H ( x , y , p ˜ 1 , p ˜ 2 ) = p X 2 ( q , s ) X 1 ( q , s )
= p ( f ( s ) q + g ( s ) ) q .
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
q ˙ = f ( s ) q g ( s ) s ˙ = q p ˙ = 1 + p f ( s ) + p p q f ( s ) + g ( s ) ,
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 overall.

3. Geometric Numerical Integration of Liénard Systems

3.1. Contact Splitting Integrators

Contact splitting integrators are a class of geometric integrators recently introduced in the context of celestial mechanics [15]. 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
H ( q i , p i , s ) = j = 1 N h j ( q i , p i , s ) .
Then, the Hamiltonian vector field associated with H is separable as well
X H = j = 1 N X h j .
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 1
(Contact splitting integrators). In the hypotheses above, let e t X h j denote the map given by the time-t exact flow of each vector field X h j , for j = 1 , , N . Then
S 2 ( τ ) = e τ 2 X h 1 e τ 2 X h 2 e τ X h N e τ 2 X h 2 e τ 2 X h 1
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 [29], 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.
Proposition 2
(Integrator with exact coefficients). If S 2 n ( τ ) is an integrator of order 2 n , then the map
S 2 n + 2 ( τ ) = S 2 n ( z 1 τ ) S 2 n ( z 0 τ ) S 2 n ( z 1 τ ) ,
with z 0 and z 1 given by
z 0 ( n ) = 2 1 2 n + 1 2 2 1 2 n + 1 , z 1 ( n ) = 1 2 2 1 2 n + 1 ;
is an integrator of order 2 n + 2 .
Proposition 3
(Integrator with approximated coefficients). There exist m N and a set of real coefficients { w j } j = 0 m such that the map
S ( m ) ( τ ) = S 2 ( w m τ ) S 2 ( w m 1 τ ) S 2 ( w 0 τ ) S 2 ( w m 1 τ ) S 2 ( w m τ ) ,
is an integrator of order 2 n .
In Table 1 we list the values of the approximated coefficients { w j } j = 0 m for three different sixth order integrators, labelled as A, B and C. Note that w 0 : = 1 2 j = 1 m w i .
Remark 3.
The splitting integrator with approximate coefficients labelled as A is the better performer among the approximate splitting integrators of sixth order presented here. This can be related to the fact that its largest coefficient is the smallest among the approximate integrators.

3.2. 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 [15] for further details on the derivation of the modified Hamiltonian in the contact case). Indeed, for an integrator of order 2 n multiple applications of the BCH formula give [15]
S 2 n ( τ ) = exp τ X H + i = n τ 2 i X 2 i ,
where all the corrections X 2 i are Hamiltonian vector fields. Therefore, S 2 n ( τ ) 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
H m o d , 2 n ( q a , p a , s ; τ ) = H ( q a , p a , s ) + i = n τ 2 i Δ H 2 i ( q a , p a , s ) ,
where the subscript 2 n in H m o d , 2 n denotes the fact that it is associated with an integrator of order 2 n , and Δ H 2 i are the Hamiltonian functions associated with the Hamiltonian vector fields X 2 i , that is,
Δ H 2 i ( q a , p a , s ) = ι X 2 i η .
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 2 n ( τ ) . 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.

3.3. 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
H = p s C + f ( q ) s A + g ( q ) B ,
and the consequent identification of the corresponding vector fields
X A = p f ( q ) + s f ( q ) p s f ( q ) s ,
X B = g ( q ) p g ( q ) s ,
X C = s q p 2 p .
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
e τ X A q i + 1 = q i p i + 1 = e τ f ( q i ) ( p i f ( q i ) s i τ ) s i + 1 = e τ f ( q i ) s i e τ X B q i + 1 = q i p i + 1 = g ( q i ) τ + p i s i + 1 = g ( q i ) τ + s i e τ X C q i + 1 = q i + s i τ p i + 1 = p i 1 + p i τ s i + 1 = s i
Example 3
(The van der Pol oscillator yet again). Applying the above splitting to the Hamiltonian (20) we obtain
H = p s C ϵ ( 1 q 2 ) s A + q B ,
and the corresponding time-τ flow maps are
e τ X A q i + 1 = q i p i + 1 = e q i 2 1 τ ϵ ( p i 2 q i s i τ ϵ ) s i + 1 = e τ ϵ ( 1 q i 2 ) s i e τ X B q i + 1 = q i p i + 1 = p i τ s i + 1 = s i τ q i e τ X C q i + 1 = q i + s i τ p i + 1 = p i 1 + p i τ s i + 1 = s i
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 ( τ ) ( C B A B C ) to indicate that we are using the second order integrator obtained using the splitting (27) of the maps (40) composed in the order indicated in parentheses.
Remark 4.
As we show in Proposition 4, the equations of motion of q and s are necessarily independent of p. In fact one could have just considered the first two rows of (21) and used the corresponding splitting integrator. The integrator itself would be identical to the restriction of (40) to the ( q , s ) plane.
By going through the route of the contact Hamiltonian, however, we show that there is more to this scheme in certain cases: It produces geometric integrators, guaranteed to preserve an underlying contact structure. Our belief is that this hidden structure is the reason the integrator is so performant in terms of the preservation of the limit cycle, even in the stiff regime.

3.4. A Remark on Variational Integrators

We can also think of Liénard systems as forced Lagrangian systems. The discrete version of the Lagrange–D’Alembert principle can then be used to construct geometric integrators for the systems without the need to introduce an ancillary variable. The theory of variational integration for forced systems is developed in detail in [30] (Part 3), we present here only the bare minimum and the integrator that results. We will use this integrator in the following sections as an extra element of comparison for our splitting integrators.
If we denote the Lagrangian by L : T Q R and the external force by F : T Q T * Q , the Lagrange–D’Alembert principle tells us that the trajectories of the system are the solutions of
δ 0 T L ( q ( t ) , q ˙ ( t ) ) d t + 0 T F ( q ( t ) , q ˙ ( t ) ) · δ q d t = 0 .
In the case of the van der Pol oscillator, L ( q , q ˙ ) = 1 2 ( q ˙ 2 q 2 ) and F ( q , q ˙ ) = ϵ ( 1 q 2 ) q ˙ .
For the discrete Lagrange–D’Alembert principle, we seek a discrete curve { q k } k = 0 N satisfying
δ k = 0 N 1 L d ( q k , q k + 1 ) + k = 0 N 1 F d ( q k , q k + 1 ) δ q k + F d + ( q k , q k + 1 ) δ q k + 1 = 0
where L d : Q × Q R and F d ± : Q × Q T * ( Q × Q ) are the discretisations of the Lagrangian and the force, respectively. For instance, when considering the midpoint rule,
L d ( q 0 , q 1 ; τ ) : = τ L q 0 + q 1 2 , q 1 q 0 τ F d ± ( q 0 , q 1 ; τ ) : = τ 2 F q 0 + q 1 2 , q 1 q 0 τ .
Applying a discrete Legendre transform with forces, one ends up with the implicit integrator
p n = D 1 L d ( q n , q n + 1 ) F d ( q n , q n + 1 ) , p n + 1 = D 2 L d ( q n , q n + 1 ) + F d + ( q n , q n + 1 ) ,
where D i denotes the derivative with respect to the ith variable.
In all cases, however, the explicit nature of the contact integrator makes it orders of magnitude faster than the variational one, as shown in Table 2.

4. Geometric Numerical Integration of the van der Pol Oscillator: Numerical vs. Analytical Results

4.1. 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 [7]: 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.

4.1.1. ϵ = 0 (Harmonic Oscillator)

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 τ = 0.785 . By increasing the time step, the typical circular orbit of the harmonic oscillator becomes more elliptic, and the period changes. In this regime, in which we are effectively integrating a Hamiltonian system, we can see the extra stability of the variational integrator compared to a splitting integrator.
Let us focus for a moment only on the contact integrator. 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 ] .

4.1.2. ϵ 1 and ϵ 1 (Non-Stiff Regime)

In Figure 3 we show the persistence of the limit cycle for different values of ϵ { 0.1 , 0.5 , 1 , 2 , 4 } and increasing values of τ .
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 , 10 , 000 ] , is an evidence of the stability of the integrator.
Moreover, in this case, the dependence of the period and the frequency of the limit cycle with respect to the time step shown in Figure 4 is very similar to that of the harmonic oscillator.
According to adiabatic perturbation theory [31] (Chapter XII), we can expect the variational integrator obtained in Section 3.4 to perform relatively well in the integration of Liénard systems with small coupling. However, already in this regime, we start observing a difference in stability with respect to the variational integrator. In Figure 5, where we replicate the simulation of Figure 3 with the variational integrator of Section 3.4, we can clearly see that the increase of the time step τ corresponds to a progressive disappearance of the limit cycle, also for small values of ϵ .

4.1.3. ϵ 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 [1,7]
q ̲ = q s ̲ = q q 3 3 s ϵ .
This change of variables transforms the dynamics into
q ̲ ˙ = ϵ q ̲ q ̲ 3 3 s ̲ s ̲ ˙ = q ̲ / ϵ .
and enables a nice geometric description of the limit cycle. Indeed, the q ̲ nullcline, which is the locus of points such that q ̲ ˙ = 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 6 shows the cubic nullcline and the numerically simulated attractor for ϵ { 25 , 50 , 100 } and for different values of the time step. As one can observe, the limit cycle is preserved also for large values of the nonlinear coupling, although it suffers from a distortion for larger values of τ : This is especially clear in the first picture of the last row of plots of Figure 6, corresponding to ϵ = 100 and τ = 0.01 .
When we replicate the previous figure also for the variational integrator, see Figure 7, we see again that there is a sudden disappearance of the limit cycle also for values of τ which are relatively small compared to ϵ . This gets more and more pronounced by increasing the coupling constant and seems to indicate a complex dependence of the error on the product τ ϵ . In contrast, we prove in Proposition 7 that the contact integrator of order 2 i has an error of order τ 2 i ( 1 + ϵ 2 i ) . Thus we can expect, provided τ ϵ < 1 , to preserve the qualitative features of the dynamics, in agreement with the results of Figure 6.

4.2. 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 m o d , i.e., that it is linear in p, as it is the original Hamiltonian (20). This is the content of the next result.
Proposition 4.
For any contact splitting integrator, the corresponding modified Hamiltonian H m o d is linear in p. It follows that the modified equations for q and s are independent of p.
Proof. 
We prove first the second part: The claim is that if H m o d 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 m o d is indeed linear in p: Considering the splitting in (39), we have that A = ϵ 1 q 2 s , B = q , and C = p s , 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. Indeed, H m o d is linear in p, because otherwise in the modified equations we would have q ˙ = 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.
Proposition 5.
Restricted to the plane ( q , s ) , the integrator S 2 ( τ ) ( C B A B C ) has a unique fixed point at ( 0 , 0 ) which is unstable. Furthermore, both eigenvalues λ 1 , 2 of the Jacobian of the mapping ( q i , s i ) ( q i + 1 , s i + 1 ) satisfy | λ 1 , 2 | > 1 for all τ > 0 and ϵ > 0 .
Proof. 
The proof is based on writing explicitly the action of the integrator on an initial condition, that is, we apply e τ / 2 X C e τ / 2 X B e τ X A e τ / 2 X B e τ / 2 X C to ( q i , s i ) , to obtain
q i + 1 = q i + τ 2 s i + τ 2 s i + 1 s i + 1 = e ϵ τ ( 1 ( q i + τ 2 s i ) 2 ) s i τ 2 q i + s i τ 2 τ 2 q i + s i τ 2 .
Now when we impose the condition for a fixed point
q i + 1 = q i s i + 1 = s i .
Using the second equation in (45) into the first equation in (44) we obtain
q i + 1 = q i + τ s i = q i ,
which is true if and only if s i = 0 .
Next, we substitute s i = 0 = s i + 1 into the second equation in (44) and we obtain
0 = s i + 1 = τ 2 q i e ϵ τ ( 1 ( q i ) 2 ) + 1 ,
which is true if and only if q i = 0 .
To prove that ( 0 , 0 ) is unstable, we compute the Jacobian of the map (44) 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 (44) at ( 0 , 0 ) are
λ 1 , 2 = 1 4 α ± β , α : = ( 2 τ 2 ) ( e ϵ τ + 1 ) , β : = α 2 16 e ϵ τ .
Depending on the sign of β we have two cases: The eigenvalues are both real or they are complex conjugates.
Case (I)
λ 1 , 2 C : The eigenvalues are complex conjugates, therefore | λ 1 | = | λ 2 | . Since det J = λ 1 λ 2 = e ϵ τ , we have | λ 1 | = | λ 2 | = e ϵ τ 2 > 1 .
Case (II)
λ 1 , 2 R : This happens when β 0 , that is
α 4 e ϵ τ 2 .
The fact that λ 1 > 1 follows from λ 1 = 1 4 α + β 1 4 α ( 47 ) e ϵ τ 2 > 1 .
Let us now focus on λ 2 . Notice that since λ 1 λ 2 > 0 and λ 1 > 0 , we necessarily have that λ 2 > 0 . Therefore, it suffices to prove that
λ 2 = 1 4 α β > 1 .
By repeatedly rearranging (48) and observing that (47) implies α 4 > 0 , we obtain that (48) is equivalent to the following inequalities
α β > 4
( α 4 ) 2 > β
α 2 8 α + 16 > α 2 16 e ϵ τ
16 8 ( 2 τ 2 ) ( e ϵ τ + 1 ) > 0 .
Since ( 2 τ 2 ) < 2 , (52) is always true, proving (48).
   □
In what follows we split the analysis into three different cases depending on the value of ϵ , as we did in the Section 4.1.

4.2.1. ϵ = 0 (Harmonic Oscillator)

In this case, we have a harmonic oscillator, for which each non-trivial 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 m o d , 2 = p s F ( τ ) + q G ( τ )
where
F ( τ ) = 1 τ 2 12 τ 4 120 τ 6 840 τ 8 5040 + O ( τ 10 ) ,
G ( τ ) = 1 + τ 2 6 + τ 4 30 + τ 6 140 + τ 8 630 + O ( τ 10 ) .
The corresponding modified system is thus
q ˙ ( t ) = s F ( τ ) s ˙ ( t ) = q G ( τ ) p ˙ ( t ) = p 2 F ( τ ) G ( τ )
which is again exactly solvable (recall that τ is fixed), and the solution in q and s is a harmonic oscillator with frequency
ω ( τ ) = F ( τ ) G ( τ ) = 1 + τ 2 24 + 3 τ 4 640 + 5 τ 6 7168 + 35 τ 8 294 , 912 + O ( τ 10 ) .
In Figure 8 we compare (57) 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 eighth order in τ and the numerical results.

4.2.2. ϵ 1 (Non-Stiff Regime)

This regime can be studied using perturbation theory and therefore there are many results (see, e.g., [32,33]). We study the persistence of the limit cycle for the contact splitting integrators in a way similar to [7], that means, we use the modified equations in order to provide some estimations on the amplitude and period of the limit cycle.
Proposition 6.
For any contact splitting integrator of order 2 n 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 ( τ 2 n ) . Moreover, the approximate radius of the limit cycle for the S 2 ( τ ) ( C B A B C ) integrator, up to order 4 in τ, is
r = 2 τ 2 4 + O τ 4 .
Proof. 
Let us consider a contact splitting integrator S 2 n ( τ ) of order 2 n ; using the BCH formula (see Section 3.2) we can argue that the modified Hamiltonian whose time- τ flow is given by S 2 n ( τ ) is of the form
H m o d , 2 n = p s ϵ 1 q 2 s + q + τ 2 n Δ H 2 n ( q , p , s ) + O ( τ 2 n + 2 ) .
Thus the modified equations read
q ˙ = s + τ 2 n Δ H 2 n p + O ( τ 2 n + 2 ) s ˙ = q ϵ ( 1 q 2 ) s + τ 2 n Δ H 2 n q p Δ H 2 n s + O ( τ 2 n + 2 ) p ˙ = 1 p 2 + ϵ ( 1 q 2 ) p 2 q s + τ 2 n p Δ H 2 n p Δ H 2 n + O ( τ 2 n + 2 ) .
We know from Proposition 4 that the equations for q ˙ and s ˙ are independent of p, and from Proposition 5 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 r ˙ reads
r ˙ = ϵ r sin 2 ( θ ) 1 r 2 cos 2 ( θ ) + τ 2 n R 2 n ( r , θ ) + O ( τ 2 n + 2 ) ,
Since the modified Hamiltonian is by construction a polynomial in the variables q and s, the dependence on θ of R 2 n is only through sums and products of trigonometric functions. In particular, this implies that the averaged dynamics of r ˙ obtained by the integration along a period has the form
1 2 π 0 2 π r ˙ d θ = 1 8 r r 2 4 ϵ + O ( τ 2 n ) .
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 proposition.
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.
To prove the second part of the statement, we concentrate on the integrator S 2 ( τ ) ( C B A B C ) (since this is the integrator that has been used throughout the simulations in the paper). The corresponding modified Hamiltonian, in this case, is
H m o d , 2 = p s + ϵ q 2 1 s + q + τ 2 12 ( q 2 1 ϵ 2 q p q s + q 2 + 4 s 2 1 p s + ϵ p q q 2 2 s 2 1 s 7 q 2 + s 2 + 1 p s + 2 q ) + O ( τ 3 ) ,
leading to the following modified equations for q ˙ and s ˙
q ˙ = s + τ 2 12 q ϵ q 2 2 s 2 1 + q 2 1 2 s ϵ 2 s s ˙ = q ϵ ( 1 q 2 ) s + τ 2 12 q q 2 1 ϵ 2 q 2 + 4 s 2 1 + s ϵ 7 q 2 + s 2 + 1 2 q ,
and to the radial equation
r ˙ = ϵ r sin 2 ( θ ) 1 r 2 cos 2 ( θ ) + τ 2 12 r ( 3 sin ( θ ) cos ( θ ) 4 r 2 ϵ 2 sin 3 ( θ ) cos ( θ ) r 2 cos 2 ( θ ) 1 + r 2 ϵ sin 4 ( θ ) + ϵ cos 2 ( θ ) r 2 cos 2 ( θ ) 1 + ϵ sin 2 ( θ ) 1 9 r 2 cos 2 ( θ ) ) .
An explicit computation then gives
1 2 π 0 2 π r ˙ d θ = 1 32 r ϵ r 2 τ 2 + 4 16 + O ( τ 4 ) ,
leading to the claimed radius r = 2 τ 2 4 + O τ 4 .   □
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., [33]). The first step consists in the time reparametrisation t = ω t , which leads to the differential equation
ω q = X H m o d q ω s = X H m o d s , .
where the derivatives are now expressed in terms of t , instead of t, and, as usual, we omit the decoupled equation for p ˙ . Noticing that the modified Hamiltonian vector field depends on the two parameters ϵ and τ , we suppose, in analogy to the traditional approach [33], that all the terms appearing in the equations can be expanded in Taylor series with respect to such parameters as follows
ω ( ϵ , τ ) = i = j = 0 + ω i , j ϵ i τ 2 j ,
q ( t , ϵ , τ ) = i = j = 0 + q i , j ( t ) ϵ i τ 2 j ,
s ( t , ϵ , τ ) = i = j = 0 + s i , j ( t ) ϵ i τ 2 j .
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 [33] and assume that
q ( 0 , ϵ , τ ) = 0 , q ( 0 , ϵ , τ ) > 0 .
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 )
whose solution is
q 0 , 0 ( t ) = A cos t ω 0 , 0 + B sin t ω 0 , 0 , s 0 , 0 ( t ) = A sin t ω 0 , 0 + B cos t ω 0 , 0 .
Since we want q 0 , 0 ( t ) and s 0 , 0 ( t ) to have period 2 π , this fixes ω 0 , 0 = 1 , while condition (70) implies A > 0 and B = 0 .
To fix A, we need to consider the order ϵ 1 , τ 0 , which gives the differential equations
ω 1 , 0 q 0 , 0 ( t ) + q 1 , 0 ( t ) = s 1 , 0 ( t ) , ω 1 , 0 s 0 , 0 ( t ) + s 1 , 0 ( t ) = q 1 , 0 ( t ) + ( 1 q 0 , 0 2 ( t ) ) s 0 , 0 ( t ) .
Inserting the solution of the previous step we can solve (73). 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
ω i , j = 1 1 24 3 640 5 7168 0 0 0 0 1 16 27 128 149 2048 559 16 , 384 0 0 0 0 17 3072 781 73 , 728 339 , 041 3 , 538 , 944 4 , 695 , 149 84 , 934 , 656 0 0 0 0 .
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 (67)) are exactly the same as for the approximation of the frequency obtained by using the modified Hamiltonian (cf. Equation (57)), which shows a remarkable consistency between the two methods. Moreover, Equation (67), with the coefficients ω i , j given in (74), allows us to extend the analytical analysis for the frequency and period of the limit cycle to the case ϵ 0 . In Figure 9 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 τ .

4.2.3. ϵ 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: An analysis of the modified Hamiltonian (see, e.g., (59)) allows now to control the error of the contact integrator of order 2 i asymptotically with τ 2 i ( 1 + ϵ 2 i ) . Provided that τ ϵ < 1 , this leads to the preservation of the qualitative features of the dynamics, in agreement with what is shown in Figure 6.
The result follows from the structure of the modified Hamiltonian itself and its explicit construction via the Jacobi brackets. It is important here to stress that this analysis holds also for large values of ϵ and not just for the adiabatic case ϵ 1 .
Proposition 7.
For any contact splitting integrator of order 2 n based on the maps (40), the truncation at order 2 m (in τ), m n , of the modified Hamiltonian differs from the original Hamiltonian by a polynomial in ϵ of order at most 2 m . Moreover, any term of order ϵ j is of the form ϵ j τ k , j k and 2 n k 2 m .
Proof. 
Let n i m , it can be proven that (see [15]), each correction Δ H 2 i in (32) is the result of taking 2 i + 1 nested Jacobi brackets. Since the Jacobi bracket is anti-symmetric, we may have at most 2 i 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 2 i is at most 2 i and the result follows from the structure of (32).   □
From Proposition 7 it follows that the each correction τ 2 i Δ H 2 i in the modified Hamiltonian contributes a term of order τ 2 i ( 1 + ϵ 2 i ) . 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 .

5. 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 behaviour [26,34].
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 [26]. We test the second order contact integrator S 2 ( τ ) ( C B A B C ) and two different sixth order integrators: S 6 e ( τ ) ( C B A B C ) , with exact coefficients, and S 6 a ( τ ) ( C B A B C ) , with approximate coefficients taken from family A in Table 1 (these are the integrators that have showed the best performance, cf. Remark 3). All the comparisons are made with respect to the LSODA solver [35] provided by SciPy [36], a robust adaptive solver with automatic and dynamic selection of the applied stiff or nonstiff methods, with a relative accuracy parameter of 10 13 and an absolute accuracy parameter of 10 15 .
Remark 5.
At this point, one may wonder why we are not also making a comparison with a variational integrator of the type described in Section 3.4. As a matter of fact, variational integrators with time-dependent forcing can be developed using [30] (Parts 3 and 4). However, the introduction of time into the equation comes with the introduction of an extra pair of equations involving time’s dual variable, the energy. As a direct consequence, one ends up with an implicit integrator in position, momentum, time and energy. In general, a fixed time step will not satisfy the extra equations involving the energy and the corresponding integrator is not going to be geometric [30] (Section 4.3.3). A naïve simulation using a fixed–step variational integrator with time-dependent forcing will immediately show this fact: Unless the time step is minuscule, the integrated trajectory will miss some evident qualitative features of the dynamics, see Figure 11. An additional simulation, reproducing the equivalent of Figure 10, can be found in the accessory simulation code [22].
The contact splitting integrator benefits from being an explicit fixed–time integrator even in presence of time-dependent terms. This results in a simpler implementation and much faster computation times, even when the number of evaluations of the vector fields is virtually the same as the variational one.

The Forced van der Pol Oscillator

Following [26,34], we consider a forced van der Pol oscillator of the following form,
x ¨ = ϵ ( 1 x 2 ) x ˙ x + A cos ( ω t ) ,
where 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 ) = p s ϵ ( 1 q 2 ) s + q A cos ( ω t ) .
Indeed, on the ( q , s ) plane, the corresponding contact Hamiltonian system reduces to
q ˙ = s s ˙ = ϵ ( 1 q 2 ) s q + A cos ( ω t ) .
The nontrivial behaviour of this example is well known [34], 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 10, 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 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 analyse the performances of the modified leapfrog methods introduced in [26]. While the numerical test in [26] uses a sixth order integrator, we will still include a test for our second order integrator.
While both integrators are geometric in nature, explicit and with fixed time step, the ones introduced in this paper present two main differences from those in [26]: 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 11 we show the trajectories computed by the aforementioned integrators. As one can see by comparing Figure 12 and Figure 4 in [26], despite the simplicity of the contact methods, their performance is comparable to the ones presented in [26], 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 [26], 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.
As we already discussed in Remark 5, a time-dependent variational integrator requires an additional implicit equation in time and energy. As one can observe in Figure 11, the variational integrator obtained by imposing a fixed time step is unable to reproduce fine details of the system unless the time step is extremely small. In addition to this, the implicit nature of the integrator renders it extremely slower than the contact integrator, see also Table 2.

6. 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 two-dimensional projections of contact Hamiltonian systems in three dimensions. This Hamiltonisation enables us to use the contact splitting integrators recently introduced in [15] 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.
While we have shown here some interesting results, several questions still remain to be addressed. For instance, we have not fully exploited the modified Hamiltonian and 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 three-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 [37]. From the numerical perspective, this could open the door to the use of contact variational integrators [16,20]. 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.

Author Contributions

Software, F.Z. and M.S.; Supervision, A.B. and M.S.; Writing-original draft, F.Z., A.B. and M.S.; Writing-review and editing, F.Z., A.B. and M.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by the second author’s starter grant and by NWO Visitor Travel Grant 040.11.698 that sponsored the visit of AB to the Bernoulli Institute. M. Seri research is partially supported by the NWO co-fund grant 613.009.10.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

The authors would like to thank Qihuai Liu, Arjan Van der Schaft and Mats Vermeeren for multiple interesting discussions and useful comments, and Edoardo Zadra for providing emergency computational facilities during the pandemic.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Liénard, A. Etude des oscillations entretenues. Revue Générale L’électricité 1996, 23, 901–912. [Google Scholar]
  2. Van der Pol, B. A theory of the amplitude of free and forced triode vibrations. Radio Rev. 1 1920, 701–710, 754–762. [Google Scholar]
  3. Nucci, M.C.; Tamizhmani, K.M. Lagrangians for Dissipative Nonlinear Oscillators: The Method of Jacobi Last Multiplier. J. Nonlinear Math. Phys. 2021, 17, 167. [Google Scholar] [CrossRef] [Green Version]
  4. Cariñena, J.F.; Guha, P. Nonstandard Hamiltonian structures of the Liénard equation and contact geometry. Int. J. Geom. Methods Mod. Phys. 2019, 16, 1940001. [Google Scholar] [CrossRef]
  5. Choi, J.S.; Tapley, B.D. An extended canonical perturbation method. Celest. Mech. 1973, 7, 77–90. [Google Scholar] [CrossRef]
  6. Shah, T.; Chattopadhyay, R.; Vaidya, K.; Chakraborty, S. Conservative perturbation theory for nonconservative systems. Phys. Rev. E 2015, 92. [Google Scholar] [CrossRef] [Green Version]
  7. Chen, Z.; Raman, B.; Stern, A. Structure-Preserving Numerical Integrators for Hodgkin–Huxley-Type Systems. SIAM J. Sci. Comput. 2020, 42, B273–B298. [Google Scholar] [CrossRef] [Green Version]
  8. Geiges, H. An Introduction to Contact Topology; Cambridge University Press: Cambridge, UK, 2008; Volume 109. [Google Scholar]
  9. Mrugala, R.; Nulton, J.D.; Schön, J.C.; Salamon, P. Statistical approach to the geometric structure of thermodynamics. Phys. Rev. A 1990, 41, 3156–3160. [Google Scholar] [CrossRef]
  10. Van der Schaft, A.; Maschke, B. Geometry of Thermodynamic Processes. Entropy 2018, 20, 925. [Google Scholar] [CrossRef] [Green Version]
  11. Bravetti, A. Contact geometry and thermodynamics. Int. J. Geom. Methods Mod. Phys. 2019, 16, 1940003. [Google Scholar] [CrossRef]
  12. Bravetti, A.; Tapias, D. Thermostat algorithm for generating target ensembles. Phys. Rev. E 2016, 93. [Google Scholar] [CrossRef] [Green Version]
  13. Bravetti, A.; Cruz, H.; Tapias, D. Contact Hamiltonian mechanics. Ann. Phys. 2017, 376, 17–39. [Google Scholar] [CrossRef] [Green Version]
  14. Liu, Q.; Torres, P.J.; Wang, C. Contact Hamiltonian dynamics: Variational principles, invariants, completeness and periodic behavior. Ann. Phys. 2018, 395, 26–44. [Google Scholar] [CrossRef]
  15. Bravetti, A.; Seri, M.; Vermeeren, M.; Zadra, F. Numerical integration in Celestial Mechanics: A case for contact geometry. Celest. Mech. Dyn. Astron. 2020, 132. [Google Scholar] [CrossRef] [Green Version]
  16. Vermeeren, M.; Bravetti, A.; Seri, M. Contact variational integrators. J. Phys. A Math. Theor. 2019, 52, 445206. [Google Scholar] [CrossRef] [Green Version]
  17. Gaset, J.; Gràcia, X.; Muñoz-Lecanda, M.C.; Rivas, X.; Román-Roy, N. New contributions to the Hamiltonian and Lagrangian contact formalisms for dissipative mechanical systems and their symmetries. Int. J. Geom. Methods Mod. Phys. 2020, 17, 2050090. [Google Scholar] [CrossRef]
  18. Gaset, J.; Gràcia, X.; Muñoz-Lecanda, M.C.; Rivas, X.; Román-Roy, N. A contact geometry framework for field theories with dissipation. Ann. Phys. 2020, 414, 168092. [Google Scholar] [CrossRef] [Green Version]
  19. Ciaglia, F.; Cruz, H.; Marmo, G. Contact manifolds and dissipation, classical and quantum. Ann. Phys. 2018, 398, 159–179. [Google Scholar] [CrossRef] [Green Version]
  20. Simoes, A.A.; Martín de Diego, D.; Lainz Valcázar, M.; de León, M. On the Geometry of Discrete Contact Mechanics. J. Nonlinear Sci. 2021, 31. [Google Scholar] [CrossRef]
  21. Goto, S.i.; Hino, H. Fast symplectic integrator for Nesterov-type acceleration method. arXiv 2021, arXiv:2106.07620. [Google Scholar]
  22. Zadra, F.; Seri, M.; Bravetti, A. Support Code for Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach (v2.0). Zenodo. Available online: https://research.rug.nl/en/publications/support-code-for-geometric-numerical-integration-of-l%C3%ACenard-syste (accessed on 10 August 2021). [CrossRef]
  23. Perko, L. Differential Equations and Dynamical Systems; Springer: New York, NY, USA, 1991. [Google Scholar] [CrossRef]
  24. Arnol’d, V.I. Mathematical Methods of Classical Mechanics; Springer: New York, NY, USA, 2010. [Google Scholar]
  25. De León, M.; Lainz Valcázar, M. Contact Hamiltonian systems. J. Math. Phys. 2019, 60, 102902. [Google Scholar] [CrossRef]
  26. Pihajoki, P. Explicit methods in extended phase space for inseparable Hamiltonian problems. Celest. Mech. Dyn. Astron. 2014, 121, 211–231. [Google Scholar] [CrossRef] [Green Version]
  27. Blair, D.E. Riemannian Geometry of Contact and Symplectic Manifolds; Birkhäuser: Boston, MA, USA, 2010. [Google Scholar] [CrossRef]
  28. Liu, Q.; (Guilin University of Electronic Technology, Guilin, China). Personal communication, 2020.
  29. Yoshida, H. Construction of higher order symplectic integrators. Phys. Lett. A 1990, 150, 262–268. [Google Scholar] [CrossRef]
  30. Marsden, J.E.; West, M. Discrete mechanics and variational integrators. Acta Numer. 2001, 10, 357–514. [Google Scholar] [CrossRef] [Green Version]
  31. Hairer, E.; Wanner, G.; Lubich, C. Geometric Numerical Integration. Springer Ser. Comput. Math. 2002. [Google Scholar] [CrossRef]
  32. Amore, P.; Boyd, J.P.; Fernández, F.M. High order analysis of the limit cycle of the van der Pol oscillator. J. Math. Phys. 2018, 59, 012702. [Google Scholar] [CrossRef] [Green Version]
  33. Andersen, C.M.; Geer, J.F. Power Series Expansions for the Frequency and Period of the Limit Cycle of the Van Der Pol Equation. SIAM J. Appl. Math. 1982, 42, 678–693. [Google Scholar] [CrossRef]
  34. Parlitz, U.; Lauterborn, W. Period-doubling cascades and devil’s staircases of the driven van der Pol oscillator. Phys. Rev. A 1987, 36, 1428–1434. [Google Scholar] [CrossRef]
  35. Hindmarsh, A. ODEPACK. A Collection of ODE System Solvers. In Scientific Computing; Stepleman, R.S., Ed.; North-Holland: Amsterdam, The Netherlands, 1983; pp. 55–64. [Google Scholar]
  36. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods 2020, 17, 261–272. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. De León, M.; Lainz Valcázar, M. Singular Lagrangians and precontact Hamiltonian systems. Int. J. Geom. Methods Mod. Phys. 2019, 16, 1950158. [Google Scholar] [CrossRef]
Figure 1. Orbit of the van der Pol oscillator with ϵ = 0 (harmonic oscillator) with initial condition ( q 0 , p 0 , s 0 ) = ( 0 , 0 , 1 ) integrated for different values of the time step τ . The dashed blue line shows the exact solution.
Figure 1. Orbit of the van der Pol oscillator with ϵ = 0 (harmonic oscillator) with initial condition ( q 0 , p 0 , s 0 ) = ( 0 , 0 , 1 ) integrated for different values of the time step τ . The dashed blue line shows the exact solution.
Mathematics 09 01960 g001
Figure 2. van der Pol oscillator with ϵ = 0 (harmonic oscillator). Dependence of the period of the numerical solution with respect to the time step. The inset plot is a close-up of the periods for τ [ 0.001 , 0.5 ] .
Figure 2. van der Pol oscillator with ϵ = 0 (harmonic oscillator). Dependence of the period of the numerical solution with respect to the time step. The inset plot is a close-up of the periods for τ [ 0.001 , 0.5 ] .
Mathematics 09 01960 g002
Figure 3. Limit cycle of the van der Pol oscillator computed with the second-order contact integrator for values of ϵ = 0.1 (blue), 0.5 (orange), 1 (green), 2 (red), 4 (purple) and with different time steps.
Figure 3. Limit cycle of the van der Pol oscillator computed with the second-order contact integrator for values of ϵ = 0.1 (blue), 0.5 (orange), 1 (green), 2 (red), 4 (purple) and with different time steps.
Mathematics 09 01960 g003
Figure 4. Dependence of the period of the numerical solution of the van der Pol limit cycle with respect to the time step for ϵ { 0.1 , 0.5 , 0.9 } increasing from left to right.
Figure 4. Dependence of the period of the numerical solution of the van der Pol limit cycle with respect to the time step for ϵ { 0.1 , 0.5 , 0.9 } increasing from left to right.
Mathematics 09 01960 g004
Figure 5. Limit cycle of the van der Pol oscillator computed using the variational approach of Section 3.4 for values of ϵ = 0.1 (blue), 0.5 (orange), 1 (green), 2 (red), 4 (purple) and with different time steps.
Figure 5. Limit cycle of the van der Pol oscillator computed using the variational approach of Section 3.4 for values of ϵ = 0.1 (blue), 0.5 (orange), 1 (green), 2 (red), 4 (purple) and with different time steps.
Mathematics 09 01960 g005
Figure 6. Orbits for the stiff van der Pol oscillator obtained with the second-order contact integrator for different values of the coupling ϵ and of the time step τ after the Liénard transformation: With ϵ { 25 , 50 , 100 } increasing from top to bottom and τ { 10 2 , 5 × 10 3 , 10 3 , 5 × 10 4 , 10 4 } decreasing from left to right.
Figure 6. Orbits for the stiff van der Pol oscillator obtained with the second-order contact integrator for different values of the coupling ϵ and of the time step τ after the Liénard transformation: With ϵ { 25 , 50 , 100 } increasing from top to bottom and τ { 10 2 , 5 × 10 3 , 10 3 , 5 × 10 4 , 10 4 } decreasing from left to right.
Mathematics 09 01960 g006
Figure 7. Orbits for the stiff van der Pol oscillator obtained with the variational approach of Section 3.4 for different values of the coupling ϵ and of the time step τ after the Liénard transformation: With ϵ { 25 , 50 , 100 } increasing from top to bottom and τ { 10 2 , 5 × 10 3 , 10 3 , 5 × 10 4 , 10 4 } decreasing from left to right.
Figure 7. Orbits for the stiff van der Pol oscillator obtained with the variational approach of Section 3.4 for different values of the coupling ϵ and of the time step τ after the Liénard transformation: With ϵ { 25 , 50 , 100 } increasing from top to bottom and τ { 10 2 , 5 × 10 3 , 10 3 , 5 × 10 4 , 10 4 } decreasing from left to right.
Mathematics 09 01960 g007
Figure 8. Dependence of the period of the numerical solution for the harmonic oscillator ( ϵ = 0 ) with respect to the time step. The numerically estimated period is compared with the period computed from the modified equations.
Figure 8. Dependence of the period of the numerical solution for the harmonic oscillator ( ϵ = 0 ) with respect to the time step. The numerically estimated period is compared with the period computed from the modified equations.
Mathematics 09 01960 g008
Figure 9. Comparison between the numerical and analytical results (using perturbation theory) for the period of the limit cycle. Each figure is an analogue of Figure 8 for the value of ϵ indicated in the top right corner.
Figure 9. Comparison between the numerical and analytical results (using perturbation theory) for the period of the limit cycle. Each figure is an analogue of Figure 8 for the value of ϵ indicated in the top right corner.
Mathematics 09 01960 g009
Figure 10. Orbit of the forced van der Pol oscillator with ( x 0 , x ˙ 0 ) = ( 2 , 2 ) . The green dots correspond to the second order integrator and the orange dots to a sixth 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 sixth order integrator.
Figure 10. Orbit of the forced van der Pol oscillator with ( x 0 , x ˙ 0 ) = ( 2 , 2 ) . The green dots correspond to the second order integrator and the orange dots to a sixth 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 sixth order integrator.
Mathematics 09 01960 g010
Figure 11. Numerical orbits of the forced van der Pol oscillator (75) with A = μ = 5 , ω = 2.463 and ( x 0 , x ˙ 0 ) = ( 2 , 2 ) , with the reference contact integrators and LSODA.
Figure 11. Numerical orbits of the forced van der Pol oscillator (75) with A = μ = 5 , ω = 2.463 and ( x 0 , x ˙ 0 ) = ( 2 , 2 ) , with the reference contact integrators and LSODA.
Mathematics 09 01960 g011
Figure 12. Maximum absolute errors in x and x ˙ up to a given time for the reference contact integrators compared to the LSODA method along the orbit in Figure 11.
Figure 12. Maximum absolute errors in x and x ˙ up to a given time for the reference contact integrators compared to the LSODA method along the orbit in Figure 11.
Mathematics 09 01960 g012
Table 1. The coefficients w i for three different 6th order integrators.
Table 1. The coefficients w i for three different 6th order integrators.
ABC
w 0 1.315186320683906 2.37635274430774 2.3894477832436816
w 1 1.17767998417887 2.13228522200144 0.00152886228424922
w 2 0.235573213359357 0.00426068187079180 2.14403531630539
w 3 0.784513610477560 1.43984816797678 1.44778256239930
Table 2. Execution time statistics for the integration of a van der Pol oscillator with initial conditions ( q 0 , p 0 , s 0 ) = ( 2 , 0 , 0 ) , ϵ = 3.5 for t [ 0 , 1000 ] .
Table 2. Execution time statistics for the integration of a van der Pol oscillator with initial conditions ( q 0 , p 0 , s 0 ) = ( 2 , 0 , 0 ) , ϵ = 3.5 for t [ 0 , 1000 ] .
Integrator Type (Order)Mean Running
Time (ms)
Standard Deviation
(over 10 Runs)
τ = 0.02
Contact hamiltonian (2nd)729±13.2
Contact hamiltonian (6th)3940±65.9
Variational (2nd)8630±149
τ = 0.2
Contact hamiltonian (2nd)74±1.6
Contact hamiltonian (6th)404±14.6
Variational (2nd)972±15.2
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zadra, F.; Bravetti, A.; Seri, M. Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach. Mathematics 2021, 9, 1960. https://doi.org/10.3390/math9161960

AMA Style

Zadra F, Bravetti A, Seri M. Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach. Mathematics. 2021; 9(16):1960. https://doi.org/10.3390/math9161960

Chicago/Turabian Style

Zadra, Federico, Alessandro Bravetti, and Marcello Seri. 2021. "Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach" Mathematics 9, no. 16: 1960. https://doi.org/10.3390/math9161960

APA Style

Zadra, F., Bravetti, A., & Seri, M. (2021). Geometric Numerical Integration of Liénard Systems via a Contact Hamiltonian Approach. Mathematics, 9(16), 1960. https://doi.org/10.3390/math9161960

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop