Numerical integration in celestial mechanics: a case for contact geometry

09/05/2019 ∙ by Alessandro Bravetti, et al. ∙ CIMAT University of Groningen Berlin Institute of Technology (Technische Universität Berlin) 0

Several dynamical systems of interest in celestial mechanics can be written in the form q̈ + ∂ V(q,t)/∂ q+f(t)q̇=0 . i=1,...,n . For instance, the modified Kepler problem, the spin--orbit model and the Lane--Emden equation all belong to this class. In this work we start an investigation of these models from the point of view of contact geometry. In particular we focus on the (contact) Hamiltonisation of these models and on the construction of the corresponding geometric integrators.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Geometric methods in the study of dynamical systems have proven to be useful both for analytical and for numerical investigations [hairer2006geometric, marsden2001discrete, sanz1992symplectic]

. When the system has a Hamiltonian form, it is possible to exploit the geometric structure of the problem and provide powerful tools to study and classify the dynamics. One of the central properties of Hamiltonian systems is the conservation of energy. Therefore, it is not surprising that finding a Hamiltonian formulation for second order differential equations with non–conservative forces is considered a very hard, if not hopeless, problem.

In this work we consider the important class of systems


This class arises from the Newtonian mechanics of systems with time–varying non–conservative forces, and includes systems of primary interest in celestial mechanics such as the modified Kepler problem, the spin–orbit model and the Lane–Emden equation, which will be analysed in more detail below.

Our first crucial remark is that all systems of the form (1) can be given a Hamiltonian formulation in the context of contact geometry: Section 2 will be devoted to review contact Hamiltonian systems and make the previous statement precise.

We also argue that this geometric reformulation of equation (1) can be useful in order to study the dynamics, both from an analytical and from a numerical perspective. Referring the analytical study to future investigations, we start in this paper a thorough investigation of numerical methods for the analysis of such systems based on contact geometry. In particular, in Section 3 we develop both Lagrangian and Hamiltonian integrators for contact systems, and then in Section 4 we consider the modified Kepler problem, the spin–orbit model and Lane–Emden equation, and show with numerical tests that the contact perspective can improve over previously presented methods to study the dynamics numerically.

Finally, in Section 5 we provide a summary of results and discuss future directions.


The authors would like to thank the Bernoulli Institute for the hospitality. This research was partially supported by MS’s starter grant and NWO Visitor Travel Grant 040.11.698 that sponsored the visit of AB at the Bernoulli Institute. MS and FZ research is supported by the NWO project 613.009.10. MV is supported by the DFG through the SFB Transregio 109 “Discretization in Geometry and Dynamics”.

2 Contact dynamics

In this section we review briefly the basic notions of contact geometry and dynamics, both from the Hamiltonian and from the Lagrangian perspective.

2.1 Contact geometry

Contact geometry is the “odd–dimensional cousin” of symplectic geometry. However, it has received much less attention in the scientific literature until very recently (with the exception of contact topology 

[geiges2008introduction]). One of the reasons may be that the definition of a contact manifold is somewhat involved, and therefore we decide here to present the general definition together with a more restrictive one, which is all is needed for the present work. In general, a contact manifold is a –dimensional manifold endowed with a contact structure

, that is, a maximally non–integrable distribution of hyperplanes, meaning that its integrable submanifolds have dimension at most

. Here, borrowing the terminology from [ciaglia2018contact], we restrict to the following important case.

Definition 1.

An exact contact manifold is a –dimensional manifold endowed with a contact structure which is given globally as the kernel of a –form satisfying the non–degeneracy condition .

It is easy to verify that any exact contact manifold is a contact manifold, but that the converse is not true (for instance, if the manifold is not orientable, then it can admit a contact structure, but it cannot admit a global –form satisfying the non–degeneracy condition needed in the definition of an exact contact manifold). In this work we will always assume that is an exact contact manifold.

Notice also that since , we have in fact an equivalence class of –forms that define the same . These are given by re–scaling by multiplication by a non–vanishing function. In the following definition we make this statement precise.

Definition 2.

Given an exact contact manifold and a representative such that , a diffeomorphism is called a contactomorphism if it satisfies , where is the pullback induced by and

is any non–vanishing function. Analogously, a vector field

is called a contact vector field if , where is the Lie derivative and is any function (the case when is called sometimes a strict contact vector field).

Example 1.

The standard example of a contact manifold is with Cartesian coordinates , for and the standard contact structure .

2.2 Contact Hamiltonian systems

We are now in the position to define Hamiltonian systems on (exact) contact manifolds. To do so, we will first fix a representative –form and then define the given Hamiltonian vector field associated to any function , where the additional dimension is used to account for time-dependence.

Definition 3.

Let be an exact contact manifold and fix a representative –form generating the contact structure. For any function we define the corresponding contact Hamiltonian vector field by


where is the Lie derivative, and represents the interior product.

Note that the first condition in (2) is the requirement that be a contact vector field, while the second condition is the association between the vector field and its Hamiltonian function.

We refer the interested reader to the classical textbook [arnold] for an introduction to contact geometry, and to [Bravetti2017a, bravetti2018contact] for an overview of physical applications. More specifically, in [bravetti2017contact, deleon2017cos, de2019singular, gaset2019contact, gaset2019new, valcazar2018contact] one can find a detailed study of the contact geometry of dissipative systems in classical mechanics and field theories, while in [ciaglia2018contact] an approach to study dissipative systems in quantum mechanics via contact systems is investigated. Here we limit ourselves to collect some relevant properties in the following proposition.

Proposition 1.

In the neighbourhood of any point , there always exist local coordinates , called Darboux coordinates, such that . In such coordinates takes the form


where summation over repeated indices is assumed. The corresponding contact Hamiltonian equations are given by


The following property guarantees that contact vector fields come with an associated Hamiltonian function, and will be essential for the construction of the contact Hamiltonian integrators.

Proposition 2 (Isomorphism between functions and contact vector fields).

Let be an exact contact manifold and fix a representative –form . Then induces an isomorphism of Lie algebras, between the Lie algebra of contact vector fields (infinitesimal generators of contactomorphisms) with the standard Lie bracket and the Lie algebra of functions on with the Jacobi bracket defined as


In particular, to every contact vector field we can associate uniquely a contact Hamiltonian by means of (cf. the second condition in (2)). In Darboux coordinates, the Jacobi bracket reads


Even though the Jacobi bracket can be understood as a generalization of the Poisson bracket, the equations of motion are not the same as for the Poisson bracket. As one can verify by direct calculation, there holds for any function .

The following proposition is the starting point for our analysis of equation (1) in terms of contact geometry.

Proposition 3.

Equation (1) corresponds to the flow of the contact Hamiltonian


To prove the above statement, observe that the Hamiltonian equations (4)–(6) in this case read


One can immediately check that the system (10)–(11), gives exactly the system (1), while equation (12) decouples from the rest. ∎

Using contact geometry, we have immediately obtained a “Hamiltonisation” of all the dynamical systems of the form (1). This fact should not be underestimated: a Hamiltonian structure for all such systems allows us to benefit from the theory of Hamiltonian systems (extended to the contact case) and its powerful analytical and numerical tools. For instance, one can apply weak–KAM theorems and variational methods, as done e.g. in [cannarsa2019herglotz, liu2018contact, wang2016implicit, wang2019aubry].

It is important to also compare the simplicity and generality of the formulation provided here against previous attempts in the literature. For instance, in [gkolias2017hamiltonian] an algorithm for the symplectic Hamiltonisation of systems of the type (1) has been provided. However, the construction suggested there is based on a non–trivial reparameterization that requires solving an additional differential equation in order to obtain the new time variable (which in many cases cannot be done exactly, cf. [gkolias2017hamiltonian]). We stress that in our analysis we do not encounter any such complication.

2.3 Herglotz’ variational principle

As for symplectic Hamiltonian systems, the dynamics of contact Hamiltonian systems can be characterised by a variational principle. This was originally published by Herglotz in a set of lecture notes [herglotz1930vorlesungen], which might explain why it has received relatively little attention. A modern discussion of Herglotz’ variational principle can be found for example in [georgieva2003generalized, vermeeren2019contact] (see also [gaset2019contact, lazo2017action, lazo2018action] for extensions to field theories).

Definition 4 (Herglotz’ variational principle).

Let be an –dimensional manifold with local coordinates and let . For any given smooth curve we consider the initial value problem


Then the value is a functional of the curve . We say that is a critical curve if is invariant under infinitesimal variations of that vanish at the boundary of .

If the Lagrange function does not depend on , then we find

which is the usual action functional from symplectic mechanics. Hence the classical formulation of Lagrangian mechanics is a special case of Herglotz’ variational principle.

Proposition 4.

Critical curves for the Herglotz’ variational principle are characterised by the following generalised Euler–Lagrange equations:


The generalised Euler–Lagrange equations can be derived by solving the differential equation

We find


where the boundary terms on the second line vanish because variations leave the endpoints fixed. ∎

If we restrict our attention to solutions to the generalised Euler–Lagrange equations, but allow variations of the endpoint, equation (15) reduces to

where . This means that the flow consists of contact transformations with respect to the 1–form .

One can verify that for the Hamiltonian , where is written in function of and , the equations (4)–(6) are equivalent to the system consisting of equation (13) and the generalised Euler–Lagrange equations (14).

There is a natural discretisation of Herglotz’ variational principle, which was introduced in [vermeeren2019contact].

Definition 5 (Discrete Herglotz’ variational principle).

Let be an –dimensional manifold with local coordinates and let . For any given discrete curve we consider the initial value problem


Then the value is a functional of the discrete curve . We say that is a critical curve if

Equivalently, we can require that for all . By elementary calculations, this gives us the discrete generalised Euler–Lagrange equations. As in the conventional discrete calculus of variations, they can be formulated as the equality of two formulas for the momentum [marsden2001discrete].

Proposition 5.


Then solutions to the discrete Herglotz variational principle are characterised by

Furthermore, the map induced by a critical discrete curve preserves the contact structure .

Without loss of generality it is possible to take the discrete Lagrange function depending on only one instance of : . Then the discrete generalised Euler–Lagrange equations read

3 Integrators

One of the advantages of having a contact Hamiltonian structure is that many ideas for geometric integrators for symplectic systems can be carried over with relatively small effort. A study of variational integrators in the contact case has been started recently in [vermeeren2019contact]. Here we develop higher order variational and Hamiltonian integrators for contact systems.

3.1 Lagrangian

Like in the symplectic case, we can construct higher order contact variational integrators using a Galerkin discretisation. In such a discretisation, the set of curves on the time interval of one step is replaced by a finite–dimensional space of polynomials

To paremetrise this space we introduce control points , where and . If for each of these control points a value is prescribed, then the polynomial is uniquely determined. We denote by the polynomial thus obtained.

Given a continuous Lagrangian , we would like to define by specifying an initial condition and setting

To approximate , and in particular , we use an explicit Runge–Kutta method of order with coefficients , and , i.e. we calculate

and set

The discrete Lagrangian is defined by finding a critical value of and subtracting to match the formulation of the Herglotz variational principle:

where denotes the critical value with respect to variations of .

Remark 1.

Based on numerical evidence from the symplectic version of this construction [ober2015construction], we expect the variational integrator defined by this Lagrangian to be of order , where is the degree of the polynomials and the order of the Runge–Kutta method. Hence the order of the integrator can be twice the degree of the polynomial approximation. A general proof of this fact is the topic of a future work.

Example 2.

We use second order polynomials and control points , and . If for , then is given by

and its derivative by

In particular we have

We use the classical fourth order Runge–Kutta method to calculate , approximating as a function of :


This gives us the discrete Lagrangian

from which a difference equation for is obtained by the discrete Herglotz variational principle (see Proposition 5).

3.2 Hamiltonian

Here we review the standard splitting procedure for the generation of higher order methods for separable flows and then apply it to the case of contact Hamiltonian systems in which the Hamiltonian can be split into the sum of different pieces that can be integrated exactly (see also [mclachlan2002splitting]).

3.2.1 Splitting methods

We begin with a brief review of the splitting method, following closely [yoshida1990construction]. First of all, we have the following definition and a related important property.

Definition 6.

We say that a vector field is exactly integrable if there exists a solution to the differential equation , given by , that can be explicitly written in closed form.

Proposition 6 (2nd–order integrator).

If a vector field can be split as a sum


where each of the vector fields is exactly integrable, then


is a second order integrator for the differential equation

Remark 2.

Proposition 6 holds also for non–integrable vector fields . However, the requirement of exact integrability is crucial to be able to implement the corresponding integrators.

Based on a repeated use of the second order integrator (18) with appropriatly changed step sizes, Yoshida [yoshida1990construction] developed two different algorithms to construct integrators of any even order; the difference between the two is that the first one uses exact coefficients to calculate the rescaled time steps, while the second one uses approximated coefficients. Although the first method in principle is more accurate, the latter is sometimes preferred because it involves fewer calculations.

We can summarise these two approaches in the following statements.

Proposition 7 (Integrator with exact coefficients).

If is an integrator of order , then the map


with and given by


is an integrator of order .

Proposition 8 (Integrator with approximated coefficients).

There exist and a set of real coefficients such that the map


is an integrator of order .

The proof of Proposition 8 is constructive [yoshida1990construction], and the coefficients are obtained as approximated solutions to an algebraic equation derived from the Baker–Campbell–Hausdorff formula (see also equation (26) below). Table 1 lists values of the coefficients for three –order integrators, labeled A, B and C. Note that .

Table 1: The coefficients for three –order integrators.

3.2.2 Contact integrators of order

Let us now apply the above splitting schemes to derive contact integrators, i.e. integrators that preserve the contact structure. The time evolution in this case is given by a contact Hamiltonian vector field, and thus the flow is a contact map, as explained in Section 2. As a direct consequence of the splitting method, we have the following result.

Proposition 9.

Let the contact Hamiltonian be separable into the sum of functions


such that each of the vector fields is exactly integrable. Then the integrator


is a second order contact integrator.


The fact that (24) is a second order integrator follows directly from Proposition 6. Moreover, each map is a contact map because by definition it is the flow of a contact Hamiltonian vector field. Being the composition of contact maps, is a contact transformation itself. ∎

Corollary 1.

We can construct contact integrators of any even order.


Such a construction is obtained by combining Proposition 9 with either Proposition 7 or Proposition 8. ∎

Remark 3.

The general question of finding all contact Hamiltonian systems admitting a splitting into exactly integrable pieces has been addressed in [mclachlan2002splitting].

3.2.3 Modified Hamiltonian and error estimation

One of the most powerful techniques to study the long–time behaviour of symplectic or contact integrators is the so–called backward error analysis, that is, the study of modified differential equations that are exactly traced by the discrete maps of the integrator.

Propositions 2 and 9 suggest that any contact integrator has an associated modified Hamiltonian, meaning that the numeric integration follows exactly the flow of a different contact Hamiltonian. Below we show the construction for a second order integrator for a time–dependent Hamiltonian of the type (9).

We consider , then according to Proposition 9, we have the second order integrator


where we stress that the first and last terms, , are needed only in the case of non–autonomous Hamiltonians.

The Baker–Campbell–Hausdorff (BCH) formula provides a closed expression to compute the product of exponentials of any two non–commutative operators and in the Lie algebra of a Lie group [varadarajan1974lie]. More precisely, let be the solution to , then


where “…” indicates terms involving higher commutators of and .

Applying (26) and the property that to the product of exponentials in (25), we obtain, up to fourth order in :


Using the property that and the time reversibility, we see that (27) can be reduced to




is the correction to the original Hamiltonian up to order two for the autonomous case.

Now we want to compute the modified Hamiltonian for a time–dependent system. The previous steps continue to hold, but we further need to include the time dependence applying again the BCH formula in (28). Using the property


we find that the modified Hamiltonian is given by


For Hamiltonians of the type (9), taking , and , we obtain the explicit form:

Remark 4.

Note that (31) is a truncation after the second order in of the modified Hamiltonian, which is an asymptotic series. This is important to keep in mind, especially if some of the terms in (32) contain a singularity (i.e. negative order terms) in . In that case the overall result is not a term, because in the first few steps we have , so the singularity in leads to an order reduction in . Similar order reductions will take place in the higher order terms of the modified Hamiltonian, which we have not written explicitly here. We will see an example of this in Section 4.3.

The modified equation is the formal differential equation


generated by the modified Hamiltonian, where is any smooth function of the dynamical variables. It has the property that solutions to the modified equation, truncated at a certain order in

, interpolate discrete solutions up to an error of the same order in


Proposition 10.

If the Hamiltonian does not contain any singularities in the time variable, then the integrator is of second order and the local error is


where denotes the exact flow after time and the numerical integrator.


We have that

In particular, if

are Darboux coordinates, we find as local error estimates in each of the coordinates:


Consider a numerical solution and an exact solution with the same initial data. We can estimate an upper bound for the error after steps by

where because any integrator is close to the identity map. Hence as long as the error is small, , we have

which gives us an estimate for the error after steps:

Remark 5.

In the proof above we obtain an upper bound for the numerical error in an asymptotic sense. For relatively large the term will not be negligible. In addition, after several integration steps, will likely be too large, violating the assumption that . This can be seen in a few instance of the examples below.

Example 3.

The contact Hamiltonian of a damped harmonic oscillator is given by [bravetti2017contact]


From equation (31) it follows readily that


In this case Proposition 10 clearly holds, and it implies that the errors in the kinematic quantities at each step are