Casimir preserving stochastic Lie-Poisson integrators

Casimir preserving integrators for stochastic Lie-Poisson equations with Stratonovich noise are developed extending Runge-Kutta Munthe-Kaas methods. The underlying Lie-Poisson structure is preserved along stochastic trajectories. A related stochastic differential equation on the Lie algebra is derived. The solution of this differential equation updates the evolution of the Lie-Poisson dynamics by means of the exponential map. The constructed numerical method conserves Casimir-invariants exactly, which is important for long time integration. This is illustrated numerically for the case of the stochastic heavy top and the stochastic sine-Euler equations.


Introduction
Many problems in physics and chemistry are described by Hamiltonian models.The most familiar Hamiltonian model is the canonical one, described in terms of momenta and positions.The phase space associated with this kind of model is even dimensional, and the dynamics preserves the Hamiltonian itself as well as the canonical symplectic form.A famous example of a canonical Hamiltonian system is the planetary n-body problem, see Wisdom and Holman [1].The equations for the n-body problem cannot generally be solved analytically, and numerical integration is necessary.However, to avoid the collapse or divergence of the planetary orbits, one requires numerical methods that respect the symplectic structure of the n-body problem.
It often happens that canonical Hamiltonian systems are formulated on the cotangent bundle of a Lie group G and have a G-invariant Hamiltonian.Examples of such systems include the rigid body, the heavy top, certain special discretizations of two-dimensional ideal hydrodynamics, and many problems in quantum mechanics.The G-invariance leads to a differentiable symmetry, which, by Noether's theorem, implies an associated conservation law.The symplectic structure in this setting is replaced by a Lie-Poisson structure.This Lie-Poisson structure is degenerate on certain functions known as Casimirs.In the context of Lie-Poisson (LP) equations, structure-preserving integration aims to preserve these Casimirs.
In this paper, we develop geometric time integration that explicitly incorporates this underlying mathematical structure and preserves important invariants to machine accuracy.The main new result is the extension to stochastic dynamics that also preserves the underlying geometric structure.
We focus on stochastic systems that have an LP Hamiltonian formulation.Before discussing the stochastic setting, we first review the literature on deterministic LP integration.Then, we discuss Lie group integration and how LP integration and Lie group integration are related.The review here is far from complete; the monograph Hairer et al. [2] and the review papers Iserles et al. [3], Marsden and West [4], Celledoni et al. [5] provide a much more complete picture of geometric integration.After the brief review of deterministic LP equations and their numerical integration, we discuss stochastic mechanics, focusing first on the canonical formulation.We then move on to stochastic LP equations and their numerical integration, a topic of recent interest and the main topic of the present work.
The notions of LP structure and LP equations arise from the Marsden-Weinstein reduction theorem, and a fundamental result proved in Marsden and Weinstein [6].Most LP equations cannot be solved analytically, and numerical methods are required.A seminal work is Zhong and Marsden [7], where discrete methods for LP equations were introduced by means of Hamilton-Jacobi theory.This method requires a coordinatisation of the group, which can be difficult to implement.This difficulty was resolved by Channell and Scovel [8] by reformulating the LP algorithm of Zhong and Marsden [7] in terms of Lie algebra variables.Variational integrators became an alternative to the Hamilton-Jacobi formulation when, in Holm et al. [9], the Lagrangian perspective of Marsden-Weinstein reduction was introduced.The Lagrangian reduction is known as Euler-Poincaré reduction.Discrete analogs of Euler-Poincaré and LP reduction theory were introduced in Marsden et al. [10] for systems on finite-dimensional Lie groups with a G-invariant Lagrangian.A review paper that summarizes discrete mechanics and variational integrators at that time is Marsden and West [4].
In the pioneering work of Crouch and Grossman [11], numerical methods for ordinary differential equations on manifolds were developed.Finite dimensional mechanical systems form an important subclass of differential equations on manifolds.The Crouch-Grossman method, which uses the notion of frames, leads to a class of algorithms with a maximum order of convergence of three.Beyond order three, the analysis of the algorithms becomes very complex.To combat this issue, a different approach was introduced by Munthe-Kaas [12], where numerical methods for differential equations on manifolds were introduced based on Lie groups.These methods are now known as the Runge-Kutta-Munthe-Kaas (RKMK) methods.The review paper Iserles et al. [3] summarises most of what was then known about Lie group methods.
In Engø and Faltinsen [13], the Lie group methods of Munthe-Kaas [12] were applied specifically to LP equations.Here, the community working on numerical methods for differential equations on manifolds and the community working on discrete mechanics with symmetries intersect.Within this intersection, Lie group methods that preserve the structures associated with mechanics were developed by Bou-Rabee [14], Bou-Rabee and Marsden [15].In Bogfjellmo and Marthinsen [16], symplectic integrators of arbitrarily high order where developed building on the methods of Crouch and Grossman [11] and Munthe-Kaas [12].A review on Lie group methods can be found in Celledoni et al. [5].We now move out of the deterministic setting to discuss stochastic mechanics.
Since the work of Bismut [17], stochastic Hamiltonian systems have entered as important modeling tools for the analysis of continuous and discrete mechanical systems with uncertainty.In Bismut [17], stochastic Hamiltonian systems driven by Brownian motion on symplectic manifolds were introduced.In Lázaro-Camı and Ortega [18], the work of Bismut was generalized to manifolds, and they showed that the stochastic Hamiltonian systems extremise a stochastic action defined on the space of manifold-valued semimartingales.However, Lázaro-Camı and Ortega [18] provided a counterexample to the converse statement: "an extremum of stochastic action satisfies stochastic Hamilton's equations".In Bou-Rabee and Owhadi [19], the focus was restricted to stochastic Hamiltonian systems that are driven by Wiener processes.In this context, Bou-Rabee and Owhadi [19] were able to prove almost surely that a curve satisfies stochastic canonical Hamiltonian equations if and only if it extremises a stochastic action.For the numerical integration of stochastic canonical Hamiltonian systems, Bou-Rabee and Owhadi [19] introduce stochastic variational integrators in the special case when the Hamiltonian is independent of the momentum variable.In Deng et al. [20], high-order symplectic methods for stochastic canonical Hamiltonian systems were developed.A unifying framework of stochastic discrete variational integrators is developed in Holm and Tyranowski [21], where the method of Bou-Rabee and Owhadi [19] is extended to general Hamiltonian systems.Drift preserving methods for stochastic Hamiltonian systems were introduced in Chen et al. [22], and variational integrators for stochastic diffusive Hamiltonian systems were developed in Kraus and Tyranowski [23] using stochastic Lagrange-D' Alembert variational principles.
In Holm [24], a stochastic variational principle was introduced with the purpose of deriving stochastic Euler-Poincaré equations.Stochastic Euler-Poincaré equations are equivalent to stochastic LP equations when the Legendre transform is a diffeomorphism, so the stochastic variational principle provides a systematic means to derive stochastic LP equations.The type of stochasticity introduced in Holm [24] is called "stochastic advection by Lie transport (SALT)", and its purpose is to model principally unknown effects influencing transport in fluid problems.In Arnaudon et al. [25], 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.This provides the desire for structure-preserving numerical methods for stochastic LP equations.
Recently, attention has been paid to the structure-preserving numerical integration of stochastic LP equations.In Bréhier et al. [26], splitting methods for Stratonovich stochastic LP 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 LP equations.Our approach is based on the RKMK methods of Munthe-Kaas [12] and Engø and Faltinsen [13], meaning it preserves Casimirs invariants associated with the LP structure.The method has a beneficial property compared to others with regard to the strong order of convergence for the kind of noise considered in Holm [24].
The work presented here is a first step in creating high-order invariant-preserving RKMK methods for solving stochastic LP systems.The use of the RKMK framework is emphasized here as it opens up the possibility for an extension toward high-order methods in the future, such as the 1.5 and 2.0 strong-Taylor schemes illustrated in Kloeden and Platen [27], Holm and Tyranowski [21].The important step of achieving lower-order invariant-preserving geometric time integration for stochastic LP systems is detailed here.The invariant-preserving property is built into the numerical method and illustrated explicitly for the well-known heavy top problem as well as for the fluid-mechanical sine-Euler model developed in Zeitlin [28], which connects the current work to applications in geophysics.Different numerical tests are presented to underpin the practical usefulness of the new time integration.We restrict the illustrations to low-dimensional dynamics-this is not a principal restriction of the new geometric time integration method.In fact, an extension to the high-performance simulation of fully resolved spatial models of Navier-Stokes type was achieved in Cifani et al. [29] recently.
This paper is structured as follows.In Sect.2, we introduce LP brackets to which the geometric structure of the problem is associated.In Sect.3, we introduce stochastic LP dynamics by introducing semimartingale Hamiltonians.In Sect.4, we introduce the new numerical stochastic LP integrator.This integrator is able to preserve the Casimirs and, upon removing the noise, recovers the results of Engø and Faltinsen [13].In Sect.5, we apply the trapezoidal rule combined with the Munthe-Kaas method (TMK) to the stochastic heavy top and show that the Casimirs are conserved, whereas the trapezoidal rule applied directly to the LP equations fails to exactly conserve the Casimirs.In Sect.6, we provide further illustration of the TMK method for the sine-Euler equations Zeitlin [28].The latter system is characterized by higher-degree (polynomial) invariants, the conservation of which is crucial for proper capturing of the intricate dynamics associated with this system.We will show that the developed stochastic geometric integrator is, in fact, able to preserve such structure.Conclusions are collected in Sect.7.

Lie-Poisson brackets
We start by recalling some notation and definitions of Marsden and Ratiu [30].We will denote Lie groups by G and the associated Lie algebra by g.An LP bracket is a linear Poisson bracket on a vector space.Let f ∈ C ∞ (V ), with V a vector space, and let the duality pairing •, • V * ×V : V * × V → R define the dual V * .One can determine the variational derivative δf /δu ∈ V * as the unique element obtained by evaluating the Gateaux derivative of the functional f as where ad : g × g → g is the adjoint representation of the action of the Lie algebra on itself, and ad * : g × g * → g * 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 Ad : G × g → g and Ad * : G × g * → g * .For matrix Lie groups, Ad g v = gvg -1 and Ad * g -1 μ, v = μ, Ad g v defines the coadjoint action.The set O μ 0 := {Ad * g -1 μ 0 |g ∈ G} defined by the LP bracket is called the coadjoint orbit of μ 0 ∈ g * .On coadjoint orbits, the Casimirs are constant.This follows from the fact that the Casimirs are in the kernel of ad * and ad * is the infinitesimal generator for Ad * .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 LP integrators that are able to reflect this property numerically, see, for instance, Zhong and Marsden [7], Channell and Scovel [8], McLachlan [31], Reich [32], McLachlan and Scovel [33], Engø and Faltinsen [13] for different developments of numerical LP integrators.More recently the intersection of isospectral methods and LP methods led to novel LP integrators, as shown in Bloch and Iserles [34] and Modin and Viviani [35].
Henceforth, we will only consider finite-dimensional Lie algebras.By Ado's theorem (Rossmann [36], page 51), every finite-dimensional Lie algebra is isomorphic to a matrix Lie algebra.This means that the natural pairing is the Frobenius pairing.For finite dimensional Lie algebras with dimension N , let e i denote the basis so that a Lie algebra element σ ∈ g can be expressed as σ = N i=1 σ i e i .Let ε i denote the induced dual basis by the Frobenius pairing.Then, an element μ ∈ g * can be expressed as μ = N i=1 μi ε i .The LP bracket (2) can be expressed in terms of the structure constants C k ij of the Lie algebra g as follows The relation between ( 2) and ( 4) is expressed in the following lemma.

Lemma 2.1 For any σ ∈ g and any
Lemma 2.1 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 J.In the canonical case, the matrix J is the familiar symplectic matrix.The coadjoint representation is important since it can be used to solve LP equations exactly at a formal level.In addition, Lemma 2.1 defines the linear operator J, which generalizes the symplectic matrix encountered in canonical Hamiltonian dynamics.
LP brackets arise naturally after reducing the dimension in mechanical systems with symmetry, as is described for the rigid body in Smale [37,38] and in general in Marsden and Weinstein [6], Marsden and Ratiu [30], Holm [39], Holm et al. [40].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 Marsden and Weinstein [6].Since chirality plays a role, the LP 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.For μ ∈ g * , LP 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 N , the dimension of the Lie algebra g.
The LP equations associated to a Hamiltonian : which implies that the momentum μ satisfies the equation where ad * : g × g * → g * is the dual of the adjoint representation ad = [•, •] : g × g → g, 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 emphasize 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 LP integrator to solve equation (7).This is what Engø and Faltinsen [13] did, and they proceeded to write a deterministic LP integrator that preserves the Hamiltonian.We will first introduce stochasticity into the LP equations in Sect. 3 and formulate an integrator in Sect. 4.

Stochastic Lie-Poisson dynamics
Following Protter [41], we introduce a filtered probability space given by the quadruplet ( , F, (F t ) t≥0 , P).Here, is a set, F is the σ -algebra, (F t ) t≥0 is a right-continuous filtration and P is the probability measure.
With respect to the filtered probability space, we define a family W 1 t , . . ., W M t of independent, identically distributed Brownian motions.In this section, we will assume that all considered stochastic processes are compatible with the continuous semimartingale S t = (S 0 t , S 1 t , . . ., S M t ), where M ≤ dim g.Compatibility is understood in the sense of Definition 2.3 in Street and Crisan [42], 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 the composition of functions for which this symbol is also adopted.The Stratonovich integral is preferred for the development of stochastic LP equations because Stratonovich processes satisfy the ordinary chain rule, meaning that a Hamiltonian plus semimartingale noise is enough to determine the equations of motion.The Itô integral can also be used, but this requires a connection on the underlying manifold, see Émery [43], Huang and Zambrini [44].
Noise will be introduced via a semimartingale Hamiltonian s : g * → R and a continuous semimartingale S t .It is possible to define stochastic LP equations with general semimartingales.We will discuss the conditions for the existence and uniqueness of strong solutions for stochastic LP equations with general semimartingales.For our examples in Sects.5 and 6, we restrict to stochastic Hamiltonian systems driven by Wiener processes with specific diffusion Hamiltonians, as dictated by the SALT framework.In this case, the semimartingale S t takes the explicit form and semimartingale Hamiltonian splits as We associate the Hamiltonian : g * → R with the drift dt such that the deterministic LP equations are recovered upon setting the diffusion Hamiltonians i : g * → R to zero.The noise is controlled by diffusion Hamiltonians (or noise Hamiltonians) i associated with the diffusions dW i t .As in the deterministic case, the stochastic LP equations distinguish a left-invariant version and a right-invariant version.The left-and right-invariant stochastic LP equations are given, respectively, by Let the initial datum be μ(0) = μ 0 ∈ g * .When S t is a general semimartinagle, the integral form of equation (10) is for which Theorems 6 and 7 in Chap. 5 of Protter [41] can be applied after the Itô correction to prove that strong solutions to (11) Hence, a sufficient condition for strong wellposedness of ( 11) is that the semimartingale Hamiltonian is a C 2,1 function, that is, the semimartingale Hamiltonian is twice differentiable with Lipschitz continuous derivatives.One derivative is required for computing the gradient of the Hamiltonian, a second derivative is required for the Itô correction, and the result needs to be Lipschitz continuous.For weak solutions to (11), one can consider the conditions discussed in Stroock and Varadhan [45], and for more general results on conditions for strong wellposedness of solutions, we refer to Jacod [46].In case the semimartingale is of the form the conditions for strong wellposedness of solutions to the stochastic LP equation reduce to the conditions for Stratonovich diffusions, which require the drift to be Lipschitz continuous and the diffusions to be C 1,1 because one needs to do the Itô correction.This means that the drift Hamiltonian needs to be in C 1,1 and the diffusion Hamiltonians need to be in C 2,1 , i.e., twice differentiable with Lipschitz continuous derivatives.
A canonical choice for a noise Hamiltonian i is through the coupling of noise to the momentum μ ∈ g * .This corresponds to the concept of "stochastic advection by Lie transport" (SALT) that was introduced in Holm [24], Bethencourt de Leon et al. [47] and is also adopted in this paper.For SALT, the semimartingale Hamiltonian in (9) takes the form with β i ∈ g being the Lie algebra-valued noise coefficient so that β i • μ ∈ R for each i.The Lie algebra-valued noise coefficient determines the amplitude of the noise in the direction of each basis vector.Since the diffusion Hamiltonians are smooth in this case, the stochastic LP equation has unique, strong solutions provided that the drift Hamiltonian is in C 1,1 .Upon inserting the semimartingale Hamiltonian ( 14) into the LP bracket (4), the resulting stochastic LP equation will have linear multiplicative noise: The adjoint and coadjoint representation theory for Lie algebras and their dual can be used to solve the LP equations (10) formally.Adjoint representation theory features the linear operators Ad : G × g → g and ad : g × g → g and coadjoint representation theory features the linear operators Ad * : G × g * → g * and ad * : g × g * → g * .The operator Ad * is the dual of Ad and ad * is the dual of ad with respect to the pairing •, • : g * × g → R.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 emphasize the role of chirality.Lemma 3.1 Let g(t) : R + → G, σ (t) : R + → g and μ(t) : R + → g * all depend almost surely continuously on t and be compatible with the semimartingale S t .Then, the following formulas hold

Right-invariant case
Denote by dζ = (dg)g -1 the right-invariant vector field, The proof is a direct computation using the relations between the vector field dζ and the group element g as well as the definitions of the representations Ad, Ad * , ad and ad * .From the fourth equation in (16), it follows that the solution to the left-invariant stochastic LP equation is given by and from the third equation in (17), it follows that the solution to the right-invariant stochastic LP equation is given by In both cases, the curve g(t) : R → G is a continuous curve that is the solution to a stochastic differential equation related to the LP equations.Namely, the curve g(t) = exp ∂H ∂μ (μ(t)) is the exponential of the derivative of the Hamiltonian with respect to μ.These relations are made more precise in the next section.However, since the dependence of the Hamiltonian on μ is nonlinear (if the Hamiltonian is hyperregular), it is not possible to generate the entire curve g(t) from the initial condition.The curve μ(t) : R + → g * that solves the left-or right-invariant stochastic LP equation lives on the set O μ 0 defined by The set O μ 0 is the coadjoint orbit generated by the initial condition.Since Casimirs are in the kernel of the operator ad * , 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 realized by computing a numerical approximation to the group element g that itself is an element of the group G.This has been pioneered by Engø and Faltinsen [13] for deterministic LP systems.We use that as a starting point for the extension toward stochastic problems.

Numerical integration
In Engø and Faltinsen [13], implicit numerical methods that preserve Casimirs as well as the Hamiltonian were introduced, based on the Runge-Kutta Munthe-Kaas (RKMK) integrators of Munthe-Kaas [12] and the discrete gradient method of Gonzalez [48] 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 the conservation of the Hamiltonian.We will employ similar methods for the integration of stochastic LP equations, but we first introduce the RKMK method in the deterministic setting.
The RKMK method is based on canonical coordinates of the first kind.Following the notation of Bou-Rabee and Marsden [15], one introduces a map τ : g → G such that τ is a local diffeomorphism of neighborhood of 0 ∈ g to a neighbourhood of the identity e ∈ G with τ (0) = e.The map τ is assumed to be analytic in the neighbourhood of the identity and is also required to satisfy τ (X)τ (-X) = e for all elements X ∈ g.This means that τ induces a local chart on G for which the left translation can be used to form an atlas.The exponential map exp : g → G is an important example of a canonical coordinate-inducing map.One can construct RKMK methods based on any τ : g → G by computing the derivative dτ of τ and its inverse dτ -1 .For a general τ , the RKMK method is defined below.Definition 4.1 Given RK coefficients b i , a ij ∈ R (i, j = 1, . . ., s), set c i = s j=1 a ij .An sstage RKMK approximant to the differential equation ġ(t) = g(t)f (t, g(t)) with initial datum g(0) = g 0 ∈ G is given by If a ij = 0 for i < j, then the RKMK method is explicit, and it is implicit otherwise.
In most cases, dτ -1 is expressed by a series expansion that will have to be truncated to be able to use the RKMK method.Theorem 4.7 in Bou-Rabee and Marsden [15] explains that given an rth order approximation to the exact exponential map τ and an RK method of order p with r ≥ p, then the RKMK method is of order p if the truncation of dτ -1 satisfies q ≥ p -2.By Ado's theorem, we have that every finite-dimensional Lie algebra is isomorphic to a matrix Lie algebra.This means that for computational purposes, the exponential map can be expressed as the matrix exponential This means that we can compute the differential d exp and its inverse d exp -1 explicitly, which will be used in the construction of the RKMK method for stochastic LP equations.By means of the linear representations of Lie algebras and their duals, one can establish the following important relationship.Lemma 4.1 For any σ ∈ g, Lemma 4.1 is the coadjoint version of the fundamental relation Ad exp(σ ) = exp(ad σ ).The proof is a straightforward calculation, see, e.g., Rossmann [36].Lemma 4.1 shows that -ad * is the infinitesimal generator of Ad * .It also implies that Casimirs, which are in the kernel of ad * , result in the identity operator.Hence, the value of the Casimir is constant on coadjoint orbits.We will now construct the RKMK method for stochastic LP equations.We will represent a group element g = exp σ (locally) by a Lie algebra element.As a consequence, the following representation of the solution μ(t) is obtained.

Left-invariant case
The solution to the left-invariant LP equation is given by The next step is to apply the stochastic product rule to (24).This requires a derivative of the exponential map since the vector field dζ is given by dζ = g -1 dg in this case.Taking the stochastic differential and using the third equation in ( 16) yields

Right-invariant case
The solution to the right-invariant LP equation is given by The next step is to apply the stochastic product rule to (25).This requires a derivative of the exponential map since the vector field dζ is given by dζ = dgg -1 in this case.Taking the stochastic differential and using the fourth equation in (17) yields Since μ satisfies the stochastic LP equation, we can deduce that Now, we compute dζ from its definition dζ = g -1 dg with g = exp(σ ) Since μ satisfies the stochastic LP equation, we can deduce that Now, we compute ξ from its definition ξ = dgg -1 with g = exp(σ ) Using ( 28) and (30), we obtain the stochastic differential equation for σ , which is given by Using ( 29) and (31), we obtain the stochastic differential equation for σ , which is given by The operator in equation ( 32) acting on the variational derivative of the Hamiltonian is defined through its power series expansion as where B + k are the Bernoulli numbers with the convention Similarly, in the right-invariant case in (33), we have where B - k are the Bernoulli numbers with the convention B - 1 = -1/2.Note that the operator d exp -1 σ (•) : g → g is a linear operator.The iterated commutators are nonlinear in σ and linear in the other variable.The SDEs ( 32) and ( 33) can be discretized using any numerical method that is consistent with Stratonovich SDEs, see Kloeden and Platen [27] for instance.As a result, the Casimirs will be preserved by construction.This is the essence of the deterministic RKMK method developed in Munthe-Kaas [12].It was shown in Munthe-Kaas [12] that one needs as many terms in the expansions (34) and (35) 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 conserve energy, one also has to put in some more work.Engø and Faltinsen [13] used the discrete derivative method of Gonzalez [48] together with the RKMK method to obtain energy and Casimir conservation.Energy conserving methods based on the discrete derivative method of Gonzalez [48] are second order in time, which implies that only the first term in the expansions (34) and ( 35) is required.
The energy and Casimir conserving method that Engø and Faltinsen [13] developed yields a discretization based on the trapezoidal rule.In terms of the variable μ on the dual of the Lie algebra, this leads to the natural discretization for the Stratonovich integral, see, e.g., Stratonovich [49], Protter [41].It is, in general, not possible to construct arbitrarily high-order methods for stochastic differential equations based solely on increments of Wiener processes, as was shown by Clark and Cameron [50].One can use iterated integrals to go to higher order, but this is a notoriously difficult problem, see Kloeden and Platen [27].Thus, in practice, only low-order methods (strong order of convergence is usually below 1.5) are used for the integration of SDEs, and just the first term in the expansions (34) and ( 35) is required.This reduces the operator d exp -1 σ (•) to the identity operator.The SDE for σ with initial condition σ (0) = 0 is then given by The SDE ( 36) 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 (36) depends on the Hamiltonian, which in turn depends on the momentum μ.This is important for the solution procedure that we will now construct.36) has additive noise in the SALT case.This facilitates the analysis of stochastic Lie-Poisson equations and implies that the strong order of convergence of the Euler-Maruyama method is one instead of one-half.
In the following, we employ the trapezoidal rule to discretize equation (36).The discretization of (36) viewed in terms of μ coincides with the midpoint rule, but when viewed solely in terms of σ , the discretization is the trapezoidal rule.The trapezoidal rule coincides with the discrete derivative approach of Gonzalez [48], which was designed to conserve integrals of Hamiltonian dynamics.In Engø and Faltinsen [13], 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 discretization of the Stratonovich integral and, without noise, recovers the results of Engø and Faltinsen [13].Applying the discrete derivative of Gonzalez [48] to both the drift and the diffusion terms in (36) yields the following discretization where t is the time-step size, ).Note that by definition, the second term in (37) converges to the Stratonovich integral as the time step size tends to zero.The Hamiltonian depends on μ(t) and is updated every time step as follows.The differential equation ( 36) is solved for a single time step.The result is mapped to a group element by means of the exponential map every update.
This means that σ n-1 always equals the initial condition, which is zero.Without stochasticity, (37) coincides with the method of Engø and Faltinsen [13] that conserves the deterministic Hamiltonian.To find an estimate for μ n+1 , we use a type of quasi-Newton method called the chord method to find the root of a function.This function will be equation ( 37) with σ n subtracted from both sides.We will use Lemma 4.1 to write this function as We now want to determine when f (σ n ) = 0, as this will yield the σ n that is required to go from μ n to μ n+1 .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 σ [0] n = 0 results in Here D 2 denotes the Hessian of a functional .Hence, the chord method is given by The chord method restricts the maximum stepsize, what the maximum stepsize is depends highly on the problem at hand, which will be particularly clear in the examples in Sects.5 and 6.We conclude this section by summarising the stochastic Lie-Poisson integrator based on the trapezoidal rule.In Sect.5, we will apply the stochastic Lie-Poisson integrator specified in Algorithm 1 above to the stochastic heavy top.

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 se(3) ⊂ R 4×4 associated with the special Euclidean group SE(3) ⊂ R 4×4 .The special Euclidean group is the group of rotations and translations.A representation of SE( 3) is as follows where R ∈ SO(3) ⊂ R 3×3 is a rotation matrix.The Lie algebra se(3) has the associated representation where ξ ∈ so(3) ⊂ R 3×3 is a traceless, skew-symmetric matrix.The matrix exponential applied to an element of se(3) yields an element of SE (3).Traceless, skew-symmetric matrices ξ ∈ R 3×3 have three degrees of freedom and can be represented by a vector in ξ ∈ R 3 via the hat map isomorphism : R 3 → so(3).The hat map isomorphism permits us to represent all variables associated with the heavy top in the sequel by vectors in R 3 .The semimartingale Hamiltonian for the stochastic heavy top was chosen to be Here, π ∈ R 3 is the angular momentum vector.The moment of inertia tensor is represented by the diagonal matrix I = diag(I 1 , I 2 , I 3 ) ∈ R 3×3 .The vector γ ∈ R 3 is called the gravity vector and tracks the direction of gravity relative to the orientation of the top.The vector χ ∈ R 3 connects the point around which the top rotates towards the center of gravity of the top.The vector α ∈ R 3 is the vector of noise amplitudes.The semimartingale Hamiltonian (44) 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 π(0) = π 0 ∈ R 3 and γ (0) = γ 0 ∈ R 3 .The heavy top is a completely integrable system in a number of situations.If I 1 = I 2 = I 3 , then the heavy top is known as the Euler top, which is completely integrable.If two moments of inertia are the same and the center of gravity lies on the symmetry axis, then we deal with a Lagrange top.This is also a completely integrable system, as shown in Ratiu [51].Two more integrable cases are known; the heavy top is a completely integrable system if I 1 = I 2 = 2I 3 (the Kovalevskaya top) or if I 1 = I 2 = 4I 3 (the Goryachev-Chaplygin top).We investigate two cases.In the first case, we introduce SALT-type stochasticity to the Goryachev-Chaplygin top, and in the second case, we focus on a nonintegrable situation.We will compare standard implicit midpoint (IM) integration applied directly to the stochastic LP equations (45) to the integrator introduced above (we will refer to this integrator as the Trapezoidal Munthe-Kaas Figure 1 A deterministic trajectory of the heavy top with a moment of inertia given by I = diag(4, 4, 1) generated with the TMK method.The angular momentum variable (dashed) π is at constant z value in the xy-plane.The integrable nature of this particular setting explains the regular periodic orbits method or TMK method for short).The dynamics of the gravity vector γ (t) is restricted to the sphere with radius |γ 0 | as a result of the so(3) action, whereas the dynamics of the angular momentum π (t) wanders more freely in R 3 .A deterministic trajectory of the symmetric heavy top with moment of inertia I = diag(4, 4, 1) is displayed in Fig. 1.
For stochastic simulations, we first take the Goryachev-Chaplygin top by setting the moment of inertia to be I = diag (4,4,1).A second simulation concerns a nonintegrable case with moment of inertia I = diag(4, 2, 1).Realizations of these two stochastic heavy tops are shown in Figs. 3. The sphere with radius |γ 0 | is shaded in grey.For the simulations, we have used time step size t = 0.01 and simulated until time T = 100.If one takes step size beyond a certain value t ≥ t c (whose precise value depends on which stochastic LP system one is considering, the parameters of said LP system, and the realization of the noise), the properties of the method deteriorate.For simulations of the Goryachev-Chaplygin top, this critical value was empirically found to be typically around t c ≈ 0.8 and for the nonintegrable case t c ≈ 0.6.The time step size t = 0.01 is chosen so that it is an order of magnitude smaller than the critical values in both cases.
From Figs. 4 and 5, it is evident that the IM is not able to conserve the Casimirs exactly, in contrast to the TMK method.In Fig. 6, the relative error of both Casimirs is plotted for the TMK method, which shows that TMK is able to represent the conservation of Casimirs to machine accuracy.The linear trend that is visible in the evolution of the Casimirs for the TMK method in Fig. 6 is a result of the accumulation of round-off errors since the exponential map is computed up to machine precision rather than exactly.For special cases where one has closed-form expressions for the exponential map, such as for the rigid body, one obtains exact conservation.Extrapolating this linear growth would mean that after approximately 10 10 time steps, the error in the TMK method would be comparable to the error of the IM method.Upon removing the noise, by setting α = 0, the TMK method gains conservation of the Hamiltonian as an additional property.This is shown in Fig. 7.In this section, the stochastic LP integrator based on the TMK method is illustrated for the sine-Euler equations introduced in Zeitlin [28].For completeness, we give a brief in- where {•, •} is the Poisson bracket, defined as and is the Hamiltonian with ψ the stream function, linked to the vorticity through ψ = ω.System ( 46) is a Lie-Poisson system on C ∞ (T 2 ), which admits infinitely many Casimir functions i.e., the integrated powers of vorticity are invariants of motion.In the Fourier space, equations (46) become where m = (m 1 , m 2 ) is an integer vector, and m ∧ n = m 1 n 2m 2 n 1 .Numerical simulation of (50) requires a truncation to a finite set of Fourier modes.The strict conservation of the Casimirs is then no longer maintained.In fact, to preserve the underlying geometric structure also upon truncation at arbitrary finite order, Zeitlin [28] put forward an approach based on the theory of geometric quantization of Hoppe [52].In this context, quantization refers to the process of constructing a Lie algebra of N × N complex matrices, which replaces the Poisson bracket by the matrix commutator.In particular, there exists a basis of su(N) (skew-Hermitian traceless matrices) with structure constants converging to those of the Fourier basis of C ∞ (T 2 ).In this basis, one can rewrite (46) in terms of the vorticity matrix W : where W , P ∈ su(N) with P the stream matrix.As pointed out in Zeitlin [28], traces of powers of W , are conserved by (51).This is the discrete analog of conservation of integrated powers of vorticity in the continuum.In Fourier coordinates, (51) gives the sine-Euler equations: with K = (N -1)/2 and ε = 2π/N .Following the SALT approach, we derive the stochastic Euler equations by introducing diffusion Hamiltonians To illustrate the TMK method specified in Sects. 2 and 3 for the sine-Euler equations, we simulate the dynamics for N = 3. Differently from the heavy top (Sect.5), this system has a quadratic as well as a cubic Casimir, being C 2 = Tr(W 2 ) and C 3 = Tr(W 3 ), respectively.By simulating the sine-Euler equations for N > 2, we demonstrate the conservation of polynomial invariants of high order.
As initial condition, we set the Fourier coefficients of vorticity to be randomly distributed over a uniform distribution between -0.5 and 0.5.The time step is t = 0.5, and the final simulation time is T = 10,000 to reach a well-developed statistically stationary solution.Noise is injected at modes n = [(1, 1), (1, -1)] with amplitude 10 -1 , representing a high-frequency disturbance.Each mode has its own independent Wiener process.At any time t, the state of the system is given by the four-dimensional state vector (ω 0,1 , ω 1,1 , ω 1,0 , ω 1,-1 ).Conservation of the Casimirs is shown in Fig. 8, and the Casimir behavior corresponding to direct integration of the sine-Euler equations using the trapezoidal rule is shown in Fig. 9.For the sine-Euler equations, the TMK method was empirically found to be stable even for time step sizes beyond t = 500.This is far beyond the time step size where the direct trapezoidal rule breaks down, which was empirically found to be around t = 2.For a fair comparison, we choose t = 0.5 since both methods are stable.
By construction, the TMK integrator conserves the Casimirs up to machine precision.Furthermore, analogously to the heavy top, the relative error appears to be robust and increases approximately linearly, as indicated by the dashed line in Fig. 8, over a rather long simulation time.

Conclusion
In this paper, we introduced numerical methods for stochastic Lie-Poisson (LP) equations that preserve the coadjoint orbit structure by extending the Runge-Kutta-Munthe-Figure 9 Relative error of the Casimirs for one stochastic realization of the sine-Euler equations simulated using the trapezoidal rule applied directly to the LP equations.The solid line refers to the conservation of C 2 , while the dash-dotted line refers to the conservation of C 3 Kaas (RKMK) method to stochastic differential equations (SDEs).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 to rigorously extend the approach to stochastic dynamics.The noise was chosen to be of the stochastic advection by Lie transport (SALT) type, i.e., multiplicative noise coupled linearly to the momentum variable.Moreover, given that the Stratonovich integral is defined at the midpoint of the integrand, the implicit midpoint rule is a natural choice for the integration of stochastic LP equations with Stratonovich noise.We applied the Munthe-Kaas approach to obtain a stochastic differential equation on the Lie algebra, which we solved 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.Exact generating coadjoint orbits guarantees the conservation of Casimirs.The implicit midpoint rule is particularly convenient because in the absence of the noise, the integrator preserves the deterministic Hamiltonian.
The implied SALT-induced SDE on the Lie algebra has additive noise, whereas the stochastic LP 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 LP equations directly and used the trapezoidal rule (TMK) within the class of stochastic LP integrators.A comparison of IM and TMK showed that TMK performs according to theory when the conservation of Casimirs is concerned, apart from small trends in the error due to round-off effects.This was considered for the examples of the heavy top and for low-order truncation of the sine-Euler equations.The numerical illustrations clarify that the implementation of the numerical integrators is fully in line with the theoretical preservation properties.By switching off the noise and using TMK to solve for the deterministic heavy top, we also showed that TMK conserves the Hamiltonian.
The developed stochastic LP integrator is invaluable for the long-time simulation of stochastic mechanical systems for which the conservation of the geometric structure is essential.Apart from the two test cases, i.e., the heavy top and the sine-Euler equations, the new LP integrator can be applied effectively to a range of applications, among which are robotics, celestial mechanics, biomechanics, rigid body mechanics, etc.These LP mechanical systems are perturbed stochastically with the SALT method from Holm [24], which also enables a basis for uncertainty quantification through stochastic forcing.An extension to stochastic affine LP equations, which arise in machine learning and mechanics on centrally extended Lie algebras, is left for future work.

Remark 4 . 2
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 (

Figure 2 AFigure 3 AFigure 4 A
Figure2A single realization of the stochastic heavy top generated with the TMK method with I = diag(4, 4, 1).The angular momentum (dashed) π wanders around the sphere, and the gravity vector γ exactly stays on the sphere

Figure 5 AFigure 6 A
Figure 5 A single realization of the Casimir |γ | 2 compared to the initial value.Results generated by the IM method applied directly to the Lie-Poisson equation and by the TMK method are compared

Figure 7
Figure 7 Upon setting α = 0, the noise is removed.In the absence of noise, the Hamiltonian of the heavy top is conserved by the TMK method

Figure 8
Figure 8Relative error of the Casimirs for one stochastic realization of the sine-Euler equations simulated using the TMK method.The solid line refers to the conservation of C 2 , while the dash-dotted line refers to the conservation of C 3 . The variational derivatives δf /δμ, δg/δμ are elements of the dual of the dual Lie algebra g * * ( g in finite dimensions), and •, • : g * × g → R is the nondegenerate pairing between the Lie algebra and its dual.This formulation is valid for both finite-and infinite-dimensional Lie algebras.In general, LP brackets are degenerate.This means that there exist functions c ∈ C ∞ (g * ) called Casimirs such that {f , c} ∓ = 0 for all f ∈ C ∞ (g * ).Casimirs are elements of the kernel of the LP bracket, which means that they are conserved by LP dynamics.
where f , g ∈ C ∞ (g * ), μ ∈ g * (46)ing s • dS t in(46)and applying the sine function as in (53), one arrives at the stochastic sine-Euler equations in the Fourier space