Metriplectic Euler-Poincaré equations: smooth and discrete dynamics

1.   Introduction

In many examples of dynamics, especially in thermodynamics, it is necessary to combine the dynamical structure of Hamiltonian systems and metric systems to produce what are called metriplectic systems, as originally discussed in the work of Morrison, see ^[1,2] (see also ^[3,4]). The dynamics is determined using a Poisson bracket for the Hamiltonian part, combined with a symmetric bracket which allows us to introduce dissipative effects (see ^[5] for a more inclusive framework using a 4-bracket to describe dissipative dynamics preserving energy and producing entropy).

After introducing the notion of metriplectic system, in this paper we study metriplectic systems derived from a perturbation of the Euler-Poincaré equations or a Lie-Poisson system by adding a special dissipation term ^[6,7]. Recall that the Euler-Poincaré equations are obtained by reduction from invariant Lagrangian systems on the tangent bundle $TG$ of a Lie group $G$. The dissipation term that we add to the equations makes the equations of motion verify two interesting properties: preservation of energy $H$ and also the existence of a Casimir function $S$ of the Lie-Poisson bracket verifying the property $\dot{S}\geq 0$. Both correspond exactly with the two laws of thermodynamics: preservation of the total energy and irreversible entropy creation.

To numerically approximate the solutions of a metriplectic system while preserving the energy and the entropy behavior it is natural to use a class of geometric integrators called discrete gradient methods. These methods are adequate when we want to preserve exactly the energy of the system. In this sense, they are quite useful for ODEs of the form $\dot{x} = \Pi(x)\nabla H(x)$ with $x\in {\mathbb R}^n$ and $\Pi(x)$ a skew-symmetric matrix (not necessarily associated to a Poisson bracket). Using a discrete gradient $\bar{\nabla} H(x, x')$ as an adequate approximation of the differential of the Hamiltonian function (see Section 4 for more details), it is possible to define a class of integrators $x'-x = \bar{\Pi}(x, x')\bar{\nabla}H(x, x')$ preserving the energy $H$ exactly, i.e. $H(x) = H(x')$. Here $\bar{\Pi}(x, x')$ is a skew-symmetric matrix approximating $\Pi(x)$. In Section 4, based on discrete gradient methods, we derive geometric integrators for metriplectic systems and in particular, the geometric derivation of the discrete dissipative term.

2.   Metriplectic systems

The theory of metriplectic systems tries to combine together the dynamics generated by Poisson brackets with additional dissipative effects. We will first review the different geometric elements that define a metriplectic system.

2.1.   Poisson structures

Consider a differentiable manifold $P$ equipped with a Poisson structure ^[8,9] given by a bilinear map

called the Poisson bracket, satisfying the following properties:

(ⅰ) Skew-symmetry, $\{g, f\} = -\{f, g\}$;

(ⅱ) Leibniz rule, $\{fg, h\} = f\{g, h\}+g\{f, h\}$;

(ⅲ) Jacobi identity, $\{\{f, g\}, h\}+\{\{h, f\}, g\}+\{\{g, h\}, f\} = 0$;

for all $f, g, h\in C^{\infty}(P)$.

Given a Poisson manifold with bracket $\{\, , \, \}$ and a function $f\in C^{\infty}(P)$ we may associate with $f$ a unique vector field $X_f\in {\mathfrak X}(P)$, the Hamiltonian vector field defined by $X_f(g) = \{g, f\}$.

Moreover, on a Poisson manifold, there exists a unique bi-vector field $\Pi$, a Poisson bivector (that is, a twice contravariant skew symmetric differentiable tensor) such that

The bivector field $\Pi$ is called the Poisson tensor and the Poisson structure is usually denoted by $(P, \{\, , \}$) or $(P, \Pi)$. The Jacobi identity in terms of the bi-vector $\Pi$ is written as $[\Pi, \Pi\, ] = 0,$ where here $[\; , \; ]$ denotes the Schouten–Nijenhuis bracket (see ^[8]).

Take coordinates $(x^i)$, $1\leq i\leq \dim P = m$, and let $\Pi^{ij}$ be the components of the Poisson bivector, that is,

Then if $f, g\in C^{\infty}(P)$ we have

and the Hamiltonian vector field is written in local coordinates as

Observe that the $m\times m$ matrix $(\Pi^{ij})$ verifies the following properties:

(ⅰ) Skew-symmetry, $\Pi^{ij} = -\Pi^{ji}$

(ⅱ) Jacobi identity,

Define ${\sharp}_{\Pi} : T^{*}P \to TP$ by

where $\alpha \in T^*P$, and $\langle \beta, \iota_\alpha \Pi\rangle = \Pi(\alpha, \beta)$ for all $\beta \in T^*P$. The rank of $\Pi$ at $p\in P$ is exactly the $\hbox{rank of }({\sharp}_\Pi)_p: T^*_pP\rightarrow T_pP$. Because of the skew-symmetry of $\Pi$, we know that the rank of $\Pi$ at a point $p\in P$ is an even integer.

Given a function $H \in C^\infty(P)$, a Hamiltonian function, we have the corresponding Hamiltonian vector field:

Therefore, on a Poisson manifold, a function $H$ determines the following dynamical system:

We say that a function $f\in C^{\infty}(P)$ is a first integral of the Hamiltonian vector field $X_H$ if for any solution $x(t)$ of Equation (2.1) we have

In other words, if $X_H(f) = 0$ or, equivalently, $\{f, H\} = 0$. In particular, the Hamiltonian function is a conserved quantity since $\{H, H\} = 0$ by the skew-symmetry of the bracket. For any Poisson manifold $(P, \Pi)$ a function $C\in C^{\infty}(P)$ is called a Casimir function of $\Pi$ if $X_C = 0$, that is, if $\{C, g\} = 0$ for all $g\in C^{\infty}(P)$.

2.2.   Positive semi-definite inner products

Assume that for each point $x\in P$ we have a positive semi-definite inner product for co-vectors

from which we can define ${\sharp}_{{\mathcal K}} : T^{*}P \to TP$ by

and the gradient vector field

for any function $S: P\rightarrow {\mathbb R}$.

${\mathcal K}$ defines a symmetric bracket given by

Take coordinates $(x^i)$, $1\leq i\leq \dim P = m$, and let $K^{ij}$ be the components of the inner product given by

Then if $f, g\in C^{\infty}(P)$, the symmetric bracket is expressed as

Observe that the $m\times m$ matrix $(K^{ij})$ verifies $K^{ij} = K^{ji}$ and all the eigenvalues are positive or zero.

2.2.1.   A construction of the positive semi-definite inner product with ${\sharp}^{\mathcal K}(dH)=0$ given a Riemannian metric

In this section, we introduce a constructive way to derive the symmetric bracket ${\mathcal K}$ which is interesting in applications. Assume that $P$ is equipped with a Riemannian metric ${\mathcal G}$ inducing a positive definite inner product ${\mathcal G}^*$ on $T^*P$ (the co-metric),

defined by ${{\mathcal G}}^*(df, dg) = {\mathcal G}(\hbox{grad}^{{\mathcal G}^*}f, \hbox{grad}^{{\mathcal G}^*} g) = \hbox{grad}^{{\mathcal G}^*}f (g)$.

Since we are interested in defining a semi-definite inner product ${\mathcal K}$ such that ${\mathcal K}(dH, \cdot) = 0$ then we define

In coordinates, we have

where $C_H = g^{ij}\frac{\partial H}{\partial x^{i}}\frac{\partial H}{\partial x^{j}}$, $(g_{ij})$ are the components of the Riemannian metric in a given coordinate system and $(g^{ij})$ denotes its inverse matrix.

By construction, ${\mathcal K}$ is positive semi-definite and ${\mathcal K}(dH, \cdot) = 0$.

Remark 2.1. Additionally we can add new functions $L_a: P\rightarrow {\mathbb R}$, $1\leq a\leq N$, to this construction in such a way that ${\mathcal K}(dL_a, \cdot) = 0$, considering

where $C_{ab} = {\mathcal G}^*(dL_a, dL_b)$, $1\leq a\leq N$ and $L_1 = H$.

2.3.   Metriplectic systems

A metriplectic system consists of a smooth manifold $P$, two smooth vector bundle maps $\sharp_{\Pi}, \sharp_{{\mathcal K}}: T^*P\rightarrow TP$ verifying that $\pi_P = \tau_P\circ \sharp_{\Pi}$ and $\pi_P = \tau_P\circ \sharp_{{\mathcal K}}$ (where $\tau_P: TP\rightarrow P$ and $\pi_P: T^*P\rightarrow P$ are the canonical projections), and two functions $H, S\in C^{\infty}(P)$ called the Hamiltonian (or total energy) and the entropy of the system, such that for all $f, g\in C^{\infty}(P)$:

● $\{f, g\} = \langle df, \sharp_{\Pi}(dg)\rangle$ is a Poisson bracket ($\Pi$ denotes the Poisson bi-vector).

● $(f, g) = \langle df, \sharp_{{\mathcal K}}(dg)\rangle$ is a positive semi-definite symmetric bracket, i.e., $(\cdot, \cdot)$ is bilinear and symmetric.

● ${\sharp}_{\mathcal K}(dH) = 0$ or, equivalently, $(H, f) = 0$, $\forall f\in C^{\infty}(P)$.

● ${\sharp}_{\Pi}(dS) = 0$ or, equivalently, $\{S, f\} = 0$, $\forall f\in C^{\infty}(P)$, that is, $S$ is a Casimir function for the Poisson bracket.

Consider the function $E = H+S: P\rightarrow {\mathbb R}$. Then, the dynamics of the metriplectic system is determined by

where $X_H = \sharp_{\Pi}(dH)$ and $\hbox{grad}^{\mathcal K} S = \sharp_{{\mathcal K}}(dS)$. From the equations of motion, it is simple to deduce the following:

● First law: conservation of energy, $\frac{dH}{dt} = \{H, H\}+(H, S) = 0$

● Second law: entropy production, $\frac{dS}{dt} = (S, S)\geq 0$.

Thus, metriplectic dynamics embodies both the first and second laws of thermodynamics.

In coordinates, the dynamics of the metriplectic system is written as

or, in matrix form, as

2.4.   Symmetry preservation

Let $\Phi: G\times P\rightarrow P$ be a smooth (left) action of a Lie group $G$ on $P$, given by $\Phi(g, x) = \Phi_g(x) = g\cdot x$ with $g\in G$ and $x\in P$. Denote by ${\mathfrak g}$ the corresponding Lie algebra. The action satisfies the following properties:

● $\Phi(e, x) = x$ where $e$ is the neutral element of $G$;

● For every $g_1, g_2\in G$ and for every $x\in P$

The infinitesimal generator of the action corresponding to a Lie algebra element $\xi \in {\mathfrak g}$ is the vector field $\xi_P$ on $P$ given by

Let $P$ be a Poisson manifold with Poisson bracket $\{\; , \; \}$ and assume that the action $\Phi$ is a Poisson action, that is,

A momentum map for the action $\Phi$ is a smooth map $J: P\rightarrow {\mathfrak g}^*$ such that for each $\xi\in {\mathfrak g}$, the associated map $J_{\xi}: P\rightarrow {\mathbb R}$ defined by $J_{\xi}(x) = \langle J(x), \xi\rangle$ satisfies that $X_{J_{\xi}} = \xi_P$ for all $\xi\in {\mathfrak g}$ where $X_{J_{\xi}}(f) = \{f, J_{\xi}\}$. As a consequence, for any function $f\in C^{\infty}(P)$

If the Lie algebra ${\mathfrak g}$ acts on the Poisson manifold $P$ and admits a momentum map $J: P\rightarrow {\mathfrak g}^*$, and if $H\circ \Phi_g = H$ (which is equivalent to $\xi_P(H) = 0$ for all $\xi\in {\mathfrak g}$ assuming that $G$ is connected), then $J_{\xi}$ is a constant of the motion of $X_H$.

Additionally, for the metriplectic system we will assume that

or, equivalently, $\sharp_{\mathcal K}(J_\xi) = 0$. Then, for the metriplectic system we have

and, therefore, $J_\xi: P\rightarrow {\mathbb R}$ is a constant of motion of the metriplectic system. See ^[10] and references therein for examples of thermodynamical systems with symmetry preservation as typically seen in the case of coupled thermomechanical problems.

As a particular case of the previous construction, we will consider in the next section the case when $P = T^*G$, where $G$ is a Lie group, and we consider as a left action $\Psi_g = T^*\mathcal{L}_{g}: T^*G\rightarrow T^*G$, where $\mathcal{L}_{g}:G\rightarrow G$ is the left action. Under the symmetry conditions, the system reduces to a metriplectic system on ${\mathfrak g}^*$, the dual of the Lie algebra of $G$.

3.   Forced Euler-Poincaré equations and metriplectic dynamics

In this section we will derive a force for the Euler-Poincaré equations in such a way that the resulting system is metriplectic. Consider a Lagrangian system $l:{\mathfrak g}\rightarrow {\mathbb R}$, where ${\mathfrak g}$ is a Lie algebra, and its corresponding Euler-Poincaré equations ^[11,12]:

where $\xi\in {\mathfrak g}$ and $\langle ad_{\xi}^*\alpha, \xi'\rangle = \langle \alpha, [\xi, \xi']\rangle$ for all $\xi'\in {\mathfrak g}$ and $\alpha\in {\mathfrak g}^*$. From this equation it is clear that the energy $E_l = \langle\frac{\delta l}{ \delta \xi}, \xi\rangle-l$ of the system is preserved, that is,

However, there are other variations of this system that are subjected to external forces that also preserve energy. This class of systems is interesting in thermodynamics when we work with a closed system, as we have seen in the Subsection 2.3 (see also ^[1,6]). For instance, if we add an external force $F: {\mathfrak g}\rightarrow {\mathfrak g}^*$ of the form

where ${\mathcal F}: {\mathfrak g}\rightarrow {\mathfrak g}^*$ is an arbitrary map, then the forced Euler-Lagrange equations are

Assume that ${\mathfrak g}$ is finite dimensional and $\{e_a\}$, $1\leq a\leq n = \dim {\mathfrak g}$ is a basis of the Lie algebra with structure constants $C_{ab}^d$, that is,

and denote by $(\xi^a(t))$ the coordinates of a curve $\xi(t)\in {\mathfrak g}$. Then, the equations (3.2) are

where ${\mathcal F}(\xi) = {\mathcal F}_{d}(\xi) e^d$ and $\{e^a\}$, $1\leq a\leq n$, is the dual basis of $\{e_a\}$.

Example 3.1. In the case of $G = SO(3)$ if we identify its Lie algebra ${\mathfrak g}$ with ${\mathbb R}^3$ with the usual vector cross product then we have

as a generalization of the equations of the rigid body also preserving the total energy of the system. In particular if

then (3.2) are

where ${\mathcal F} = ({\mathcal F}_1, {\mathcal F}_2, {\mathcal F}_3)$ and $\mathbf \Omega = (\Omega_1, \Omega_2, \Omega_3)$.

Using the Legendre transformation (that we assume in the sequel that is a local diffeomorphism) we can write the forced Euler-Lagrange equations as

where $H$ is defined by $H(\mu) = \langle \mu, \xi(\mu)\rangle -L(\xi(\mu))$ and $\mu = \frac{\delta l}{\delta \xi}(\xi)$. This is a particular case of the forced Lie-Poisson equations ^[6].

Now, if $C: {\mathfrak g}^*\rightarrow {\mathbb R}$ is a Casimir function for the Lie-Poisson bracket of ${\mathfrak g}^*$, then along the evolution of the system (3.4), we have

Example 3.2. (See ^[1] and ^[13] for applications to control theory and stabilization of a rigid body). For the rigid body equations with Hamiltonian and Casimir functions given by

in induced coordinates $(\Pi_1 = I_1 \Omega_1, \Pi_2 = I_2\Omega_2, \Pi_3 = I_3\Omega_3)$ on $({\mathbb R}^3)^*\equiv {\mathbb R}^3$, Equation (3.5) is

For instance, if we take ${\mathcal F}: {\mathfrak g}\equiv {\mathbb R}^3\rightarrow {\mathfrak g}^*\equiv {\mathbb R}^3$ as

then we get

As in the case of metriplectic systems, we have a system verifying the first and second laws of thermodynamics:

These are the equations of the relaxed rigid body ^[1].

Motivated by this example and ^[1,13], we want to study the possible families of functions ${\mathcal F}: {\mathfrak g}\rightarrow {\mathfrak g}^*$ such that

and then our systems will automatically verify the second law of thermodynamics where the Casimir function $C$ will play the role of the entropy.

Given an arbitrary positive semidefinite scalar product on ${\mathfrak g}$

we can define ${\mathcal F}$ by

for all $\eta\in {\mathfrak g}$.

With this definition it is obvious that

Remark 3.3. The force term ${\mathcal F}$ given in 3.7 is only a proposal that generalizes the construction given for the metriplectic rigid body. An interesting application of this construction is for the design of controlled Euler-Poincaré-equations preserving the metriplectic properties. Moreover, the map in (3.6) is related to the covariant 3-tensor

which, in the case of the Killing form ${\mathcal K}(\xi, \eta) = -\hbox{trace} (ad_\xi\circ ad_\eta)$, is skew-symmetric and this is related with the construction given by ^[2].

4.   Generic integrators

In this section, we will derive a second order integrator preserving some of the properties of a metriplectic system. The typical methods for the type of thermodynamics evolution equations that we are studying in this paper are known as generic integrators (GENERIC = general equation for the non-equilibrium reversible–irreversible coupling, see ^[14,15]). The methods that we are proposing in this paper are of the generic type since our construction is based on the discrete gradient methods that are typically used for systems defined by an almost-Poisson bracket, and, in this case, the methods guarantee the exact preservation of the energy and good behavior of the entropy production. We will start with the classical methods where $P = {\mathbb R}^n$, and after this, we will discuss the case of $P$ being a general differentiable manifold.

4.1.   Discrete gradient systems

For ODEs in skew-gradient form, i.e., $\dot{x} = \Pi(x)\nabla H(x)$ with $x\in \mathbb{R}^n$ and $\Pi(x)$ a skew-symmetric matrix, it is immediate to check that $H$ is a first integral. Indeed,

due to the skew-symmetry of $\Pi$. Using discretizations of the gradient $\nabla H(x)$, it is possible to define a class of integrators which preserve the first integral $H$ exactly.

Definition 4.1. Let $H:\mathbb{R}^n\longrightarrow \mathbb{R}$ be a differentiable function. Then, $\bar{\nabla}H:\mathbb{R}^{2n}\longrightarrow \mathbb{R}^n$ is a discrete gradient of $H$ if it is continuous and satisfies

Some well-known examples of discrete gradients are:

● The mean value (or averaged) discrete gradient introduced in ^[16] and given by

● The midpoint (or Gonzalez) discrete gradient, introduced in ^[17] and given by

● The coordinate increment discrete gradient, introduced in ^[18], with each component given by

when $x_i'\not = x_i$, and $\bar{\nabla}_{3}H(x, x')_i = \frac{\partial H}{\partial x_i}(x'_1, \ldots, x'_{i-1}, x'_i = x_{i}, x_{i+1}, \ldots, x_n)$ otherwise.

4.2.   Construction of Metriplectic or Generic integrators

The idea is to construct a geometric integrator preserving as much as possible the properties of the continuous metriplectic Euler-Poincaré equations and, in particular, preserving the two laws of thermodynamics. We are in the category of generic integrators ^[14,15] since we will use a discretization of the differential of $H$ using a discrete gradient, and a discretization of the positive semi-definite inner product ${\mathcal K}$.

Consider a Gonzalez' discrete gradient $\bar{\nabla}_2H: {\mathbb R}^{2n}\rightarrow {\mathbb R}^n$, the Poisson tensor $\Pi(z)$ where $z = \frac{x+x'}{2}$, and a discretization ${\mathcal K}_d$ of the inner product ${\mathcal K}$ which is also positive semi-definite. Then the generic integrator is constructed as a discretization of equation (2.3) as follows:

For any $x\in P$, we assume that the numerical scheme (4.5) generates a local evolution in a neighborhood $U$ of $x$, in the sense that there exist real numbers $\bar{h}, T > 0$, and a discrete flow map $\varphi: U\times [0, \bar{h}]\rightarrow P$ such that for any $x_0\in U$ and $h\in [0, \bar{h}]$, the sequence $\{x_k\}$ generated by

satisfies equation (4.5) for all $k$ such that $kh\in [0, T]$.

Proposition 4.2. [Second law] The generic integrator verifies

and $k(x_{k+1/2})\geq 0$, where $x_{k+1/2} = \frac{x_k+x_{k+1}}{2}$.

Proof. Using Taylor's expansion we have that

Remark 4.3. Observe that for an arbitrary second order integrator, we will uniquely obtain that $S(x_{k+1})-S(x_k) = \bar{k}(x_{k+1/2})h+ {\mathcal O}(h^3)$, but in general $\bar{k}(x_{k+1/2})$ is not necessarily always positive. Therefore, it is crucial to discretize ${\mathcal K}$ while maintaining semi-definite positiveness.

Now, for the exact preservation of the energy it is necessary to construct a discretization ${\mathcal K}_d$ of ${\mathcal K}$ given in (2.2) such that $\bar{\nabla}_2H$ is an element of the kernel of ${\mathcal K}_d$.

As in (2.2), we consider

With the semi-definite positive inner product (4.6), we deduce the following.

Proposition 4.4. [First law] The generic integrator preserves exactly the energy function $H$, that is,

Proof.

since $K_d(x_{k+1/2})\bar{\nabla}_2H(x_k, x_{k+1}) = 0$.

Remark 4.5. Observe that Proposition 4.2 guarantees the correct behavior of the entropy, avoiding the numerical dissipation which occurs with arbitrary general methods ^[19]. Moreover, the conditions for deriving this metriplectic integrator are less restrictive than the usual ones for constructing generic integrators since, in this last case, it is necessary to also construct a discrete gradient $\bar{\nabla} C(x, x')$ for the entropy verifying the property that $\tilde{\Pi}(z)\bar{\nabla} C(x, x') = 0$, where $\tilde{\Pi}$ is a discretization of the Poisson tensor.

Next, we will study the example of the relaxed rigid body to show the constructability of our method but of course it is possible to apply to various thermodynamical examples as in ^[15] and to coupled thermodynamical systems ^[10].

4.3.   Example: Numerical integration of the relaxing rigid body

The rigid body equations are given by

These equations are the Euler-Poincaré equations for the Lagrangian $l: {\mathfrak g}\rightarrow {\mathbb R}$

Now, using the Legendre transformation, we define the associated momenta:

Then, the equations of motion of the system become

This is a Lie-Poisson system, and the equations are written in matrix form as

where

and $H(p_1, p_2, p_3) = \frac{1}{2}\left(\frac{p_1^2}{I_1}+\frac{p_2^2}{I_2}+\frac{p_3^2}{I_3}\right)$. Consider now the positive semi-definite inner product defined by

where ${\mathcal G}$ is the canonical metric of ${\mathbb R}^3$. After some straightforward computations, we derive that ${\mathcal K}$ is defined by the matrix

The entropy is defined by the Casimir function

and the dynamics of the metriplectic system is given by

(see ^[20] for numerical integration of Lie-Poisson systems using the midpoint rule without the dissipative term) or

From construction, we get $\Pi\nabla S = 0$ and $K\nabla H = 0$.

Using the notation $(P_1, P_2, P_3) = \varphi(p_1, p_2, p_3, h)$, the generic integrator is constructed taking

and the discrete semi-definite scalar product

The metriplectic integrator is given in this case by the midpoint rule:

In this case, since the Casimir is quadratic, we have also that

the system verifies that

and also $H(P_1, P_2, P_3) = H(p_1, p_2, p_3)$ for the discrete flow $\varphi(p_1, p_2, p_3, h) = (P_1, P_2, P_3)$. Then, in this case, the midpoint method preserves the energy exactly and, moreover, the entropy production rate has the correct behavior.

Remark 4.6. Observe that in this case, the numerical method corresponds exactly to the midpoint discretization because the Hamiltonian is quadratic, leading to results similar to those in ^[20]. Of course, the method can be applied to general metriplectic systems, and our results would differ from the standard midpoint discretization.

Remark 4.7. The previous construction is based on the discretization ${\mathcal K}_d$ of ${\mathcal K}$ given in (2.2). The proposed method can be applied to a general metriplectic system (see Section 2.3):

where ${\mathcal K}$ is a semi-definite inner product with $\nabla H\in \ker {\mathcal K}$. Then take an arbitrary complement ${\mathcal D}$ (assuming local constant rank) such that

and consider the positive definite inner product ${\mathcal G}^*$ by

where $\{Y_i\}, 1\leq i\leq n-l$ and $\{Z_j\}, 1\leq k\leq l$ are bases of ${\mathcal D}$ and $\ker {\mathcal K}$, respectively, and where $\nabla H = Z_1$. Obviously, ${\mathcal G}^*$ is positive definite, and if we denote by ${\mathcal P}_{{\mathcal G}^*}$ the orthogonal projector onto ${\mathcal D}$, then

for all $X_1, X_2\in T^*P$. Now, given a discrete gradient $\widetilde{\nabla} H$, it is only necessary to take the new decomposition

where a basis of $\widetilde{\ker {\mathcal K}}$ is now given by $\{\widetilde{\nabla} H, Z_2, \ldots, Z_l\}$ and $\widetilde{\ker {\mathcal K}}^{\perp}$ is the corresponding ${\mathcal G}^*$-orthogonal complement. If we denote the modified orthogonal projector $\widetilde{{\mathcal P}_{{\mathcal G}^*}}$ onto $\widetilde{\ker {\mathcal K}}^{\perp}$, then the corresponding semi-definite inner product with $\widetilde{\nabla} H$ in its kernel is precisely

4.4.   Extension to differentiable manifolds

We can extend this construction to the case where we are working on $P$ a general manifold. To start, we will need to introduce a finite difference map or retraction map $R_h: U\subset TP\rightarrow P\times P$ and its inverse map $R^{-1}_h: \bar{U}\subset P\times P\rightarrow TP$ ^[21]. For any $(x, x')\in \bar{U}$ we denote by $z = \tau_P(R_h^{-1}(x, x'))\in TP$. We can use a type of retraction that is constructed using an auxiliary Riemannian metric ${\mathcal G}$ on $P$ with associated geodesic spray $\Gamma_{\mathcal G}$ ^[22]. The associated Riemannian exponential for a small enough $h > 0$ is constructed as

where we have the standard exponential map on a Riemannian manifold defined by

where $t\rightarrow \gamma_v(t)$ is the unique geodesic such that $\gamma_v'(0) = v$. Another interesting possibility related to the midpoint rule is

Both maps are local diffeomorphisms, and then we can consider the corresponding inverse maps that we generically denote by $R^{-1}_h$ at the beginning of this section.

Define a discrete gradient as a map $\bar{\nabla}{ H}: \bar{U}\subseteq P\times P \longrightarrow T^{*}P$ such that the following diagram commutes

and verifies the following two properties:

In the case, when we have a Riemannian metric ${\mathcal G}$ on $P$, we construct the following midpoint discrete gradient

where $\flat_{\mathcal G}:TP\rightarrow T^*P$ is given by $\flat_{\mathcal G}(u)(v) = {\mathcal G}(u, v)$ for $u, v\in TP$ and $z = \tau_P(R^{-1}_h(x, x'))\in P$.

The metriplectic integrator that we propose is written as

where $z = \tau_P(R^{-1}_h(x_k, x_{k+1}))$ and ${\mathcal K}_d$ is constructed as in (4.6).

Author contribution

All the authors contribute equally to this work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

Anthony Bloch was partially supported by NSF grant DMS-2103026, and AFOSR grants FA 9550-22-1-0215 and FA 9550-23-1-0400. Marta Farré Puiggal acknowledges financial support from FWO grant 1282021N. David Martín de Diego acknowledges financial support from the Spanish Ministry of Science and Innovation, under grants PID2022-137909NB-C21, TED2021-129455B-I00 and CEX2019-000904-S funded by MCIN/AEI/10.13039/501100011033 and BBVA Foundation via the project "Mathematical optimization for a more efficient, safer and decarbonized maritime transport". We also thank the reviewers for useful remarks which improved the exposition.

Conflict of interest

The authors declare there is no conflict of interest.