Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Metriplectic Euler-Poincaré equations: smooth and discrete dynamics

  • In this paper we will introduce a discrete version of systems obtained by modifications of the Euler-Poincaré equations when we add a special type of dissipative force, so that the equations of motion can be described using the metriplectic formalism. The metriplectic representation of the dynamics allows us to describe the conservation of energy, as well as to guarantee entropy production. For deriving the discrete equations we use discrete gradients to numerically simulate the evolution of the continuous metriplectic equations preserving their main properties: preservation of energy and correct entropy production rate.

    Citation: Anthony Bloch, Marta Farré Puiggalí, David Martín de Diego. Metriplectic Euler-Poincaré equations: smooth and discrete dynamics[J]. Communications in Analysis and Mechanics, 2024, 16(4): 910-927. doi: 10.3934/cam.2024040

    Related Papers:

    [1] Isaac Neal, Steve Shkoller, Vlad Vicol . A characteristics approach to shock formation in 2D Euler with azimuthal symmetry and entropy. Communications in Analysis and Mechanics, 2025, 17(1): 188-236. doi: 10.3934/cam.2025009
    [2] Jonas Schnitzer . No-go theorems for $ r $-matrices in symplectic geometry. Communications in Analysis and Mechanics, 2024, 16(3): 448-456. doi: 10.3934/cam.2024021
    [3] Shuyue Ma, Jiawei Sun, Huimin Yu . Global existence and stability of temporal periodic solution to non-isentropic compressible Euler equations with a source term. Communications in Analysis and Mechanics, 2023, 15(2): 245-266. doi: 10.3934/cam.2023013
    [4] Henryk Żołądek . Normal forms, invariant manifolds and Lyapunov theorems. Communications in Analysis and Mechanics, 2023, 15(2): 300-341. doi: 10.3934/cam.2023016
    [5] Yao Sun, Pan Wang, Xinru Lu, Bo Chen . A boundary integral equation method for the fluid-solid interaction problem. Communications in Analysis and Mechanics, 2023, 15(4): 716-742. doi: 10.3934/cam.2023035
    [6] Yining Yang, Cao Wen, Yang Liu, Hong Li, Jinfeng Wang . Optimal time two-mesh mixed finite element method for a nonlinear fractional hyperbolic wave model. Communications in Analysis and Mechanics, 2024, 16(1): 24-52. doi: 10.3934/cam.2024002
    [7] Panyu Deng, Jun Zheng, Guchuan Zhu . Well-posedness and stability for a nonlinear Euler-Bernoulli beam equation. Communications in Analysis and Mechanics, 2024, 16(1): 193-216. doi: 10.3934/cam.2024009
    [8] Vladimir Rovenski . Generalized Ricci solitons and Einstein metrics on weak $ K $-contact manifolds. Communications in Analysis and Mechanics, 2023, 15(2): 177-188. doi: 10.3934/cam.2023010
    [9] Zhen Wang, Luhan Sun . The Allen-Cahn equation with a time Caputo-Hadamard derivative: Mathematical and Numerical Analysis. Communications in Analysis and Mechanics, 2023, 15(4): 611-637. doi: 10.3934/cam.2023031
    [10] Chen Yang, Chun-Lei Tang . Sign-changing solutions for the Schrödinger-Poisson system with concave-convex nonlinearities in $ \mathbb{R}^{3} $. Communications in Analysis and Mechanics, 2023, 15(4): 638-657. doi: 10.3934/cam.2023032
  • In this paper we will introduce a discrete version of systems obtained by modifications of the Euler-Poincaré equations when we add a special type of dissipative force, so that the equations of motion can be described using the metriplectic formalism. The metriplectic representation of the dynamics allows us to describe the conservation of energy, as well as to guarantee entropy production. For deriving the discrete equations we use discrete gradients to numerically simulate the evolution of the continuous metriplectic equations preserving their main properties: preservation of energy and correct entropy production rate.



    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 ˙S0. 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 ˙x=Π(x)H(x) with xRn and Π(x) a skew-symmetric matrix (not necessarily associated to a Poisson bracket). Using a discrete gradient ˉ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 xx=ˉΠ(x,x)ˉH(x,x) preserving the energy H exactly, i.e. H(x)=H(x). Here ˉΠ(x,x) is a skew-symmetric matrix approximating Π(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.

    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.

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

    C(P)×C(P)C(P)(f,g){f,g}

    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,hC(P).

    Given a Poisson manifold with bracket {,} and a function fC(P) we may associate with f a unique vector field XfX(P), the Hamiltonian vector field defined by Xf(g)={g,f}.

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

    {f,g}:=Π(df,dg),f,gC(P).

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

    Take coordinates (xi), 1idimP=m, and let Πij be the components of the Poisson bivector, that is,

    Πij={xi,xj},Π=Πijxixj.

    Then if f,gC(P) we have

    {f,g}=mi,j=1{xi,xj}fxigxj=mi,j=1Πijfxigxj,

    and the Hamiltonian vector field is written in local coordinates as

    Xf=Πijfxjxi.

    Observe that the m×m matrix (Πij) verifies the following properties:

    (ⅰ) Skew-symmetry, Πij=Πji

    (ⅱ) Jacobi identity,

    ml=1(ΠilΠjkxl+ΠklΠijxl+ΠjlΠkixl)=0,i,j,k=1,,m.

    Define Π:TPTP by

    Π(α)=ιαΠ=Π(,α),

    where αTP, and β,ιαΠ=Π(α,β) for all βTP. The rank of Π at pP is exactly the rank of (Π)p:TpPTpP. Because of the skew-symmetry of Π, we know that the rank of Π at a point pP is an even integer.

    Given a function HC(P), a Hamiltonian function, we have the corresponding Hamiltonian vector field:

    XH=Π(dH).

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

    dxdt(t)=XH(x(t)). (2.1)

    We say that a function fC(P) is a first integral of the Hamiltonian vector field XH if for any solution x(t) of Equation (2.1) we have

    dfdt(x(t))=0.

    In other words, if XH(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,Π) a function CC(P) is called a Casimir function of Π if XC=0, that is, if {C,g}=0 for all gC(P).

    Assume that for each point xP we have a positive semi-definite inner product for co-vectors

    Kx:TxP×TxPR

    from which we can define K:TPTP by

    K(αx)=Kx(αx,)

    and the gradient vector field

    gradKS=K(dS)

    for any function S:PR.

    K defines a symmetric bracket given by

    (df,dg)=K(df,dg).

    Take coordinates (xi), 1idimP=m, and let Kij be the components of the inner product given by

    Kij=(xi,xj).

    Then if f,gC(P), the symmetric bracket is expressed as

    (f,g)=mi,j=1(xi,xj)fxigxj=mi,j=1Kijfxigxj.

    Observe that the m×m matrix (Kij) verifies Kij=Kji and all the eigenvalues are positive or zero.

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

    G:TxP×TxPR

    defined by G(df,dg)=G(gradGf,gradGg)=gradGf(g).

    Since we are interested in defining a semi-definite inner product K such that K(dH,)=0 then we define

    K(df,dg)=1G(dH,dH)G(G(dH,dH)dfG(dH,df)dH,G(dH,dH)dgG(dH,dg)dH)=G(dH,dH)G(df,dg)G(dH,df)G(dH,dg). (2.2)

    In coordinates, we have

    Kij=CHgijgiˉjHxˉjgˉijHxˉi,

    where CH=gijHxiHxj, (gij) are the components of the Riemannian metric in a given coordinate system and (gij) denotes its inverse matrix.

    By construction, K is positive semi-definite and K(dH,)=0.

    Remark 2.1. Additionally we can add new functions La:PR, 1aN, to this construction in such a way that K(dLa,)=0, considering

    K(df,dg)=G(dfCabG(dLa,df)dLb,dgCabG(dLa,dg)dLb)

    where Cab=G(dLa,dLb), 1aN and L1=H.

    A metriplectic system consists of a smooth manifold P, two smooth vector bundle maps Π,K:TPTP verifying that πP=τPΠ and πP=τPK (where τP:TPP and πP:TPP are the canonical projections), and two functions H,SC(P) called the Hamiltonian (or total energy) and the entropy of the system, such that for all f,gC(P):

    {f,g}=df,Π(dg) is a Poisson bracket (Π denotes the Poisson bi-vector).

    (f,g)=df,K(dg) is a positive semi-definite symmetric bracket, i.e., (,) is bilinear and symmetric.

    K(dH)=0 or, equivalently, (H,f)=0, fC(P).

    Π(dS)=0 or, equivalently, {S,f}=0, fC(P), that is, S is a Casimir function for the Poisson bracket.

    Consider the function E=H+S:PR. Then, the dynamics of the metriplectic system is determined by

    dxdt=Π(dE(x(t)))+K(dE(x(t)))=Π(dH(x(t)))+K(dS(x(t)))=XH(x(t))+gradKS(x(t)),

    where XH=Π(dH) and gradKS=K(dS). From the equations of motion, it is simple to deduce the following:

    ● First law: conservation of energy, dHdt={H,H}+(H,S)=0

    ● Second law: entropy production, dSdt=(S,S)0.

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

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

    ˙xi=ΠijHxj+KijSxj,1in

    or, in matrix form, as

    ˙x=ΠH+KS,1in. (2.3)

    Let Φ:G×PP be a smooth (left) action of a Lie group G on P, given by Φ(g,x)=Φg(x)=gx with gG and xP. Denote by g the corresponding Lie algebra. The action satisfies the following properties:

    Φ(e,x)=x where e is the neutral element of G;

    ● For every g1,g2G and for every xP

    Φ(g1,Φ(g2,x))=Φ(g1g2,x).

    The infinitesimal generator of the action corresponding to a Lie algebra element ξg is the vector field ξP on P given by

    ξP(x)=ddt|t=0(exp(ξt)x).

    Let P be a Poisson manifold with Poisson bracket {,} and assume that the action Φ is a Poisson action, that is,

    Φg{f,h}={Φgf,Φgh},f,hC(P)gG.

    A momentum map for the action Φ is a smooth map J:Pg such that for each ξg, the associated map Jξ:PR defined by Jξ(x)=J(x),ξ satisfies that XJξ=ξP for all ξg where XJξ(f)={f,Jξ}. As a consequence, for any function fC(P)

    {f,Jξ}=ξP(f).

    If the Lie algebra g acts on the Poisson manifold P and admits a momentum map J:Pg, and if HΦg=H (which is equivalent to ξP(H)=0 for all ξg assuming that G is connected), then Jξ is a constant of the motion of XH.

    Additionally, for the metriplectic system we will assume that

    (f,Jξ)=0,ξg and fC(P),

    or, equivalently, K(Jξ)=0. Then, for the metriplectic system we have

    dJξdt={Jξ,H}+(Jξ,S)=0

    and, therefore, Jξ:PR 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=TG, where G is a Lie group, and we consider as a left action Ψg=TLg:TGTG, where Lg:GG is the left action. Under the symmetry conditions, the system reduces to a metriplectic system on g, the dual of the Lie algebra of G.

    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:gR, where g is a Lie algebra, and its corresponding Euler-Poincaré equations [11,12]:

    ddt(δlδξ)=adξδlδξ, (3.1)

    where ξg and adξα,ξ=α,[ξ,ξ] for all ξg and αg. From this equation it is clear that the energy El=δlδξ,ξl of the system is preserved, that is,

    dEldt=ddt(δlδξ,ξl)=0.

    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:gg of the form

    F(ξ)=adξF(ξ),ξg

    where F:gg is an arbitrary map, then the forced Euler-Lagrange equations are

    ddt(δlδξ)=adξδlδξ+F=adξ[δlδξ+F]. (3.2)

    Assume that g is finite dimensional and {ea}, 1an=dimg is a basis of the Lie algebra with structure constants Cdab, that is,

    [ea,eb]=Cdabed,

    and denote by (ξa(t)) the coordinates of a curve ξ(t)g. Then, the equations (3.2) are

    ddt(δlδξb(ξ(t)))=Cdabξa(t)(δlδξd(ξ(t))+Fd(ξ(t))), (3.3)

    where F(ξ)=Fd(ξ)ed and {ea}, 1an, is the dual basis of {ea}.

    Example 3.1. In the case of G=SO(3) if we identify its Lie algebra g with R3 with the usual vector cross product then we have

    ddt(δlδΩ)=δlδΩ×Ω+F×Ω

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

    l(Ω1,Ω2,Ω3)=12(I1Ω21+I2Ω22+I3Ω23)

    then (3.2) are

    I1˙Ω1=(I2I3)Ω2Ω3+Ω3F2(Ω)Ω2F3(Ω),I2˙Ω2=(I3I1)Ω3Ω1+Ω1F3(Ω)Ω3F1(Ω),I3˙Ω3=(I1I2)Ω1Ω2+Ω2F1(Ω)Ω1F2(Ω),

    where F=(F1,F2,F3) and Ω=(Ω1,Ω2,Ω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

    ˙μ=adδH/δμ(μ+F(δHδμ)), (3.4)

    where H is defined by H(μ)=μ,ξ(μ)L(ξ(μ)) and μ=δlδξ(ξ). This is a particular case of the forced Lie-Poisson equations [6].

    Now, if C:gR is a Casimir function for the Lie-Poisson bracket of g, then along the evolution of the system (3.4), we have

    dCdt=F(δHδμ),[δHδμ,δCδμ]. (3.5)

    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

    H(Π1,Π2,Π3)=12(Π21I1+Π22I2+Π23I3),C(Π1,Π2,Π3)=12(Π21+Π22+Π23),

    in induced coordinates (Π1=I1Ω1,Π2=I2Ω2,Π3=I3Ω3) on (R3)R3, Equation (3.5) is

    dCdt=(1I21I3)Π2Π3F1+(1I31I1)Π1Π3F2+(1I11I2)Π1Π2F3.

    For instance, if we take F:gR3gR3 as

    F(Ω)=((I3I2)Ω2Ω3,(I1I3)Ω1Ω3,(I2I1)Ω1Ω2),

    then we get

    dCdt0.

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

    ˙Π1=(I2I3)I2I3Π2Π3+(I1I3)I1I23Π1Π23(I2I1)I1I22Π1Π22,˙Π2=(I3I1)I3I1Π3Π1+(I2I1)I21I2Π21Π2(I3I2)I2I23Π2Π23,˙Π3=(I1I2)I1I2Π1Π2+(I3I2)I22I3Π22Π3(I1I3)I21I3Π21Π3.

    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 F:gg such that

    dCdt0

    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 g

    K:g×gR (3.6)

    we can define F by

    F(ξ),η=K(η,[ξ,Cμ]) (3.7)

    for all ηg.

    With this definition it is obvious that

    dCdt=F(δHδμ),[δHδμ,δCδμ]=K([δHδμ,δCδμ],[δHδμ,Cμ])0.

    Remark 3.3. The force term 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

    cK:g×g×gR(ξ1,ξ2,ξ3)K(ξ1,[ξ2,ξ3])

    which, in the case of the Killing form K(ξ,η)=trace(adξadη), is skew-symmetric and this is related with the construction given by [2].

    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=Rn, and after this, we will discuss the case of P being a general differentiable manifold.

    For ODEs in skew-gradient form, i.e., ˙x=Π(x)H(x) with xRn and Π(x) a skew-symmetric matrix, it is immediate to check that H is a first integral. Indeed,

    ˙H=H(x)T˙x=H(x)TΠ(x)H(x)=0,

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

    Definition 4.1. Let H:RnR be a differentiable function. Then, ˉH:R2nRn is a discrete gradient of H if it is continuous and satisfies

    ˉH(x,x)T(xx)=H(x)H(x),forallx,xRn, (4.1a)
    ˉH(x,x)=H(x),forallxRn. (4.1b)

    Some well-known examples of discrete gradients are:

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

    ˉ1H(x,x):=10H((1ξ)x+ξx)dξ, for xx. (4.2)

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

    ˉ2H(x,x):=H(12(x+x))+H(x)H(x)H(12(x+x))T(xx)|xx|2(xx), for xx. (4.3)

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

    ˉ3H(x,x)i:=H(x1,,xi,xi+1,,xn)H(x1,,xi1,xi,,xn)xixi,1in, (4.4)

    when xixi, and ˉ3H(x,x)i=Hxi(x1,,xi1,xi=xi,xi+1,,xn) otherwise.

    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 K.

    Consider a Gonzalez' discrete gradient ˉ2H:R2nRn, the Poisson tensor Π(z) where z=x+x2, and a discretization Kd of the inner product K which is also positive semi-definite. Then the generic integrator is constructed as a discretization of equation (2.3) as follows:

    xxh=Π(z)ˉ2H(x,x)+Kd(z)S(z) (4.5)

    For any xP, 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 ˉh,T>0, and a discrete flow map φ:U×[0,ˉh]P such that for any x0U and h[0,ˉh], the sequence {xk} generated by

    xk=φ(xk1,h)=φk(x0,h)

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

    Proposition 4.2. [Second law] The generic integrator verifies

    S(xk+1)S(xk)=k(xk+1/2)h+O(h3)wherexk+1=φ(xk,h)

    and k(xk+1/2)0, where xk+1/2=xk+xk+12.

    Proof. Using Taylor's expansion we have that

    S(xk+1)S(xk)+O(|xk+1xk|3)=S(xk+1/2)T(xk+1xk)=hS(xk+1/2)TΠ(xk+1/2)ˉ2H(xk,xk+1)+hS(xk+1/2)TKd(xk+1/2)S(xk+1/2)=hS(xk+1/2)TKd(xk+1/2)S(xk+1/2)0.

    Remark 4.3. Observe that for an arbitrary second order integrator, we will uniquely obtain that S(xk+1)S(xk)=ˉk(xk+1/2)h+O(h3), but in general ˉk(xk+1/2) is not necessarily always positive. Therefore, it is crucial to discretize K while maintaining semi-definite positiveness.

    Now, for the exact preservation of the energy it is necessary to construct a discretization Kd of K given in (2.2) such that ˉ2H is an element of the kernel of Kd.

    As in (2.2), we consider

    Kd(df,dg)=G(ˉ2H,ˉ2H)[G(ˉ2H,ˉ2H)G(df,dg)G(ˉ2H,df)G(ˉ2H,dg)]. (4.6)

    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,

    H(xk+1)H(xk)=0.

    Proof.

    H(xk+1)H(xk)=ˉ2HT(xk,xk+1)(xk+1xk)=hˉ2HT(xk,xk+1)Π(xk+1/2)ˉ2H(xk,xk+1)+hˉ2HT(xk,xk+1)Kd(xk+1/2)S(xk+1/2)=0

    since Kd(xk+1/2)ˉ2H(xk,xk+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 ˉC(x,x) for the entropy verifying the property that ˜Π(z)ˉC(x,x)=0, where ˜Π 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].

    The rigid body equations are given by

    I1˙Ω1=(I2I3)Ω2Ω3,I2˙Ω2=(I3I1)Ω1Ω3,I3˙Ω3=(I1I2) Ω1Ω2.

    These equations are the Euler-Poincaré equations for the Lagrangian l:gR

    l(Ω1,Ω2,Ω3)=12I1Ω21+I2Ω22+I3Ω23.

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

    p1=lΩ1=I1Ω1,p2=lΩ2=I2Ω2,p3=lΩ3=I3Ω3.

    Then, the equations of motion of the system become

    ˙p1=I2I3I2I3p2p3,˙p2=I3I1I1I3p1p3,˙p3=I1I2I1I2p1p2.

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

    ˙p=ΠH=Π(dH),

    where

    Π=(0I3p3I2p2I3p30I1p1I2p2I1p10)

    and H(p1,p2,p3)=12(p21I1+p22I2+p23I3). Consider now the positive semi-definite inner product defined by

    K(df,dg)=[G(dH,dH)G(df,dg)G(dH,df)G(dH,dg)].

    where G is the canonical metric of R3. After some straightforward computations, we derive that K is defined by the matrix

    K=(p22I22+p23I23p1p2I1I2p1p3I1I3p1p2I1I2p21I21+p23I23p2p3I2I3p1p3I1I3p2p3I2I3p21I21+p22I22).

    The entropy is defined by the Casimir function

    S(p1,p2,p3)=12(p21+p22+p23)

    and the dynamics of the metriplectic system is given by

    ˙p=ΠH+KS,

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

    ˙p1=I2I3I2I3p2p3+(1I221I1I2)p1p22+(1I231I1I3)p1p23,˙p2=I3I1I1I3p1p3+(1I211I1I2)p2p21+(1I231I2I3)p2p23,˙p3=I1I2I1I2p1p2+(1I211I1I3)p3p21+(1I221I2I3)p22p3.

    From construction, we get ΠS=0 and KH=0.

    Using the notation (P1,P2,P3)=φ(p1,p2,p3,h), the generic integrator is constructed taking

    ˉ2H(P1+p12,P2+p22,P3+p32)=(P1+p12I1,P2+p22I2,P3+p32I3)=(z1/I1,z2/I2,z3/I3)

    and the discrete semi-definite scalar product

    Kd=(z22I22+z23I23z1z2I1I2z1z3I1I3z1z2I1I2z21I21+z23I23z2z3I2I3z1z3I1I3z2z3I2I3z21I21+z22I22).

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

    P1p1h=I2I3I2I3z2z3+(1I221I1I2)z1z22+(1I231I1I3)z1z23,P2p2h=I3I1I1I3z1z3+(1I211I1I2)z2z21+(1I231I2I3)z2z23,P3p3h=I1I2I1I2z1z2+(1I211I1I3)z3z21+(1I221I2I3)z22z3.

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

    S=ˉ2S,

    the system verifies that

    S(P1,P2,P3)S(p1,p2,p3)0,

    and also H(P1,P2,P3)=H(p1,p2,p3) for the discrete flow φ(p1,p2,p3,h)=(P1,P2,P3). Then, in this case, the midpoint method preserves the energy exactly and, moreover, the entropy production rate has the correct behavior.

    Figure 1.  I1=10,I2=5,I3=1,h=0.1, initial conditions p1=0.001,p2=1,p3=0.001.

    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 Kd of K given in (2.2). The proposed method can be applied to a general metriplectic system (see Section 2.3):

    dxdt=Π(x(t))H(x(t)))+K(x(t))S(x(t)

    where K is a semi-definite inner product with HkerK. Then take an arbitrary complement D (assuming local constant rank) such that

    TP=DkerK

    and consider the positive definite inner product G by

    G(Yi,Yi)=K(Yi,Yi),G(Yi,Zj)=0,G(Zj,Zj)=δjj

    where {Yi},1inl and {Zj},1kl are bases of D and kerK, respectively, and where H=Z1. Obviously, G is positive definite, and if we denote by PG the orthogonal projector onto D, then

    K(X1,X2)=G(PG(X1),PG(X2))

    for all X1,X2TP. Now, given a discrete gradient ˜H, it is only necessary to take the new decomposition

    TP=~kerK~kerK (4.7)

    where a basis of ~kerK is now given by {˜H,Z2,,Zl} and ~kerK is the corresponding G-orthogonal complement. If we denote the modified orthogonal projector ~PG onto ~kerK, then the corresponding semi-definite inner product with ˜H in its kernel is precisely

    Kd(X1,X2)=G(~PG(X1),~PG(X2))

    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 Rh:UTPP×P and its inverse map R1h:ˉUP×PTP [21]. For any (x,x)ˉU we denote by z=τP(R1h(x,x))TP. We can use a type of retraction that is constructed using an auxiliary Riemannian metric G on P with associated geodesic spray ΓG [22]. The associated Riemannian exponential for a small enough h>0 is constructed as

    exph(v)=(τQ(v),expτQ(v)(hv)),

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

    expτQ(v)(v)=γv(1),

    where tγv(t) is the unique geodesic such that γv(0)=v. Another interesting possibility related to the midpoint rule is

    ~exph(v)=(expτQ(v)(hv/2),expτQ(v))(hv/2)). (4.8)

    Both maps are local diffeomorphisms, and then we can consider the corresponding inverse maps that we generically denote by R1h at the beginning of this section.

    Define a discrete gradient as a map ˉH:ˉUP×PTP such that the following diagram commutes

    and verifies the following two properties:

    ˉH(x,x),R1h(x,x)=H(x)H(x), for all (x,x)ˉU, (4.9a)
    ˉH(x,x)=dH(x), for all xP. (4.9b)

    In the case, when we have a Riemannian metric G on P, we construct the following midpoint discrete gradient

    ˉ2H(x,x):=dH(z)+H(x)H(x)dH(z)(R1h(x,x))G(R1h(x,x),R1h(x,x))G(R1h(x,x)), for xx, (4.10)

    where G:TPTP is given by G(u)(v)=G(u,v) for u,vTP and z=τP(R1h(x,x))P.

    The metriplectic integrator that we propose is written as

    R1h(xk,xk+1)=Π(zτ)ˉ2H(xk,xk+1)+Kd(z)C(z)

    where z=τP(R1h(xk,xk+1)) and Kd is constructed as in (4.6).

    All the authors contribute equally to this work.

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

    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.

    The authors declare there is no conflict of interest.



    [1] P. J. Morrison, A paradigm for joined Hamiltonian and dissipative systems, Phys. D, 18 (1986), 410–419. https://doi.org/10.1016/0167-2789(86)90209-5 doi: 10.1016/0167-2789(86)90209-5
    [2] A. M. Bloch, P. J. Morrison, T. S. Ratiu, Gradient flows in the normal and Kähler metrics and triple bracket generated metriplectic systems, in Recent trends in dynamical systems, Springer Basel, (2013), 371–415. https://doi.org/10.1007/978-3-0348-0451-6_15
    [3] A. M. Bloch, P. J. Morrison, Algebraic structure of the plasma quasilinear equations, Phys. Lett. A, 88 (1982), 405–406. https://doi.org/10.1016/0375-9601(82)90664-8 doi: 10.1016/0375-9601(82)90664-8
    [4] M. Grmela, H. C.Öttinge, Dynamics and thermodynamics of complex fluids. Ⅰ. Development of a general formalism, Phys. Rev. E, 56 (1987), 6620–6632. https://doi.org/10.1103/PhysRevE.56.6620 doi: 10.1103/PhysRevE.56.6620
    [5] P. J. Morrison, M. H. Updike, Inclusive curvaturelike framework for describing dissipation: metriplectic 4-bracket dynamics, Phys. Rev. E, 109 (2024), 045202. https://doi.org/10.1103/physreve.109.045202 doi: 10.1103/physreve.109.045202
    [6] A. Bloch, P. S. Krishnaprasad, J. E. Marsden, T. S. Ratiu, The Euler-Poincaré equations and double bracket dissipation, Commun. Math. Phys., 175 (1996), 1–42. https://doi.org/10.1007/BF02101622 doi: 10.1007/BF02101622
    [7] O. Esen, G. Özcan, S. Sütlü, On extensions, Lie-Poisson systems, and dissipation, J. Lie Theory, 32 (2022), 327–382.
    [8] R. Abraham, J. E. Marsden, Foundations of Mechanics, AMS Chelsea Publishing, 1978.
    [9] P. Libermann, C. Marle, Symplectic geometry and analytical mechanics, volume 35 of Mathematics and its Applications. Springer Dordrecht, 1987. https://doi.org/10.1007/978-94-009-3807-6
    [10] I. Romero, Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics Part Ⅰ: Monolithic integrators and their application to finite strain thermoelasticity, Comput. Methods Appl. Mech. Engrg., 199 (2010), 1841–1858. https://doi.org/10.1016/j.cma.2010.02.014 doi: 10.1016/j.cma.2010.02.014
    [11] D. D. Holm, J. E. Marsden, T. S. Ratiu, The Euler–Poincaré equations and semidirect products with applications to continuum theories, Adv Math, 137 (1998), 1–81. https://doi.org/10.1006/aima.1998.1721 doi: 10.1006/aima.1998.1721
    [12] D. D. Holm, T. Schmah, C. Stoica, Geometric mechanics and symmetry: from finite to infinite dimensions, Oxford University Press, 2009.
    [13] M. Materassi, P. J. Morrison, Metriplectic torque for rotation control of a rigid body, Cybernetics and Physics, 7 (2018), 78–86. https://doi.org/10.35470/2226-4116-2018-7-2-78-86 doi: 10.35470/2226-4116-2018-7-2-78-86
    [14] H. C. Öttinger, GENERIC Integrators: Structure Preserving Time Integration for Thermodynamic Systems, J Non-Equil Thermody, 43 (2018), 89–100. https://doi.org/10.1515/jnet-2017-0034 doi: 10.1515/jnet-2017-0034
    [15] X. Shang, H. C. Öttinger, Structure-preserving integrators for dissipative systems based on reversible-irreversible splitting, Proc. A., 476 (2020), 20190446. https://doi.org/10.1098/rspa.2019.0446 doi: 10.1098/rspa.2019.0446
    [16] A. Harten, P. D. Lax, B. Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev., 25 (1983), 35–61. https://doi.org/10.1137/1025002 doi: 10.1137/1025002
    [17] O. Gonzalez, Time integration and discrete Hamiltonian systems, J. Nonlinear Sci., 6 (1996), 449–467. http://dx.doi.org/10.1007/s003329900018 doi: 10.1007/s003329900018
    [18] T. Itoh, K. Abe, Hamiltonian-conserving discrete canonical equations based on variational difference quotients, J. Comput. Phys., 76 (1988), 85–102. http://dx.doi.org/10.1016/0021-9991(88)90132-5 doi: 10.1016/0021-9991(88)90132-5
    [19] E. Hairer, C. Lubich, G. Wanner, Geometric numerical integratio, volume 31 of Springer Series in Computational Mathematics. Springer, Heidelberg, 2010.
    [20] M. A. Austin, P. S. Krishnaprasad, L. Wang, Almost Poisson integration of rigid body systems, J. Comput. Phys., 107 (1993), 105–117. https://doi.org/10.1006/jcph.1993.1128 doi: 10.1006/jcph.1993.1128
    [21] P. A. Absil, R. Mahony, R. Sepulchre, Optimization algorithms on matrix manifolds, Princeton University Press, 2008. https://doi.org/10.1515/9781400830244
    [22] E. Celledoni, S. Eidnes, B. Owren, T. Ringholm, Energy-preserving methods on Riemannian manifolds, Math. Comp., 89 (2020), 699–716. https://doi.org/10.1090/mcom/3470 doi: 10.1090/mcom/3470
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(632) PDF downloads(35) Cited by(0)

Figures and Tables

Figures(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog