Many problems in physics and chemistry are described by Hamiltonian models. The most common Hamiltonian model is the canonical one, described in terms of momenta and positions. The phase space associated to this kind of model is even dimensional and the dynamics preservers the Hamiltonian itself as well as the canonical symplectic form. A famous example of a canonical Hamiltonian system is the planetary -body problem, see wisdom1991symplectic. The equations for the -body problem cannot be solved analytically in general and numerical integration is necessary. However, for long time stability of the planetary orbits, one requires numerical methods that respect the symplectic structure of the -body problem. Many canonical Hamiltonian systems have symmetries. By symmetry reduction one can obtain a lower dimensional equivalent system of equations.
These systems have a Lie-Poisson Hamiltonian formulation. The numerical solution of deterministic Lie-Poisson systems was discussed in engo2001numerical based on the Lie group methods of munthe1999high
. From a mathematical point of view, Lie-Poisson systems are interesting because of their differential geometric properties. Lie-Poisson brackets are degenerate on functions called Casimirs. Hence, Lie-Poisson dynamics preserves the value of the Casimirs. A numerical method for the integration of ordinary differential equations (ODEs) on manifolds was introduced inmunthe1999high with arbitrarily high order of convergence. This is now known as the Runge-Kutta-Munthe-Kaas (RKMK) method and is used in engo2001numerical to integrate Lie-Poisson equations such that the the Casimirs and the Hamiltonian are exactly preserved. This seminal work established numerical methods for deterministic systems. In this paper we extend this work to stochastic problems.
In this paper we consider structure-preserving numerical integration of stochastic Hamiltonian systems of the Lie-Poisson type. In particular, we construct a numerical integrator which is able to preserve Casimir invariants that are associated with the Lie-Poisson structure of the underlying differential equations. A natural choice for discretisation of the stochastic model is the implicit midpoint rule. In absence of noise, the midpoint formulation recovers the energy preserving method of engo2001numerical.
Since the work of bismut1982mecanique, stochastic Hamiltonian systems have entered as important modelling tools for the analysis of continuous and discrete mechanical systems with uncertainty. In bismut1982mecanique, stochastic Hamiltonian systems driven by Brownian motion on symplectic manifolds were introduced. With regards to the integration of stochastic canonical Hamiltonian systems, in bou2009stochastic stochastic variational integrators are introduced for a specific type of noise. In deng2014high high order symplectic methods for stochastic canonical Hamiltonian systems were developed. Stochastic discrete variational integrators were developed in holm2018stochastic. Drift preserving methods for stochastic Hamiltonian systems were introduced in chen2020drift and variational integrators for stochastic diffusive Hamiltonian systems were developed in kraus2021variational.
Stochastic Hamiltonian systems of the Lie-Poisson type have recently received attention. In brehier2021splitting splitting methods for Stratonovich stochastic Lie-Poisson equations are introduced. These methods are explicit and preserve Casimirs as well as the Poisson map property. We will take a different approach to the numerical solution of Stratonovich stochastic Lie-Poisson equations. Our approach is based on the RKMK methods of munthe1999high and engo2001numerical, meaning that an arbitrarily high order of convergence can be obtained.
A major theoretical contribution to stochastic mechanics is lazaro2008stochastic, which extended the results of bismut1982mecanique by describing stochastic Hamiltonian systems on Poisson manifolds driven by continuous semimartingales. Many canonical Hamiltonian models have Hamiltonians that are invariant under the action of a particular Lie group. Such systems are formulated on the cotangent bundle of a Lie group. If the Hamiltonian itself is invariant under the (left or right) action of that Lie group, then one can halve the dimension of the original problem. The resulting system is referred to as a (left or right) Lie-Poisson Hamiltonian system. The rigid body, the heavy top, ideal fluid dynamics with and without advected quantities, the Poisson-Vlasov equations in plasma physics and the Maxwell-Bloch problem in optics are all examples of Lie-Poisson systems. These examples are all treated in holm1985nonlinear.
In holm2015variational, a stochastic variational principle was introduced with the purpose of deriving stochastic Euler-Poincaré equations. This type of stochasticity is called ”stochastic advection by Lie transport (SALT)” and is introduced to model principally unknown effects influencing transport in fluid problems. In arnaudon2018noise, this stochastic variational principle is used to derive equations for finite dimensional mechanical systems. The SALT noise belongs to a class of “Hamiltonian noise”, which does not affect the Poisson bracket. This means that the Casimirs for a system that is perturbed with Hamiltonian noise are the same as for the unperturbed system. Following the SALT approach, we develop a stochastic Lie-Poisson integrator for finite dimensional mechanical systems perturbed by Hamiltonian noise.
In section 2 we introduce Lie-Poisson brackets to which the geometric structure of the problem is associated. In section 3 we introduce stochastic Lie-Poisson dynamics by introducing semimartingale Hamiltonians. In section 4 we introduce the new numerical stochastic Lie-Poisson integrator. This integrator is able to preserve the Casimirs and upon removing the noise, recovers the results of engo2001numerical. In section 5 we apply the implicit midpoint rule combined with the Munthe-Kaas method (IMMK) to the stochastic heavy top and show that the Casimirs are conserved, whereas the implicit midpoint rule applied directly to the Lie-Poisson equations fails to exactly conserve the Casimirs. Conclusions are collected in section 6.
2 Lie-Poisson brackets
We start by recalling some notation and definitions of marsden2013introduction. We will denote Lie groups by and the associated Lie algebra by
. A Lie-Poisson bracket is a linear Poisson bracket on a vector space. Let, with a vector space, and let the duality pairing define the dual . One can determine the variational derivative as the unique element defined by
with and arbitrary. Vector spaces with linear Poisson brackets are duals of Lie algebras. This follows from the fact that the dual of any Lie algebra carries a linear Poisson bracket
where , . The variational derivatives are elements of the dual of the dual Lie algebra and is the nondegenerate pairing between the Lie algebra and its dual. This formulation is valid for both finite and infinite dimensional Lie algebras. Lie-Poisson brackets are degenerate. This means that there exist functions called Casimirs such that for all .
Casimirs are elements of the kernel of the Lie-Poisson bracket, which means that they are conserved by Lie-Poisson dynamics. In the adjoint and coadjoint representation theory of the Lie algebra, the Lie bracket can be expressed as
where is the adjoint representation of the action of the Lie algebra on itself and is the coadjoint representation of the action of the Lie algebra on its dual. There are also representations of the action of the Lie group on the Lie algebra and its dual, given by and . For matrix Lie groups, and defines the coadjoint action.
The set defined by the Lie-Poisson bracket is called the coadjoint orbit of . On coadjoint orbits, the Casimirs are constant. This follows from the fact that the Casimirs are in the kernel of and is the infinitesimal generator for . The coadjoint orbits will be key in the construction of the Casimir-preserving numerical integrator. For deterministic systems, there is a large class of different Lie-Poisson integrators that are able to reflect this property numerically, see for instance zhong1988lie, channell1991integrators, mclachlan1993explicit, reich1994momentum, mclachlan1995equivariant, engo2001numerical for different developments of numerical Lie-Poisson integrators. More recently the intersection of isospectral methods and Lie-Poisson methods led to novel Lie-Poisson integrators as shown in bloch2006isospectral and modin2020lie.
In the sequel, we will only consider finite dimensional Lie algebras. For finite dimensional Lie algebras with dimension , let denote the basis so that a Lie algebra element can be expressed as . Let denote the induced dual basis. Then an element can be expressed as . The Lie-Poisson bracket (2) can be expressed in terms of the structure constants of the Lie algebra as follows
For any and any
establishes a useful relation between the coadjoint representation of the Lie algebra on its dual, the structure constants of the Lie algebra and the skew-symmetric matrix. In the canonical case, the matrix is the familiar symplectic matrix. The coadjoint representation is important since it can be used to solve Lie-Poisson equations exactly at a formal level. In addition lemma 2.1 defines the linear operator , which generalises the symplectic matrix encountered in canonical Hamiltonian dynamics.
Lie-Poisson brackets arise naturally after reducing the dimension in mechanical systems with symmetry, as is described for the rigid body in smale1970topologya, smale1970topologyb and in general in marsden1974reduction, marsden2013introduction, holm2008geometric2, holm2009geometric. Such mechanical systems can be formulated on Lie groups. When the Hamiltonian is invariant under the left or right action of that Lie group, one can perform symmetry reduction and obtain left or right Lie-Poisson equations as is established in marsden1974reduction. Since chirality plays a role, the Lie-Poisson brackets (2) and (4) feature , with the minus sign corresponding to the left invariant situation and the plus sign corresponding to the right invariant situation. In holm1998euler, the reduction theory is established for Lagrangian systems. This theory is called Euler-Poincaré reduction and is equivalent to Lie-Poisson reduction if the Legendre transform is a diffeomorphism, as shown in holm1998euler. For , Lie-Poisson equations can be formulated in several equivalent forms. For the sake of compact notation, we use Einstein’s summation convention, i.e., the summation will be understood over lower and upper pairs of indices, where the summation runs to , the dimension of the Lie algebra .
The Lie-Poisson equations associated to a Hamiltonian are given for a general by
which implies that the momentum satisfies the equation
where is the dual of the adjoint representation and we have used lemma 2.1.
The left- and right-invariant cases indicate that chirality introduces a sign difference at this stage. However, chirality expresses itself more subtly in several subsequent derivations. To emphasise the (small) differences, we have occasionally split the text into two columns, one column corresponding to the left-invariant case and the other column corresponding to the right-invariant case. At this stage, one can formulate a deterministic Lie-Poisson integrator to solve equation (7). This is what engo2001numerical did and they proceeded to write a deterministic Lie-Poisson integrator that preserves the Hamiltonian. We will first introduce stochasticity into the Lie-Poisson equations in section 3 and formulate an integrator in section 4.
3 Stochastic Lie-Poisson dynamics
, we introduce a filtered probability space given by the quadruplet. Here is a set, is the -algebra, is a right-continuous filtration and is the probability measure.
With respect to the filtered probability space, one can define a family of independent, identically distributed Brownian motions. In this section, we will assume that all considered stochastic processes are compatible with the continuous semimartingale , where . For computations, we will use . Compatibility is understood in the sense of definition 2.3 in street2021semi, which says that all stochastic processes have Radon-Nikodym derivatives with respect to each other. The symbol
will indicate that the stochastic integral is to be understood in the Stratonovich sense, instead of composition of functions for which this symbol is also adopted. The Stratonovich integral is essential for the development of stochastic Lie-Poisson equations, because Stratonovich processes satisfy the ordinary chain rule. This is required for the preservation of the Casimirs. Noise will be introduced via a semimartingale Hamiltonian, that is constructed in the following general way
We associate the Hamiltonian with the drift such that the deterministic Lie-Poisson equations are recovered upon setting the diffusion Hamiltonians to zero. The amplitude of the noise is controlled by diffusion Hamiltonians (or noise Hamiltonians) which are associated with the diffusions . The noise Hamiltonians do not affect the Lie-Poisson structure since in this paper these are integrated against a Stratonovich integral. A canonical choice for a noise Hamiltonian is through coupling of noise to the momentum . This corresponds to the concept of “stochastic advection by Lie transport” (SALT) that was introduced in holm2015variational, de2020implications and is also adopted in this paper. For SALT, the semimartingale Hamiltonian in (8) takes the form
with being the noise amplitude. Upon inserting the Hamiltonian (9) into the Lie-Poisson bracket, the resulting SDEs will have linear multiplicative noise. This means that the linear growth condition is satisfied. This is important for wellposedness of the underlying SDEs.
As in the deterministic case, the stochastic Lie-Poisson equations distinguish a left-invariant version and a right-invariant version. We define the semimartingale Hamiltonian as
Then the left- and right-invariant stochastic Lie-Poisson equations are given by
The adjoint and coadjoint representation theory for Lie algebras and their dual can be used to solve the Lie-Poisson equations (11) formally. Adjoint representation theory features the linear operators and and coadjoint representation theory features the linear operators and . The operator is the dual of and is the dual of with respect to the pairing . These operators play a fundamental role in the solution of equations on Lie algebras. The following lemma shows the differential relations between the operators, where we have split the text into columns to emphasise the role of chirality.
Let , and all depend almost surely continuously on and be compatible with the semimartingale . Then the following formulas hold
The proof is a direct computation using the relations between the vector field and the group element as well as the definitions of the representations and . From the fourth equation in (LABEL:seq:leftcoadorbitrelations) it follows that the solution to the left-invariant stochastic Lie-Poisson equation is given by
and from the third equation in (LABEL:seq:rightcoadorbitrelations) it follows that the solution to the right-invariant stochastic Lie-Poisson equation is given by
In both cases, the curve is a continuous curve that is the solution to a stochastic differential equation related to the Lie-Poisson equations. The curve that solves the left- or right-invariant stochastic Lie-Poisson equation lives on the set defined by
The set is the coadjoint orbit generated by the initial condition. Since Casimirs are in the kernel of the operator , this implies that Casimirs are constant on coadjoint orbits. To construct a numerical method that preserves the Casimirs exactly, we need to numerically compute the coadjoint orbits exactly. This can be realised by computing a numerical approximation to the group element that itself is an element of the group . This has been pioneered by engo2001numerical for deterministic Lie-Poisson systems. We use that as starting point for the extension toward stochastic problems.
4 Numerical integration
In engo2001numerical, implicit numerical methods that preserve Casimirs as well as the Hamiltonian were introduced, based on the Runge-Kutta Munthe-Kaas integrators of munthe1999high and the discrete gradient method of gonzalez1996time for integration of ordinary differential equations on manifolds. The discrete gradient method is designed for the conservation of first integrals, which can be used to obtain conservation of the Hamiltonian. We will employ similar methods for the integration of stochastic Lie-Poisson equations, which rely on the following lemma.
For any ,
Lemma 4.1 is the coadjoint version of the fundamental relation . For details of the proof rossmann2006lie can be consulted. Lemma 4.1 shows that is the infinitesimal generator of . It also implies that Casimirs, which are in the kernel of result in the identity operator. Hence the value of the Casimir is constant on coadjoint orbits. In the integrator, lemma 4.1 will be used to represent a group element (locally) by a Lie algebra element. As a consequence, the following representation of the solution is obtained.
The operator in equation (22) acting on the variational derivative of the Hamiltonian is defined through its power series expansion as
where are the Bernoulli numbers with the convention and is defined recursively as . Similarly, in the right-invariant case in (27), we have
where are the Bernoulli numbers with the convention .
Note that the operator is a linear operator. The iterated commutators are nonlinear in and linear in the other variable. The SDEs (22) and (27) can be discretised using any numerical method that is consistent with Stratonovich SDEs, see kloeden1992stochastic for instance. As a result, the Casimirs will be preserved by construction.
This is the essence of the deterministic RKMK method developed in munthe1999high. It was shown in munthe1999high that one needs as many terms in the expansions (28) and (29) as the order of convergence of the method minus one to preserve the order of convergence of the overall method. Casimirs are guaranteed to be conserved numerically, but to also conserve energy one has to put in some more work. engo2001numerical used the discrete derivative method of gonzalez1996time together with the RKMK method to obtain energy and Casimir conservation. Energy conserving methods based on the discrete derivative method of gonzalez1996time are second order in time, which implies that only the first term in the expansions (28) and (29) is required.
The energy and Casimir conserving method that engo2001numerical developed yields an implicit midpoint discretisation. This is also a natural discretisation for the Stratonovich integral, see e.g., stratonovich1966new, protter2005stochastic. It is possible to construct arbitrarily high order methods for stochastic differential equations, but such methods often require integrating iterated integrals with respect to Brownian motion, which is notoriously challenging. Thus, in practice only the first term in the expansions (28) and (29) is required. This reduces the operator to the identity operator. The SDE for with initial condition is then given by
The SDE (30) can be solved with any appropriate method such that the diffusion term converges towards the Stratonovich integral. Since the Stratonovich integral is the limit of the midpoint approximation to this integral, the Heun method (explicit midpoint) and the implicit midpoint method are natural choices. Note that the SDE (30) depends on the Hamiltonian, which in turn depends on the momentum . This is important for the solution procedure that we will now construct.
Let us consider the case where the noise Hamiltonians are chosen according to SALT. This means that the noise Hamiltonians are linear in the momentum, which implies that the stochastic Lie-Poisson equations have linear multiplicative noise of the Stratonovich type. Equation (30) has additive noise in the SALT case. This may facilitate the analysis of stochastic Lie-Poisson equations.
In the following we employ the implicit midpoint method. The implicit midpoint method coincides with the discrete derivative approach of gonzalez1996time, which was designed to conserve integrals of Hamiltonian dynamics. In engo2001numerical, the discrete derivative method was used to design a class of energy conserving Lie-Poisson integrators. Here we will use the discrete derivative approach because it corresponds to the midpoint discretisation of the Stratonovich integral and in absence of noise recovers the results of engo2001numerical. Applying the discrete derivative of gonzalez1996time to both the drift and the diffusion terms in (30) yields the following discretisation
where is the time-step size, and . Note that by definition the second term in (31) converges to the Stratonovich integral as the time step size tends to zero. The Hamiltonian is updated every time step, which means that the differential equation (30) is solved for only one time step every update. This means that is always equal to the initial condition, which is zero. Without stochasticity, (31) coincides with the method of engo2001numerical
that conserves the deterministic Hamiltonian. To find an estimate for, we use a type of quasi-Newton method called the chord method to find the root of a function. This function will be equation (31) with subtracted from both sides. We will use lemma 4.1 to write this function as
We now want to determine when , as this will yield the that is required to go from to . In the Newton-Raphson method, the Jacobian has to be updated every time step and then inverted
This can be very expensive. Instead we use the chord method, which freezes the Jacobian on the initial condition. Computing the Jacobian and evaluating on results in
Here denotes the Hessian of a functional . Hence, the chord method is given by
We conclude this section by summarising the stochastic Lie-Poisson integrator based on the implicit midpoint method.
In the situation where noise is absent, the above algorithm also conserves the Hamiltonian. Lie-Poisson systems can also be integrated with variational integrators to preserve structure, as shown in bogfjellmo2016high. Stochastic generalisations of variational integrators are left for future research.
5 Heavy top
In this section we will apply the stochastic Lie-Poisson integrator to the stochastic heavy top. The stochastic heavy top is a Lie-Poisson system on the Lie algebra associated to the special Euclidean group . The special Euclidean group is the group of rotations and translations. A representation of is as follows
where is a rotation matrix. The Lie algebra has the associated representation
where is a traceless, skew-symmetric matrix. The matrix exponential applied to an element of yields an element of . Traceless, skew-symmetric matrices
have three degrees of freedom and can be represented by a vector invia the hat map isomorphism . The hat map isomorphism permits us to represent all variables associated with the heavy top in the sequal by vectors in . The semimartingale Hamiltonian for the stochastic heavy top is given by
Here is the angular momentum vector. The diagonal matrixis called the gravity vector and tracks the direction of gravity relative to the orientation of the top. The vector connects the point around which the top rotates towards the centre of gravity of the top. The vector is the vector of noise amplitudes. The semimartingale Hamiltonian (38) involves the usual Hamiltonian for the heavy top and a single noise Hamiltonian which couples noise to the momentum. The equations of motions of the stochastic heavy top are the following stochastic Lie-Poisson equations
with initial conditions and . We will compare standard implicit midpoint (IM) integration applied directly to the stochastic Lie-Poisson equations (39) to the integrator introduced above (we will refer to this integrator as implicit midpoint Munthe-Kaas method or IMMK method for short). The dynamics of the gravity vector is restricted to the sphere with radius as a result of the action, whereas the dynamics of the angular momentum wanders more freely in . This can be seen in figure 1, which depicts a trajectory of the heavy top, generated with the IMMK method. The sphere with radius is shaded in grey. For the simulations, we have used time step size and simulated until time . The moment of inertia tensor is , the initial conditions are and . The centre of gravity lies above the point around which the top rotates, so . The noise amplitude is .
In figures 2 and 3 it is evident that the IM is not able to conserve the Casimirs exactly, in contrast to the IMMK method. In figure 4 the relative error of both Casimirs is plotted for the IMMK method, which shows that IMMK is able to represent conservation of Casimirs to machine accuracy. The linear trend that is visible in the evolution of the Casimirs for the IMMK method in figure 4 is a result of accumulation of round-off errors. Extrapolating this linear growth would mean that after approximately time steps, the error in the IMMK method would be comparable to the error of the IM method. Upon removing the noise, by setting , the IMMK method gains conservation of the Hamiltonian as an additional property. This is shown in figure 5.
In this paper we introduced structure preserving numerical methods for stochastic Lie-Poisson equations by extending the RKMK method to stochastic differential equations. Specifically, this method is able to preserve Casimirs exactly. The stochasticity is defined with respect to the Stratonovich integral, since the ordinary chain rule is required. The noise was chosen to be of the SALT type, i.e., multiplicative noise coupled linearly to the momentum variable. Moreover, given that the Stratonovich integral uses the midpoint of the integrand, the implicit midpoint rule is a natural choice for the integration of stochastic Lie-Poisson equations with Stratonovich noise. We applied the Munthe-Kaas approach to obtain a stochastic differential equation on the Lie algebra which we solve using the implicit midpoint rule. The implicit midpoint method is in the class of RKMK methods. The solution of the SDE on the Lie algebra is used to generate a group element that together with the coadjoint representation of the group on the dual of the Lie algebra generates the coadjoint orbit associated with the initial condition. Exactly generating coadjoint orbits guarantees the conservation of Casimirs. The implicit midpoint rule is particularly convenient, because in absence of the noise, the integrator preserves the deterministic Hamiltonian.
This SDE on the Lie algebra has additive noise, whereas the stochastic Lie-Poisson equation on the dual of the Lie algebra has multiplicative noise. To illustrate the theoretical results with numerical experiments, we applied the implicit midpoint (IM) rule to the Lie-Poisson equations directly and used the implicit midpoint rule (IMMK) within the class of stochastic Lie-Poisson integrators. A comparison of IM and IMMK showed that IMMK performs according to theory when conservation of Casimirs is concerned, apart from small trends in the error due to round-off effects. By switching off the noise and using IMMK to solve for the deterministic heavy top, we also showed that IMMK conserves the Hamiltonian.
The developed stochastic Lie-Poisson integrator is invaluable for the long time simulation of stochastic mechanical systems for which the conservation of the geometric structure is essential. This is the case for Lie-Poisson mechanical systems that are perturbed with the SALT method from holm2015variational.
We are grateful to Raffaele D’Ambrosio, Arnout Franken, Darryl Holm and Ruiao Hu for their valuable input in the preparation of this paper. EL and SE wish to thank the Gran Sasso Science Institute for its warm hospitality during a visit in 2021. This work was performed in the project SPReStO (structure-preserving regularization and stochastic forcing for nonlinear hyperbolic PDEs), supported by a NWO TOP 1 grant.