# Conformal variational discretisation of infinite dimensional Hamiltonian systems with gradient flow dissipation

Damiano Lombardi, Cecilia Pagliantini

► To cite this version:

Damiano Lombardi, Cecilia Pagliantini. Conformal variational discretisation of infinite dimensional Hamiltonian systems with gradient flow dissipation. 2024. hal-04826188v2

HAL Id: hal-04826188

https://inria.hal.science/hal-04826188v2

Preprint submitted on 9 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Conformal variational discretisation of infinite dimensional Hamiltonian systems with gradient flow dissipation

Damiano Lombardi*

Cecilia Pagliantini†

## Abstract

Nonconservative evolution problems describe irreversible processes and dissipative effects in a broad variety of phenomena. Such problems are often characterised by a conservative part, which can be modelled as a Hamiltonian term, and a nonconservative part, in the form of gradient flow dissipation. Traditional numerical approximations of this class of problem typically fail to retain the separation into conservative and nonconservative parts hence leading to unphysical solutions. In this work we propose a mixed variational method that gives a semi-discrete problem with the same geometric structure as the infinite-dimensional problem. As a consequence the conservation laws and the dissipative terms are retained. Convergence results on the solution are established. Numerical tests of the Korteweg-de Vries equation and of the two-dimensional Navier-Stokes equations on the torus and on the sphere are presented to corroborate the theoretical findings.

## 1 Introduction

We consider time-dependent partial differential equations characterized by a conservative part, expressed through a Hamiltonian term, and a gradient flow dissipation. The problem is as follows: we look for $u$ in a Hilbert space $V$ such that

$$ \partial_t u = \mathcal{J}(u)\delta\mathcal{H}(u) + \mathcal{G}(u)\delta\mathcal{S}(u), \quad t > 0 \qquad (1.1) $$

supplied by the initial condition $u(0) = u_0 \in V$ and suitable boundary conditions encoded in $V$. In the evolution equation (1.1), $\mathcal{J}(u)$ characterizes the conservative term and is a Poisson operator, for any $u \in V$, while $\mathcal{G}(u)$ characterizes the dissipation and is a positive/negative semi-definite self-adjoint operator. The function $\mathcal{H}: V \to \mathbb{R}$ plays the role of Hamiltonian, while $\mathcal{S}: V \to \mathbb{R}$ is a produced/dissipated quantity. We postpone to Section 2 a detailed description of the model problem (1.1).

Systems of the form (1.1) describe irreversible processes and dissipative effects in a broad variety of phenomena. Examples are dissipative systems characterised by thermal conduction, electric

<footnote>*Sorbonne Université, Inria and CNRS, UMR 7598 Laboratoire Jacques-Louis Lions, Paris, France. (dami-ano.lombardi@inria.fr).</footnote>

<footnote>†Dipartimento di Matematica, Università di Pisa, Pisa, Italy. (cecilia.pagliantini@unipi.it).</footnote>

<footnote>C. Pagliantini acknowledges the MIUR Excellence Department Project awarded to the Department of Mathematics, University of Pisa, CUP I57G22000700001, and the INDAM/GNCS 2024 project CUP E53C23001670001.</footnote>
resistivity, viscous dissipation, and phenomena modelled by gradient flows [10, 15]. Problems (1.1) can also model phenomena where energy is conserved but entropy production and viscous heating are allowed. These processes are typical of a large variety of time-dependent systems out of thermodynamic equilibrium such as complex fluids [19], geophysical fluids [6], climate models [25], port-based networks [24], relativistic hydrodynamics and cosmology models [18], etc.. The latter class comprises also dissipative Hamiltonian flows and evolution problems ensuing from the Navier-Stokes equations, and kinetic models, like the Boltzmann equation and the Vlasov equation with collisions [11], where energy, mass, and momentum are conserved and entropy is produced. The formulation (1.1) includes also, as limit cases, Hamiltonian PDEs, obtained when $G(u)$ vanishes, and gradient flows, when $J(u)$ is zero.

In the numerical discretisation of problem (1.1) it is of crucial importance to retain the separation into a conservative part and a dissipative part to ensure a physically consistent approximation of the solution. It is by now well established that formulating numerical approximations compatible with the geometric structure underlying the model problem and not just approximating it yields numerical discretisations with superior stability properties. The geometric discretisation of canonical Hamiltonian PDEs has received considerable attention, see [13] for a review and references, and, similarly, structure-preserving approximations of dissipative evolution problems have been proposed both for finite dimensional [20] and infinite dimensional problems [5]. By contrast the numerical approximation of noncanonical Hamiltonian systems with gradient flow dissipation is, to the best of our knowledge, largely unexplored.

The bulk of geometric discretisations for problems of the form (1.1) has focused on specific examples such as the Navier-Stokes equations [21, 27] or a specific subclass, such as port-Hamiltonian systems [8]. A recent contribution to this topic is the finite element in time discretisation proposed in [2] where, although the geometric structure of the problem is not preserved, dissipation laws and conservation of general invariants are enforced via a particular choices of test functions.

In this work we derive a variational method for the numerical spatial approximation of evolution equations of the form (1.1) resulting in a semi-discrete problem with the same geometric structure as the continuous problem, in particular the separation into a conservative and a dissipative part. The approach we propose in this work is to first introduce a mixed formulation that, under natural regularity assumptions on the differential operators, is equivalent to the classical variational formulation of (1.1). Then a semi-discrete problem is obtained from a conformal spatial discretisation on a finite dimensional space $V_N \subset V$ of the aforementioned mixed variational formulation. The resulting approximation is shown to be composed of a conservative part characterized by a finite dimensional Poisson operator, and a dissipative part associated to a finite dimensional symmetric positive/negative semi-definite operator. The mixed formulation, at the expense of having two extra unknowns, enables a structure-preserving discretisation with classical low-order finite elements. Indeed, we prove that: (i) the discretisation of the invariants of motion of (1.1) yields invariants of motion of the semi-discrete problem; (ii) the kernel of the operators $J$ contained in the approximation space $V_N$ is preserved; (iii) the equilibria belonging to $V_N$ are equilibria of the semi-discrete problem; and (iv) the semi-discretisation of a metriplectic system, as defined in Section 2.3.1, remains metriplectic. In terms of accuracy of the approximation, we establish convergence of the solution of the semi-discrete system to the solution of the continuous mixed formulation.

The remainder of the paper is organised as follows. In Section 2 we introduce the mathematical framework of evolution equations of the form (1.1) by first characterizing the properties of Hamil-
tonian PDEs and, then, of gradient flows. Section 3 pertains to the derivation of the variational formulation of problem (1.1) and its mixed variational formulation. Section 4 is the core part of the paper where we introduce the proposed conformal variational discretisation based on the mixed formulation of the problem, we prove the preservation of the geometric structure of the problem and establish convergence results. In Section 5 we discuss the temporal discretisation of the semi-discrete problem via established symplectic methods and discrete gradients. Several numerical tests on the Korteweg-de Vries equation and on the 2D Navier-Stokes equation are presented in Section 6. Some concluding remarks are drawn in Section 7.

# 2 Hamiltonian systems with gradient flow dissipation

## 2.1 Infinite dimensional Hamiltonian dynamics

The definition of Hamiltonian partial differential equations require the introduction of function spaces endowed with a symplectic manifold structure and a characterization of the velocity field describing the Hamiltonian flow. To this end we recall some important concepts and refer to, e.g. [12], for a detailed introduction to the topic. First we introduce the concept of symplectic Hilbert scales.

**Definition 2.1 (Hilbert scale).** Let $X_0$ be a separable real Hilbert space with scalar product $\langle \cdot, \cdot \rangle_0$ and a Hilbert basis $\{\psi_k\}_{k \in \mathbb{Z}}$. Let $\{\theta_k\}_{k \in \mathbb{Z}}$ be a positive sequence such that $\theta_k \to \infty$ as $|k| \to \infty$. For any $s \in \mathbb{Z}$, let $X_s$ be a Hilbert space with norm and scalar product denoted by $\|\cdot\|_s$ and $\langle \cdot, \cdot \rangle_s$, respectively, and defined as

$$ \langle u, u \rangle_s = \|u\|_s^2 = \sum_{k \in \mathbb{Z}} |\langle u, \psi_k \rangle_0|^2 \theta_k^{2s}, \quad \forall u \in X_s. $$

The collection $\{X_s\}_{s \in \mathbb{Z}}$ is called Hilbert scale.

Note that $X_s$ is compactly embedded in $X_r$ if $s > r$, and is dense there. Moreover, the spaces $X_s$ and $X_{-s}$ are conjugated with respect to the scalar product $\langle \cdot, \cdot \rangle_0$, that is

$$ \|u\|_s = \sup\{\langle u, v \rangle_0 : v \in X_{-s} \cap X_0, \|v\|_{-s} = 1\} \quad \forall u \in X_s \cap X_0, $$

so that the scalar product $\langle \cdot, \cdot \rangle_0$ extends to a bilinear pairing $X_s \times X_{-s} \to \mathbb{R}$. With a small abuse of notation we will denote this pairing as $X_0$-scalar product.

For the characterization of a symplectic structure on Hilbert scales we need to introduce the concept of morphism of order $d$.

**Definition 2.2 (Morphism of order $d$).** Let $s_0, s_1 \in \mathbb{Z}$, $s_0 \le s_1$, be fixed and let $\{X_s\}_{s \in \mathbb{Z}}$ be a Hilbert scale. Given a linear map $L: \bigcap_{s \in \mathbb{Z}} X_s \to \bigcup_{s \in \mathbb{Z}} X_s$ we denote by $\|L\|_{L(X_{s_1}, X_{s_2})}$ its operator norm as a map $L: X_{s_1} \to X_{s_2}$. The map $L$ defines a morphism of order $d$ for $s \in [s_0, s_1]$ if $\|L\|_{L(X_s, X_{s-d})} < \infty$ for any $s \in [s_0, s_1]$. If the inverse map $L^{-1}$ exists and defines a morphism of order $-d$ for $s \in [s_0 + d, s_1 + d]$, then $L$ is an automorphism of order $d$ for $s \in [s_0, s_1]$. The adjoint $L^*: X_{d-s} \to X_{-s}$ of a morphism of order $d$ for $s \in [s_0, s_1]$ is a morphism of order $d$ for $s \in [d - s_1, d - s_0]$.
In particular, a bounded linear operator $L: X_{s_0} \to X_{s_0-d}$ can be regarded as a morphism of order $d$ for $s_0$.

Let $d \ge 0$ and $O_d$ be a non-empty open set in a Hilbert space $X_d$ from a Hilbert scale $\{X_s\}_{s \in \mathbb{Z}}$. For any $u \in O_d$ we identify the tangent space $T_u O_d$ with $X_d$, and treat differential $k$-forms as continuous functions

$$ \alpha_u : \bigwedge^k X_d \to \mathbb{R} $$

that are multi-linear and skew-symmetric. In particular, 1-forms can be written as $a(u)du(\xi) := \langle a(u), \xi \rangle_0$, for any $\xi \in X_d$, where $a: O_d \to X_{-d}$ and $d$ is the exterior derivative; while 2-forms can be written as $\omega_u = A(u)du \wedge du$ such that $\omega_u(\xi, \eta) := \langle A(u)\xi, \eta \rangle_0$, for any $\xi, \eta \in X_d$, where $A(u): X_d \to X_{-d}$ is a bounded self-adjoint operator. This implies that $\omega_u: X_\sigma \times X_\sigma \to \mathbb{R}$ defines a continuous skew-symmetric bilinear form if $\sigma \ge -d/2$.

Given a $C^k$-smooth function $\mathcal{F}: O_d \subset X_d \to \mathbb{R}$, $d \ge 1$, we consider the gradient map with respect to the pairing $\langle \cdot, \cdot \rangle$ as the $C^{k-1}$-smooth function

$$ \delta \mathcal{F}: O_d \to X_{-d} \quad \text{such that} \quad \langle \delta \mathcal{F}(u), \xi \rangle = d \mathcal{F}(u) \xi \quad \forall \xi \in X_d. \qquad (2.1) $$

A symplectic structure on the Hilbert scale can be introduced via the map $\mathcal{J}(u): X_{s+d_J} \to X_s$, with $u \in O_d$, a skew-adjoint morphism of order $d_J$ for $s \in [-d-d_J, d]$ called operator of the Poisson structure.

**Definition 2.3** (Poisson bracket). Let $\mathcal{F}, \mathcal{L}$ be $C^1$-smooth functions on $O_d$. Assume that they define gradient maps of orders $d_1$ and $d_2 \le 2d$ such that $d_1 + d_2 + d_J \le 2d$. Then, the Poisson bracket of $\mathcal{F}$ and $\mathcal{L}$ is the continuous (on $O_d$) function

$$ \{\mathcal{F}, \mathcal{L}\}(u) = \langle \mathcal{J}(u)\delta\mathcal{F}(u), \delta\mathcal{L}(u) \rangle_0 \quad \forall u \in O_d. $$

The Poisson bracket is skew-symmetric since $\mathcal{J}$ is skew-adjoint. Moreover, it satisfies the Leibniz rule $\{\mathcal{F}\mathcal{L}, \mathcal{H}\} = \{\mathcal{F}, \mathcal{H}\}\mathcal{L} + \mathcal{F}\{\mathcal{L}, \mathcal{H}\}$ and the Jacobi identity $\{\mathcal{F}, \{\mathcal{L}, \mathcal{H}\}\} + \{\mathcal{H}, \{\mathcal{F}, \mathcal{L}\}\} + \{\mathcal{L}, \{\mathcal{H}, \mathcal{F}\}\} = 0$, for any $\mathcal{F}, \mathcal{L}, \mathcal{H} \in C^2(O_d)$.

The set $O_d$ endowed with the bracket $\{\cdot, \cdot\}$ forms a Poisson manifold, while the space $C^\infty(O_d)$ of real-valued smooth functions over $(O_d, \{\cdot, \cdot\})$ together with the bracket $\{\cdot, \cdot\}$ forms a Lie algebra [1, Proposition 3.3.17], called the Poisson algebra of $O_d$.

Hilbert scales endowed with a Poisson structure provide the scaffolding for the definition of Hamiltonian dynamics in infinite dimension. Such dynamics is described by a Hamiltonian vector field as follows. By the bilinearity of $\{\cdot, \cdot\}$ and the Leibniz rule, given a $C^1$-smooth function $\mathcal{H}$ on $O_d$, there exists a locally unique vector field $V_{\mathcal{H}}$ defined as the continuous map $V_{\mathcal{H}}: O_d \to X_{-d-d_J}$ such that

$$ \forall u \in O_d \quad L_{V_{\mathcal{H}}(u)}\mathcal{F} = \{\mathcal{H}, \mathcal{F}\}(u), \quad \forall \mathcal{F} \in C^1(O_d), $$

where $L_X$ denotes the Lie derivative with respect to the velocity field $X$. The vector $V_{\mathcal{H}}(u)$ is called the Hamiltonian vector field of $\mathcal{H}$. Let $i_X$ denote the contraction operator with the vector field $X$. Since, for any $\mathcal{F} \in C^1(O_d)$,

$$ L_{V_{\mathcal{H}}(u)}\mathcal{F} = i_{V_{\mathcal{H}}(u)}d\mathcal{F} = d\mathcal{F}V_{\mathcal{H}}(u) = \langle V_{\mathcal{H}}(u), \delta\mathcal{F}(u) \rangle_0 $$

and $\{\mathcal{H}, \mathcal{F}\}(u) = \langle \mathcal{J}(u)\delta\mathcal{H}(u), \delta\mathcal{F}(u) \rangle_0$, it follows that $V_{\mathcal{H}}(u) = \mathcal{J}(u)\delta\mathcal{H}(u)$.
If the operator of the Poisson structure $\mathcal{J}(u) : X_{s+d_J} \to X_s$ is an isomorphism, the operator $\overline{\mathcal{J}}(u) := -\mathcal{J}^{-1}(u) : X_s \to X_{s+d_J}$ defines an skew-adjoint automorphism of order $-d_J$. Let us consider the closed 2-form $\omega_u := \overline{\mathcal{J}}(u)du \wedge du$ where $\overline{\mathcal{J}}(u) : X_d \to X_{-d}$ depends $C^1$-smoothly on $u \in O_d$ and defines a linear isomorphism $\overline{\mathcal{J}}(u) : X_d \to X_{d+d_J}$ for any $d_J \ge 0$. The 2-form $\omega$ supplies $O_d$ with a symplectic structure and $\overline{\mathcal{J}}(u)$ is called operator of the symplectic structure. The Hilbert scale $\{X_s\}_{s \in \mathbb{Z}}$ endowed with the 2-form $\omega$ forms a so-called *symplectic Hilbert scale*, denoted as $(\{X_s\}_{s, \omega})$. The Poisson bracket on the symplectic space $(O_d, \omega)$ is defined as

$$ \{\mathcal{F}, \mathcal{L}\}(u) = \omega_u(V_{\mathcal{F}}, V_{\mathcal{L}}) \quad \forall u \in O_d. $$

To a $C^1$-smooth function $\mathcal{H}$ on $O_d$ the symplectic structure $\omega$ corresponds the Hamiltonian vector field $V_{\mathcal{H}}$ defined as the continuous map $V_{\mathcal{H}} : O_d \to X_{-d-d_J}$ such that

$$ \forall u \in O_d \quad \omega_u(V_{\mathcal{H}}(u), \xi) = -d\mathcal{H}(u)\xi \quad \forall \xi \in T_u O_d \simeq X_d. $$

With this definition, the Hamiltonian vector field $V_{\mathcal{H}}$ satisfies $\langle \overline{\mathcal{J}}(u)V_{\mathcal{H}}(u), \xi \rangle_0 = -\langle \delta\mathcal{H}(u), \xi \rangle_0$. In particular, for any $C^1$-smooth function $\mathcal{H} : O_d \times \mathbb{R} \to \mathbb{R}$, we denote by $V_{\mathcal{H}}$ the non-autonomous vector field $V_{\mathcal{H}}(u, t) = \mathcal{J}(u)\delta\mathcal{H}(u, t)$. The corresponding Hamiltonian equation is

$$ \partial_t u(t) = V_{\mathcal{H}}(u, t) = \mathcal{J}(u)\delta\mathcal{H}(u, t). \qquad (2.2) $$

A partial differential equation, supplemented by appropriate boundary conditions, is called a Hamiltonian PDE if, under a suitable choice of a symplectic Hilbert scale $(\{X_s\}_{s, \omega})$, a domain $O_d \subset X_d$ and a Hamiltonian $\mathcal{H}$, it can be written in the form (2.2).

A vector field $V_{\mathcal{H}}$ on a manifold $O_d$ determines a phase flow, namely a one-parameter group of diffeomorphisms $\Phi_{V_{\mathcal{H}}}^t : O_d \to O_d$ satisfying $d_t \Phi_{V_{\mathcal{H}}}^t(u) = V_{\mathcal{H}}(\Phi_{V_{\mathcal{H}}}^t(u))$ for all $t \in \mathcal{T}$ and $u \in O_d$, with $\Phi_{V_{\mathcal{H}}}^0(u) = u$.

Hamiltonian dynamics is characterized by the existence of differential invariants, and symmetry-related conservation laws.

**Definition 2.4 (Invariants of motion).** A function $\mathcal{I} \in C^\infty(O_d)$ is an invariant of motion of the dynamical system with flow map $\Phi_{V_{\mathcal{H}}}^t$ if $\{\mathcal{I}, \mathcal{H}\}(u) = 0$ for all $u \in O_d$. Consequently, $\mathcal{I}$ is constant along the orbits of $V_{\mathcal{H}}$.

The Hamiltonian function, if time-independent, is an invariant of motion. A particular subset of the invariants of motion of a dynamical system is given by the *Casimir invariants*, i.e. functions on $O_d$ which $\{\cdot, \cdot\}$-commute with every other function on $O_d$.

**Definition 2.5 (Casimir invariants).** If $\mathfrak{g}$ is a Lie algebra with Lie product $\{\cdot, \cdot\}$, the *centralizer* of a subset $\mathcal{U}$ of $\mathfrak{g}$ is defined as $C_{\mathfrak{g}}(\mathcal{U}) := \{\mathcal{C} \in \mathfrak{g} : \{\mathcal{C}, \mathcal{F}\} = 0 \text{ for all } \mathcal{F} \in \mathcal{U}\}$. The centralizer $C_{\mathfrak{g}}(\mathfrak{g})$ of the Lie algebra itself is called the *centre* of $\mathfrak{g}$ and its elements are called Casimir functions.

## 2.2 Gradient flow dissipation

In this work we consider dissipation mechanisms in form of gradient flows. Since the seminal work [10], numerous contributions in the literature investigated the geometry and the behaviour
of dissipative partial differential equations. Far from being exhaustive, we refer to [22, 16, 3] for a detailed discussion of the topic.

Leveraging the existence theory for gradient flows in Banach spaces, the Rothe method [14] in particular, we formulate the following gradient flow. Let us consider a Hilbert scale {$X_s$}$_{s \in \mathbb{Z}}$. Note that each $X_s$ is reflexive and, since the injection $X_0 \hookrightarrow X_s$, $s \ge 0$, is continuous, we can identify the Gelfand triple $X_s \subset X_0 \subset X_{-s}$. Let us formulate a discrete in time gradient flow (whose limit is the continuous time gradient flow). Let $\mathcal{E}: X_s \to \mathbb{R}$ be a Fréchet differentiable convex functional. Let $\tau > 0$, and $v \in X_s$, we define

$$u = \arg \inf_{\tilde{u} \in X_s} \left( \frac{1}{2} \langle \tilde{u} - v, \tilde{u} - v \rangle_0 + \tau \mathcal{E}(\tilde{u}) \right).$$

The functional augmented by the norm of the Hilbert space is clearly convex and Fréchet differentiable. The Euler-Lagrange equations read as follows:

$$\langle u - v + \tau \delta \mathcal{E}(u), z \rangle_0 = 0 \quad \forall z \in X_s.$$

Since the above expression holds for any $z \in X_s$, we take $z = u - v$ and get

$$\|u - v\|_0^2 + \tau \langle \delta \mathcal{E}(u), u - v \rangle_0 = 0,$$

which gives

$$\langle \delta \mathcal{E}(u), v - u \rangle_0 = \frac{1}{\tau} \|u - v\|_0^2 \ge 0. \qquad (2.3)$$

Let us now exploit the convexity of the functional $\mathcal{E}$ to show that the energy is dissipated. The convexity of $\mathcal{E}$ gives, for $\vartheta \in [0, 1]$ and $u, v \in X_s$, that $\mathcal{E}((1 - \vartheta)u + \vartheta v) \le (1 - \vartheta)\mathcal{E}(u) + \vartheta\mathcal{E}(v)$. This entails $\vartheta^{-1}[\mathcal{E}(u + \vartheta(v - u)) - \mathcal{E}(u)] \le \mathcal{E}(v) - \mathcal{E}(u)$ for any $\vartheta \in [0, 1]$. In the limit $\vartheta \to 0$,

$$\langle \delta \mathcal{E}(u), v - u \rangle_0 \le \mathcal{E}(v) - \mathcal{E}(u).$$

Using (2.3) yields

$$0 \le \frac{1}{\tau} \|u - v\|_0^2 \le \mathcal{E}(v) - \mathcal{E}(u),$$

which clearly shows that $\mathcal{E}(u) \le \mathcal{E}(v)$. The continuous gradient flow is the limit, for $\tau \to 0$, of the time discrete gradient flow and it enjoys the same property of energy dissipation. Its weak formulation can be stated as follows: for all $v \in X_s$

$$\langle \partial_t u, v \rangle_0 + \langle \delta \mathcal{E}(u), v \rangle_0 = 0.$$

It is possible to show [14] that there exists a unique solution $u \in L^p(\mathcal{T}, V)$, $1 \le p \le \infty$, of this problem in the temporal interval $\mathcal{T}$.

In the present work we consider gradient flows with respect to a (potentially singular) metric. Let $\mathcal{S} = -\mathcal{E}$ be a concave functional and $\mathcal{G}(u)$ be a positive semi-definite self-adjoint operator. This gives rise to the following formulation:

$$\partial_t u = \mathcal{G}(u) \delta \mathcal{S}(u).$$
Let us consider the case in which $\mathcal{G}(u)$ has an inverse $\mathcal{G}(u)^{-1}$. In such a case, the gradient flow can be considered as the continuous limit of the following discrete in time gradient flow

$$u = \arg \inf_{\tilde{u} \in X_s} \left( \frac{1}{2} \langle \tilde{u} - v, \mathcal{G}(v)^{-1}(\tilde{u} - v) \rangle_0 - \tau \mathcal{S}(\tilde{u}) \right).$$

We can show that the quantity $\mathcal{S}$ is non-decreasing.

Let us consider the case in which $\mathcal{S}: X_0 \to \mathbb{R}$ is the entropy, and $\mathcal{G}(u): X_{-s} \to X_s$ is a linear self-adjoint operator. The equation for the time evolution of the entropy reads

$$d_t \mathcal{S}(u) = \langle \delta \mathcal{S}(u), \mathcal{G}(u) \delta \mathcal{S}(u) \rangle_0 \geq 0.$$

It is possible to introduce a positive semi-definite bilinear form $(\cdot, \cdot): X_s \times X_s \to \mathbb{R}$ defined as

$$(\mathcal{F}, \mathcal{L})(u) = \langle \mathcal{G}(u)\delta\mathcal{F}(u), \delta\mathcal{L}(u) \rangle_0 \quad \forall u \in X_s, \mathcal{F}, \mathcal{L}: X_s \to \mathbb{R}. \qquad (2.4)$$

We present some examples of dissipation terms, obtained as a gradient flow.

**Example 2.1 (Heat equation).** Let $\Omega \subset \mathbb{R}^n$ be a Lipschitz continuous domain. We consider the heat flow in the following setting: $X_0 = L^2(\Omega)$ and $u \in H_0^1(\Omega) = X_1$. By Rellich-Kondrakov theorem, the injection $X_1 \hookrightarrow X_0$ is compact. Let $\mathcal{G}(u) = -\Delta$ with domain $H_0^1(\Omega)$ and let

$$\mathcal{S}(u) = -\frac{1}{2} \int_{\Omega} u^2 dx.$$

Since its Fréchet derivative is $\delta\mathcal{S}(u) = -u$ one gets

$$d_t \mathcal{S}(u) = - \int_{\Omega} u \Delta u \, dx = \int_{\Omega} \nabla u \cdot \nabla u \, dx.$$

We obtain, as entropy evolution, the classical evolution of the $L^2$ norm of the solution of the heat equation.

The classical formulation of the heat equation as a gradient flow is the following. Let $X_0 = L^2(\Omega)$ and $X_1 = H_0^1(\Omega)$ as above. The entropy is the negative Dirichlet energy:

$$\mathcal{S}(u) = -\frac{1}{2} \int_{\Omega} |\nabla u|^2 dx.$$

Its Fréchet derivative is $\delta\mathcal{S}(u) = \Delta u \in X_{-1} = X_1'$. The operator $\mathcal{G}(u)$ is the identity, and for all $z \in H_0^1(\Omega)$:

$$\langle \Delta u, z \rangle_0 = - \langle \nabla u, \nabla z \rangle_0.$$

Remark that when we take $z = u$ we get: $\langle \delta\mathcal{S}(u), u \rangle_0 \leq 0$.

**Example 2.2 (Double bracket dissipation).** Another example of dissipation structure is provided by the so-called double bracket, see e.g. [4, 23]. In this formulation, the dissipation operator is given by $\mathcal{G} = -\mathcal{J}^2$ or $\mathcal{G} = -\overline{\mathcal{J}}^2$, where $\mathcal{J}$ and $\overline{\mathcal{J}}$ are the Poisson and symplectic operators defined
in Section 2.1. Let us show that the operator $G = -\mathcal{J}^2$ is positive definite and self-adjoint. Let $v \in X_{s+d_j}$, it holds:

$$ \langle G(u)v, v \rangle_0 = -\langle \mathcal{J}^2(u)v, v \rangle_0 = -\langle \mathcal{J}(u)v, \mathcal{J}^*(u)v \rangle_0 = \langle \mathcal{J}(u)v, \mathcal{J}(u)v \rangle_0 \geq 0. $$

The self-adjointness follows by observing that, for any $v, \xi \in X_{s+d_j}$,

$$ -\langle \mathcal{J}^2(u)v, \xi \rangle_0 = -\langle \mathcal{J}(u)v, \mathcal{J}^*(u)\xi \rangle_0 = \langle \mathcal{J}(u)v, \mathcal{J}(u)\xi \rangle_0 = -\langle \mathcal{J}^*(u)v, \mathcal{J}(u)\xi \rangle_0 = -\langle v, \mathcal{J}^2(u)\xi \rangle_0. $$

This is enough to show that the operator $\mathcal{J}^2$ induces an increase in the entropy functional and, therefore, it represents a candidate to describe dissipation mechanisms. Remark that, if we can consider, as the entropy, $S = -\mathcal{H}$, the double bracket dissipates the total energy (Hamiltonian) of the system.

## 2.3 Hamiltonian systems with gradient flow dissipation

The focus of this work is on Hamiltonian partial differential equations that are supplied by a gradient flow dissipation. The setting is as follows: let us consider the Gelfand triple $X_s \subset X_0 \subset X_{-s}$ and let $\mathcal{T} := (0, T] \subset \mathbb{R}$ be a bounded temporal interval. We consider the evolution equation (1.1), namely: find $u \in L^p(\mathcal{T}, X_s)$ such that

$$ \partial_t u = \mathcal{J}(u)\delta\mathcal{H}(u) + \mathcal{G}(u)\delta\mathcal{S}(u), \quad t \in \mathcal{T}, $$

supplied by the initial condition $u(0) = u_0 \in X_s$ and suitable boundary conditions that we assume to be encoded in the definition of the space $X_s$. Moreover, we assume that the functions $\mathcal{H}, \mathcal{S}: X_s \to \mathbb{R}$ satisfy $\mathcal{H}(u), \mathcal{S}(u) < +\infty$ for any $u \in X_s$. Further, the Poisson structure $\mathcal{J}$ and the operator $\mathcal{G}$ satisfy the conditions of Sections 2.1 and 2.2 and are such that $\mathcal{J}(u)\delta\mathcal{H} \in X_{-s}$ and $\mathcal{G}(u)\delta\mathcal{S} \in X_{-s}$ for any $u \in X_s$.

With a small abuse of notation, we will refer to the function $\mathcal{S}$ appearing in the generic evolution equation (1.1) as “entropy” even if the physical interpretation of such quantity might differ based on the particular example considered.

We present some examples of Hamiltonian systems with a dissipative term.

**Example 2.3 (Advection-diffusion equation).** Let us consider the advection-diffusion equation

$$ \partial_t u + \mathbf{v} \cdot \nabla u - \Delta u = 0 $$

in a Lipschitz continuous spatial domain $\Omega \subset \mathbb{R}^n$, with constant velocity $\mathbf{v}$ and homogeneous Dirichlet boundary conditions. Let $X_0 = L^2(\Omega)$ and $X_1 = H_0^1(\Omega)$. Let us define $\mathcal{H}: X_1 \to \mathbb{R}$ as

$$ \mathcal{H}(u) = \frac{1}{2} \int_{\Omega} |u|^2 dx. $$

Then, the Fréchet derivative of $\mathcal{H}$ is $\delta\mathcal{H}(u) = u$. We also introduce the operators $\mathcal{G}(u) = \mathcal{G} = \Delta$, with domain $X_1$ and $\mathcal{J}(u) = \mathcal{J} = -\mathbf{v} \cdot \nabla$. It holds $\mathcal{J}\delta\mathcal{H}(u) = -\mathbf{v} \cdot \nabla u \in X_0 \subset X_{-1}$, and $\mathcal{G}\delta\mathcal{H}(u) = \Delta u \in X_{-1}$. Thus, the advection-diffusion problem can then be cast as

$$ \partial_t u = (\mathcal{J} + \mathcal{G})\delta\mathcal{H}(u). $$
The kinetic energy of the system dissipates. Indeed,

$$d_t \mathcal{H}(u) = \int_{\Omega} u \Delta u \, dx = - \int_{\Omega} |\nabla u|^2 \, dx.$$

**Example 2.4** (Isentropic KdV equation). Let us consider the Korteweg-de Vries equation

$$\partial_t u + \alpha u \partial_x u + \eta \partial_x^3 u = \nu \partial_x^2 u \qquad (2.5)$$

in a one-dimensional spatial domain $\Omega \subset \mathbb{R}$ with periodic boundary conditions. We assume that $\alpha, \eta, \nu$ are positive bounded constants.

Let $X_0 = L^2(\Omega)$ and $X_1 = H^1(\Omega)$, we introduce the quantities

$$\mathcal{J}(u) = \mathcal{J} = -\partial_x, \qquad \mathcal{H}(u) = \int_{\Omega} \left( \frac{\alpha}{6} u^3 - \frac{\eta}{2} (\partial_x u)^2 \right) dx.$$

The Fréchet derivative of the Hamiltonian reads

$$\delta \mathcal{H}(u) = \frac{\alpha}{2} u^2 + \eta \partial_x^2 u.$$

For the KdV equation in Gardner's bracket formulation [7], a possible entropy is given by

$$S(u) = \int_{\Omega} u \, dx. \qquad (2.6)$$

This can be interpreted as the mass of the wave (especially when considering solitons). When having periodic boundary conditions, or a closed system, it is therefore natural to consider an evolution equation which preserves the mass, $d_t S = 0$, and dissipates an energy.

Let the operator $\mathcal{G}(u) : X_1 \to X_{-1}$ be defined as $\mathcal{G}(u)v = \nu v \partial_x^2 u$ for any $v \in X_1$. Then, the KdV equation (2.5) can be cast as in (1.1). In this setting the entropy is conserved since it holds

$$d_t S(u) = \langle \delta S(u), \partial_t u \rangle_0 = \langle 1, \mathcal{J} \delta \mathcal{H}(u) \rangle_0 + \nu \langle 1, \partial_x^2 u \rangle_0 = 0,$$

in view of the periodic boundary conditions and the fact that $1 \in \ker(\mathcal{J})$.

The kinetic energy of the system

$$\mathcal{E}(u) = \int_{\Omega} \frac{u^2}{2} \, dx$$

instead dissipates since

$$d_t \mathcal{E}(u) = \langle \delta \mathcal{E}(u), \partial_t u \rangle_0 = \langle u, -\alpha u \partial_x u - \eta \partial_x^3 u \rangle_0 + \nu \langle u, \partial_x^2 u \rangle_0 = -\nu \|\partial_x u\|_0^2 \le 0.$$

Although the Hamiltonian evolution does not have a definite sign, it is possible to show that what dissipates is the Hamiltonian of the so-called boosted soliton, that is the energy given by the Hamiltonian augmented by the kinetic energy of a reference frame moving at a certain velocity $\beta > 0$, that is

$$\hat{\mathcal{E}}(u) := \mathcal{H}(u) + \beta \int_{\Omega} \frac{u^2}{2} \, dx.$$
It is possible to show that there exists $\beta_0$ such that, for all $\beta > \beta_0$, we have $d_t \hat{\mathcal{E}} < 0$. Indeed,

$$d_t \hat{\mathcal{E}}(u) = \langle \delta \hat{\mathcal{E}}(u), \partial_t u \rangle_0 = \langle \delta \mathcal{H}(u), \partial_t u \rangle_0 + \beta \langle u, \partial_t u \rangle_0 = \nu \eta \|\partial_x^2 u\|_0^2 + \frac{\nu \alpha}{2} \langle \partial_x^2 u, u^2 \rangle_0 - \nu \beta \|\partial_x u\|_0^2.$$

Moreover, an equilibrium solution is given by

$$u_* = \varrho_0 = \frac{1}{|\Omega|} \int_{\Omega} u \, dx,$$

and the mass (2.6) is simply $S = \varrho_0|\Omega|$. Indeed, we have

$$\delta \mathcal{H}|_{u_*} = \frac{\alpha}{2} \varrho_0^2, \quad \text{and} \quad \partial_x^2 u_0 = 0,$$

and, hence,

$$\partial_t u_* = \frac{\alpha}{2} \mathcal{J} \varrho_0^2 = 0,$$

since the constant function is in the kernel of the operator $\mathcal{J}$.

### 2.3.1 The case of metriplectic or GENERIC systems

A special subclass of problems that can be described via (1.1) is the one where the dynamics is generated by both a Hamiltonian and an entropy function. The corresponding geometric formalism has been introduced under the name of metriplectic structure [17] in plasma physics and GENERIC formalism [9] in the context of non-equilibrium thermodynamics.

Metriplectic systems are characterized by an evolution equation as in (1.1), and further equipped with so-called *mutual degeneracy conditions*, namely

$$\mathcal{G}(u)\delta\mathcal{H}(u) = 0 \quad \text{and} \quad \mathcal{J}(u)\delta\mathcal{S}(u) = 0 \quad \forall u \in X_s. \qquad (2.7)$$

Note that, in view of the skew-adjointness of the operator $\mathcal{J}$ and the self-adjointness of $\mathcal{G}$, these conditions imply that $\{\mathcal{H}, \mathcal{S}\}(u) = (\mathcal{H}, \mathcal{S})(u) = 0$ for any $u \in X_s$. This, in turn, ensures that the Hamiltonian $\mathcal{H}$ is a conserved quantity of motion while the entropy $\mathcal{S}$ increases in time. Indeed,

$$d_t \mathcal{H}(u) = \langle \delta \mathcal{H}(u), \partial_t u \rangle_0 = \{\mathcal{H}, \mathcal{H}\}(u) + (\mathcal{H}, \mathcal{S})(u) = 0.$$

$$d_t \mathcal{S}(u) = \langle \delta \mathcal{S}(u), \partial_t u \rangle_0 = \{\mathcal{H}, \mathcal{S}\}(u) + (\mathcal{S}, \mathcal{S})(u) \geq 0.$$

## 3 The variational formulation of the problem

The variational formulation of the evolution problem (1.1) reads: Given $u_0 \in X_s$, find $u \in L^p(\mathcal{T}, X_s)$ such that

$$\langle \partial_t u, v \rangle_0 = \langle \mathcal{J}(u) \delta \mathcal{H}(u), v \rangle_0 + \langle \mathcal{G}(u) \delta \mathcal{S}(u), v \rangle_0, \quad \forall v \in X_s.$$

Since $\mathcal{J}(u)$ is skew-adjoint and $\mathcal{G}(u)$ is self-adjoint, we can rewrite the weak formulation as

$$\langle \partial_t u, v \rangle_0 = - \langle \delta \mathcal{H}(u), \mathcal{J}(u)v \rangle_0 + \langle \delta \mathcal{S}(u), \mathcal{G}(u)v \rangle_0, \quad \forall v \in X_s. \qquad (3.1)$$

Equation (3.1) is the starting point to derive a mixed variational formulation, detailed in the next section.
## 3.1 Mixed formulation

Let $s, d, e \in \mathbb{N} \setminus \{0\}$. Let us consider the Hilbert spaces $X_{s-d}$ and $X_{s-e}$ and remark that $X_s$ is dense in $X_{s-d}$ and in $X_{s-e}$. In the following, we are going to assume that $\delta\mathcal{H} \in X_{-s+d}$ and $\delta\mathcal{S} \in X_{-s+e}$.

The mixed formulation reads: Given $u_0 \in X_s$, find $(u, z, y) \in C^1(\mathcal{T}; X_s) \times C^0(\mathcal{T}; X_{-s+d}) \times C^0(\mathcal{T}; X_{-s+e})$ such that

$$ egin{cases} 
\langle \partial_t u, v \rangle_0 = -\langle z, \mathcal{J}(u)v \rangle_0 + \langle y, \mathcal{G}(u)v \rangle_0, & \forall v \in X_s, \\ 
\langle z - \delta\mathcal{H}, v \rangle_0 = 0, & \forall v \in X_s, \\ 
\langle y - \delta\mathcal{S}, v \rangle_0 = 0, & \forall v \in X_s. 
\end{cases} 
\quad (3.2) $$

### 3.1.1 Equivalence of the variational and mixed formulations

Let us consider the subspaces $K_{\mathcal{J}}$, $K_{\mathcal{G}} \subset X_s$ defined as

$$ K_{\mathcal{J}} := \{u \in X_s : \mathcal{J}(u)v \in X_{s-d} \text{ for all } v \in X_s\}, $$

$$ K_{\mathcal{G}} := \{u \in X_s : \mathcal{G}(u)v \in X_{s-e} \text{ for all } v \in X_s\}. $$

**Lemma 3.1.** Let $\xi \in X_{-s+d}$ and $u \in K_{\mathcal{J}}$. If $\langle \xi, v \rangle_0 = 0$ for all $v \in X_s$, then

$$ \langle \xi, \mathcal{J}(u)v \rangle_0 = 0, \quad \forall v \in X_s. $$

*Proof.* Let $\tilde{v} \in X_s$. Adding and subtracting this element to the second component of the inner product gives

$$ |\langle \xi, \mathcal{J}(u)v \rangle_0| = |\langle \xi, \tilde{v} \rangle_0 + \langle \xi, \mathcal{J}(u)v - \tilde{v} \rangle_0| $$

The first term of the right-hand side vanishes by hypothesis, since $\tilde{v} \in X_s$. The second term can be bounded, using Hölder's inequality, as

$$ |\langle \xi, \mathcal{J}(u)v - \tilde{v} \rangle_0| \leq \|\xi\|_{-s+a} \|\mathcal{J}(u)v - \tilde{v}\|_{s-d}. $$

Since $X_s$ is dense in $X_{s-d}$ we can choose $\tilde{v}$ in such a way that the second term becomes arbitrarily small, and the conclusion follows. $\square$

Remark that, by using the same density argument, we can prove that, given $\zeta \in X_{-s+e}$ and $u \in K_{\mathcal{G}}$, if $\langle \zeta, v \rangle_0 = 0$ for all $v \in X_s$ then $\langle \zeta, \mathcal{G}(u)v \rangle_0 = 0$.

**Proposition 3.2.** Assume that $\delta\mathcal{H} \in X_{-s+d}$ and $\delta\mathcal{S} \in X_{-s+e}$. If $t \in \mathcal{T} \mapsto u(t) \in X_s$ solution of the mixed formulation (3.2) belongs to $K_{\mathcal{J}} \cap K_{\mathcal{G}} \subset X_s$, then the mixed formulation is equivalent to the weak formulation (3.1) of the evolution equation (1.1).

*Proof.* By virtue of the result of Lemma 3.1, the solution $(u(t), z(t), y(t)) \in X_s \times X_{-s+d} \times X_{-s+e}$ of the mixed formulation (3.2) satisfies, for any $t \in \mathcal{T}$, $\langle z - \delta\mathcal{H}, \mathcal{J}(u)v \rangle_0 = 0$, and $\langle y - \delta\mathcal{S}, \mathcal{G}(u)v \rangle_0 = 0$ for any $v \in X_s$. From the first equation of (3.2), we thus get, for any $v \in X_s$,

$$ \langle \partial_t u, v \rangle_0 = - \langle z, \mathcal{J}(u)v \rangle_0 + \langle y, \mathcal{G}(u)v \rangle_0 = - \langle \delta\mathcal{H}(u), \mathcal{J}(u)v \rangle_0 + \langle \delta\mathcal{S}(u), \mathcal{G}(u)v \rangle_0. $$

Remark that, under the assumption that $X_s$ is dense in both $X_{-s+d}$ and $X_{-s+e}$, that is $2s \leq d, e$, then $\langle \xi, v \rangle_0 = 0$ for any $v \in X_s$ directly implies $\xi = 0 \in X_{-s+d}$. A similar result holds for $X_{-s+e}$.
# 4 The semi-discrete problem

## 4.1 Conformal variational discretisation

Under the assumption that $X_s$ is dense in both $X_{-s+d}$ and $X_{-s+e}$ ($2s \ge d, e$), we build a conformal variational discretisation of the mixed formulation (3.2) by introducing a finite dimensional subspace $V_N \subset X_s$, closed in $X_s$, and such that, for $N \to \infty$, $V_N$ is dense in $X_s$. Let $N \in \mathbb{N}^*$ and $\{v_i\}_{i=1}^N$ be a basis of $V_N$. We restrict the mixed formulation by taking test functions which belong to $V_N$. Moreover, by considering a Galerkin method, we look for numerical approximations $u_N, z_N, y_N$ in $V_N$. Let $a, b, c : \mathcal{T} \to \mathbb{R}^N$, we write the numerical approximation of the solution, for any $(t, x) \in \mathbb{R}^+ \times \Omega$, as:

$$u_N(t, x) = \sum_{i=1}^{N} a_i(t)v_i(x), \quad z_N(t, x) = \sum_{i=1}^{N} b_i(t)v_i(x), \quad y_N(t, x) = \sum_{i=1}^{N} c_i(t)v_i(x).$$

This leads to the following discretisation of the mixed formulation (3.2):

$$\left\{
\begin{array}{ll}
\displaystyle \sum_{j=1}^{N} \langle v_j, v_i \rangle_0 d_t a_j = - \sum_{j=1}^{N} (\langle v_j, \mathcal{J}(u_N)v_i \rangle_0 b_j + \langle v_j, \mathcal{G}(u_N)v_i \rangle_0 c_j), & 1 \le i \le N, \\
\\
\displaystyle \sum_{j=1}^{N} \langle v_j, v_i \rangle_0 b_j - \langle \delta\mathcal{H}(u_N), v_i \rangle_0 = 0, & 1 \le i \le N, \\
\\
\displaystyle \sum_{j=1}^{N} \langle v_j, v_i \rangle_0 c_j - \langle \delta\mathcal{S}(u_N), v_i \rangle_0 = 0, & 1 \le i \le N.
\end{array}
\right.
\quad (4.1)$$

In order to rewrite the discrete formulation of the problem, let us state the following result.

**Lemma 4.1.** Let $u_N \in V_N$ and let $a \in \mathbb{R}^N$ be the coefficients representing $u_N$ in the basis $\{v_i\}_{i=1}^N$ of $V_N$. Let $H, S: \mathbb{R}^N \to \mathbb{R}$ be defined as

$$H(a) = \mathcal{H}(u_N) = \mathcal{H}\left(\sum_{i=1}^{N} a_i v_i\right), \quad \text{and} \quad S(a) = \mathcal{S}(u_N) = \mathcal{S}\left(\sum_{i=1}^{N} a_i v_i\right).$$

Then, the Fréchet derivatives of $\mathcal{H}$ and $\mathcal{S}$, evaluated in $u_N \in V_N$, are such that

$$\partial_{a_i} H(a) = \langle \delta \mathcal{H}(u_N), v_i \rangle_0, \quad \partial_{a_i} S(a) = \langle \delta \mathcal{S}(u_N), v_i \rangle_0, \quad \forall 1 \le i \le N.$$

*Proof.* The proof is performed by direct computation starting from the definition of the Fréchet derivative. For all $v \in X_s$, it holds

$$\langle \delta \mathcal{H}(u), v \rangle_0 = \lim_{\vartheta \to 0} \frac{\mathcal{H}(u + \vartheta v) - \mathcal{H}(u)}{\vartheta}.$$

When $u = u_N \in V_N$, the variation is restricted to $v \in V_N$, $v = \sum_{j=1}^N q_j v_j$, and we can write

$$\langle \delta \mathcal{H}(u_N), v \rangle_0 = \lim_{\vartheta \to 0} \frac{\mathcal{H}(u_N + \vartheta v) - \mathcal{H}(u_N)}{\vartheta} = \lim_{\vartheta \to 0} \frac{\mathcal{H}\left(\sum_{j=1}^N (a_j + \vartheta q_j) v_j\right) - \mathcal{H}(u_N)}{\vartheta}.$$
Since {$q_j$}$_{1\le j\le N}$ is arbitrary, we can take $q_j = \delta_{ij}$ and get:

$$ \langle \delta \mathcal{H}(u_N), v_i \rangle_0 = \lim_{\vartheta \to 0} \frac{\mathcal{H}\left(\sum_{j=1}^N (a_j + \vartheta \delta_{ij})v_j\right) - \mathcal{H}\left(\sum_{j=1}^N a_j v_j\right)}{\vartheta} = \partial_{a_i} H(a). $$

The proof for the Fréchet derivative of $\mathcal{S}$ is analogous.

Let $u_N \in V_N$ and let $a \in \mathbb{R}^N$ be the coefficients of the expansion of $u_N$ in the basis {$v_i$}$_{i=1}^N$ of $V_N$. We introduce the following matrices: the mass matrix $M \in \mathbb{R}^{N \times N}$, and $J(a) \in \mathbb{R}^{N \times N}$ and $G(a) \in \mathbb{R}^{N \times N}$ which are the representation of the symplectic and metric operators, respectively. Their components read, for any $1 \le i, j \le N$,

$$ M_{i,j} = \langle v_j, v_i \rangle_0, \quad J(a)_{i,j} = \langle \mathcal{J}(u_N)v_j, v_i \rangle_0, \quad G(a)_{i,j} = \langle \mathcal{G}(u_N)v_j, v_i \rangle_0. \quad (4.2) $$

With this notation and the results of Lemma 4.1, the semi-discrete (in space) mixed formulation of the problem can be written, for any $t \in \mathcal{T}$, as:

$$ \left\{ \begin{array}{l} M d_t a(t) = J(a(t)) b(t) + G(a(t)) c(t) \\ M b(t) - \nabla H(a(t)) = 0, \\ M c(t) - \nabla S(a(t)) = 0. \end{array} \right. \qquad (4.3) $$

## 4.2 Geometric structure of the semi-discrete problem

In the next Lemma we prove that the discretisation of the continuous Poisson bracket from Definition 2.3 yields a bilinear form that is still a Poisson bracket.

**Lemma 4.2.** Let

$$ \hat{J}(a) := M^{-1}J(a)M^{-1} \quad \forall a \in \mathbb{R}^N. \qquad (4.4) $$

Let the bilinear form {$\cdot, \cdot$}$_N: C^1(\mathbb{R}^N) \times C^1(\mathbb{R}^N) \to C^1(\mathbb{R}^N)$ be defined as:

$$ \{F, L\}_N(a) = \hat{J}(a)\nabla F(a) \cdot \nabla L(a), \quad \forall a \in \mathbb{R}^N, $$

where $\cdot$ denotes the Euclidean product in $\mathbb{R}^N$ and $F, L: \mathbb{R}^N \to \mathbb{R}$ are real valued functions in $C^1(\mathbb{R}^N)$. Let $\mathcal{F}, \mathcal{L}: X_s \to \mathbb{R}$ be $C^1$-functionals, such that, for all elements $u_N \in X_s$ with $u_N = \sum_{i=1}^N a_i v_i$, it holds $\mathcal{F}(u_N) = F(a)$ and $\mathcal{L}(u_N) = L(a)$. Then,

1. {$F, L$}$_N$ is the restriction to $V_N \subset X_s$ of the continuous Poisson bracket {$\mathcal{F}, \mathcal{L}$}.

2. {$\cdot, \cdot$}$_N$ is a Poisson bracket, namely it is skew-symmetric, it satisfies the Leibniz rule and the Jacobi identity.

*Proof.* Let $\mathcal{F}, \mathcal{L}: X_s \to \mathbb{R}$ and $u = u_N \in V_N$. Then $\delta\mathcal{F}(u_N) = \sum_{j=1}^N f_j v_j$ where $f_j$, by virtue of Lemma 4.1, is the jth component of the vector $f = M^{-1}\nabla F(a)$. Similar considerations can be
made for $\delta\mathcal{L}$. Using Definition 2.3 of the Poisson bracket, it holds

$$ \begin{aligned} \{\mathcal{F}, \mathcal{L}\}(u_N) &= \langle \mathcal{J}(u_N)\delta\mathcal{F}(u_N), \delta\mathcal{L}(u_N) \rangle_0 = \sum_{i,j=1}^{N} \langle \mathcal{J}(u_N)f_j v_j, \ell_i v_i \rangle_0 \\ &= \sum_{i,j,k,m=1}^{N} M_{ik}^{-1} M_{jm}^{-1} \langle \mathcal{J}(u_N)v_j, v_i \rangle_0 \partial_{a_k} L(a) \partial_{a_m} F(a) \\ &= \sum_{k,m=1}^{N} (\nabla L)_k(a) \left( \sum_{i,j=1}^{N} M_{k,i}^{-1} J(a)_{i,j} M_{j,m}^{-1} \right) (\nabla F)_m(a) = \{F, L\}_N(a). \end{aligned} $$

It is easy to show that $\{\cdot, \cdot\}_N$ is skew-symmetric and satisfies the Leibniz rule. The second item follows from the fact that the continuous bracket satisfies the Jacobi identity. $\square$

Observe that an analogous result holds for the bilinear form associated with $\mathcal{G}$. Indeed, the bilinear form (2.4) restricted to $V_N$ gives

$$ (\mathcal{F}, \mathcal{L})(u_N) = \langle \mathcal{G}(u_N)\delta\mathcal{F}(u_N), \delta\mathcal{L}(u_N) \rangle_0 = \sum_{k,m=1}^{N} (\nabla L)_k(a) \left( \sum_{i,j=1}^{N} M_{k,i}^{-1} G(a)_{i,j} M_{j,m}^{-1} \right) (\nabla F)_m(a). $$

The discrete positive semi-definite bilinear form associated with (2.4) is then

$$ (F, L)_N(a) = \hat{G}(a)\nabla F(a) \cdot \nabla L(a), \quad \hat{G}(a) := M^{-1}G(a)M^{-1}, \quad \forall a \in \mathbb{R}^N. $$

We are now in the position of proving the main results of the present work. In particular, we are going to show that the semi-discrete in space mixed formulation (4.3) is a finite dimensional Hamiltonian system with a dissipation mechanism in the form of a finite dimensional gradient flow. We are going to separate the proof in two propositions, treating the cases $\mathcal{J}(u) = 0$ and $\mathcal{G}(u) = 0$ separately.

**Proposition 4.3.** Let $\mathcal{G}(u) = 0$. Then, the mixed formulation (4.3):

1. It is a consistent semi-discretised form of the continuous Hamiltonian system obtained by setting $\mathcal{G}(u) = 0$ in (1.1).

2. It is Hamiltonian, with the discrete Poisson bracket defined in Lemma 4.2.

*Proof.* Since $\mathcal{G}(u) = 0$ we can restrict to the first and second equations of the mixed formulation (4.3) and take $c = 0$. From the second equation we have $b = M^{-1}\nabla H(a)$, that we can substitute into the first equation and get

$$ d_t a = M^{-1}J(a)M^{-1}\nabla H(a) = \hat{J}(a)\nabla H(a), $$

and this concludes the first part of the proof.

Moreover, Lemma 4.2 ensures that $\hat{J}(a)$ is a discrete Poisson tensor. This is enough to conclude that the system is a finite dimensional Hamiltonian system. $\square$
Let us now prove that the discretisation of the gradient flow leads to a finite dimensional gradient flow.

**Proposition 4.4.** Let $\mathcal{J}(u) = 0$. Then, the mixed formulation (4.3):

1. It is a consistent semi-discretised form of the continuous gradient flow system obtained by setting $\mathcal{J}(u) = 0$ in (1.1).

2. It is a finite dimensional gradient flow for the entropy $S$.

*Proof.* Since $\mathcal{J}(u) = 0$ we can simply set $b = 0$ and consider the first and third equations in the mixed formulation (4.3). Replacing $c = M^{-1}\nabla S(a)$ into the mixed formulation yields

$$d_t a = M^{-1}G(a)M^{-1}\nabla S(a) = \hat{G}(a)\nabla S(a).$$

Let us remark that $G(a)$ is a symmetric positive semi-definite matrix by construction and $M$ is symmetric and positive definite. Then, the discrete dissipation tensor $\hat{G}(a)$ is symmetric positive semi-definite. This implies that the discrete system entropy $S(a)$ cannot decrease:

$$d_t S(a) = \nabla S^\top(a) \hat{G}(a) \nabla S(a) \geq 0,$$

and this concludes the proof.

□

## 4.3 Conservation properties of the semi-discrete problem

Another interesting question concerns the quantities which are conserved at discrete level by the mixed discretisation we have described. Concerning the Hamiltonian part, a major role is played by the Casimirs of the system, as introduced in Definition 2.5.

**Proposition 4.5.** Let $u \in X_s$ and $\mathcal{G}(u) = 0$. The following holds:

1. If, for any $u \in X_s$, $\ker(\mathcal{J}(u))$ has finite dimension $r \leq N$, and $\ker(\mathcal{J}(u)) \subset V_N$, then, for any $a \in \mathbb{R}^N$, $\hat{\mathcal{J}}(a)$ has a kernel of dimension $r$, and the discretisations of the continuous Casimirs are discrete Casimirs.

2. If, for any $u \in X_s$, $\ker(\mathcal{J}(u)) \perp_{X_s} V_N$, then the restriction of any Casimir to $V_N$ is a constant function.

*Proof.* Let $\mathcal{F}: X_s \to \mathbb{R}$ be a Casimir of $\mathcal{J}$ that is

$$\{\mathcal{F}, \mathcal{C}\}(u) = \langle \mathcal{J}(u)\delta\mathcal{F}(u), \delta\mathcal{C}(u) \rangle_0 = 0 \quad \forall u \in X_s, \mathcal{C}: X_s \to \mathbb{R}.$$

If $\ker(\mathcal{J}(u)) \subset V_N$ for any $u \in X_s$, then the Fréchet derivative $\delta\mathcal{F}$ of the Casimir $\mathcal{F}$ belongs to $V'_N$. This means that we can write

$$\delta\mathcal{F}(u_N) = \sum_{j=1}^{N} f_j v_j \quad \forall u_N = \sum_{j=1}^{N} a_j v_j,$$

where $f_j = f_j(a) = (M^{-1}\nabla F(a))_j$, and $F(a) = \mathcal{F}(u_N)$.
Moreover, since $\mathcal{F}$ is a Casimir, it holds

$$0 = \langle \mathcal{J}(u_N)\delta\mathcal{F}(u_N), v \rangle_0 = \sum_{j=1}^{N} f_j(a) \langle \mathcal{J}(u_N)v_j, v \rangle_0, \quad \forall u_N \in V_N, v \in X_s.$$

Taking $v = v_i$, we get $(J(a)f(a))_i = 0$. This implies that the vector $f(a) = M^{-1}\nabla F(a)$ belongs to the kernel of $J(a)$, for any $a \in \mathbb{R}^N$, and, hence, $\nabla F(a)$ belongs to the kernel of $\hat{\mathcal{J}}(a) = M^{-1}J(a)M^{-1}$. This concludes the first part of the proof.

Concerning the second item, let us proceed by contradiction. Let us assume that the discrete Casimir $F$ depends on $a$ so that $\nabla F(a)$ is, in general, different from zero. It follows that there exists a function $\delta\mathcal{F} \in V_N$ such that its coefficients in the basis $\{v_i\}_{i=1}^N$ are given by $f_i = \sum_{l=1}^N M_{il}^{-1}\partial_{a_l}F(a)$. This function is the Fréchet derivative of a functional $\mathcal{F}$ restricted to $V_N$. For $\mathcal{F}$ to be a Casimir, it must hold: $\delta\mathcal{F} \in \ker(\mathcal{J}(u))$. This could be possible only if at least part of the kernel is in $V_N$, hence contradicting the hypothesis. $\square$

**Proposition 4.6 (Conservation of invariants).** Let $K \subseteq X_s$ be the phase space of problem (1.1). Assume that $K$ is convex and contains the origin. Let $\mathcal{F}: K \to \mathbb{R}$ be an invariant of motion of problem (1.1) as in Definition 2.4. Let $F: \mathbb{R}^N \to \mathbb{R}$ be a real valued function belonging to $C^1(\mathbb{R}^N)$ defined as $F(a) = \mathcal{F}(u_N)$ for any $u_N = \sum_{i=1}^N a_i v_i \in V_N \subset K$. Then, $F$ is an invariant of motion of the semi-discrete problem (4.3).

*Proof.* Let $u_N \in V_N$. Since $\mathcal{F}$ is an invariant of (3.2) it holds $\{\mathcal{F}, \mathcal{H}\}(u_N) + (\mathcal{F}, \mathcal{S})(u_N) = 0$. This is equivalent to $\{F, H\}_N(a) + (F, S)_N(a) = 0$. Since $u_N$ is chosen arbitrarily, $F$ is an invariant of motion of the semi-discrete problem (4.3). $\square$

Next result shows that the equilibria of (1.1) belonging to the approximation space $V_N$ are equilibria of the semi-discrete problem.

**Proposition 4.7.** Let $u^* \in X_s$ be an equilibrium point of problem (1.1) that is

$$\mathcal{J}(u^*)\delta\mathcal{H}(u^*) + \mathcal{G}(u^*)\delta\mathcal{S}(u^*) = 0.$$

If $u^*$ belongs to $V_N$, then it is equilibrium point of the semi-discrete system (4.1).

*Proof.* Since $u^* \in V_N$ it can be expressed as $u^* = \sum_{i=1}^N a_i^* v_i$. Moreover, since $u^*$ is equilibrium point it holds $\langle f(u^*), v \rangle_0 = 0$ for any $v \in X_s$, where $f(u) := \mathcal{J}(u)\delta\mathcal{H}(u) + \mathcal{G}(u)\delta\mathcal{S}(u)$ for any $u \in X_s$. Taking $v = v_i \in V_N$ yields

$$0 = \langle f(u^*), v_i \rangle_0 = \sum_{j=1}^{N} b_j \langle \mathcal{J}(u^*)v_j, v_i \rangle_0 + \sum_{j=1}^{N} c_j \langle \mathcal{G}(u^*)v_j, v_i \rangle_0 \\ = \sum_{j=1}^{N} b_j J(a^*)_{i,j} + \sum_{j=1}^{N} c_j G(a^*)_{i,j} = (J(a^*)b + G(a^*)c)_i.$$

Since $b = M^{-1}\nabla H(a^*)$ and $c = M^{-1}\nabla S(a^*)$, one has $J(a^*)M^{-1}\nabla H(a^*) + G(a^*)M^{-1}\nabla S(a^*) = 0$ and the conclusion follows. $\square$
### 4.3.1 The case of metriplectic or GENERIC systems

The mixed variational formulation (3.2) holds also for the metriplectic systems introduced in Section 2.3.1. The mutual degeneracy conditions are imposed weakly as

$$ \langle \mathcal{G}(u)\delta\mathcal{H}(u), w \rangle_0 = 0 \quad \text{and} \quad \langle \mathcal{J}(u)\delta\mathcal{S}(u), w \rangle_0 = 0 \quad \forall u, w \in X_s. \quad (4.5) $$

The semi-discretisation of problem (3.2) in $V_N \subset X_s$ can be performed as in (4.3) for metriplectic systems as well. Moreover, the mutual degeneracy conditions (4.5) are discretised by imposing $w = v_i \in V_N$ and taking $u_N \in V_N$; this yields

$$ 0 = \langle \mathcal{G}(u_N)\delta\mathcal{H}(u_N), v_i \rangle_0 = \sum_{j,k=1}^{N} M_{jk}^{-1} \partial_{a_k} H(a) \langle \mathcal{G}(u_N)v_j, v_i \rangle_0 = (G(a)M^{-1}\nabla H(a))_i, \quad (4.6) $$

and analogously $J(a)M^{-1}\nabla S(a) = 0$ for all $a \in \mathbb{R}^N$. The resulting semi-discrete system remains metriplectic, as shown in the next result.

**Proposition 4.8.** The semi-discrete in space mixed formulation (4.3) together with the conditions (4.6) is a finite dimensional metriplectic system.

*Proof.* Following the reasoning in the proofs of Propositions 4.3 and 4.4, the semi-discrete system reads

$$ d_t a(t) = \hat{J}(a(t))\nabla H(a(t)) + \hat{G}(a(t))\nabla S(a(t)). $$

Moreover, the discretisation of the mutual degeneracy conditions (4.6) yields

$$ \hat{G}(a)\nabla S(a) = 0 \quad \text{and} \quad \hat{J}(a)\nabla S(a) = 0. \quad (4.7) $$

□

As a consequence of the previous result, the Hamiltonian $H$ is a conserved quantity and $S$ is increasing. Indeed, from (4.7), it holds

$$ \begin{aligned} (H, S)_N(a) &= \nabla H(a)^\top \hat{G}(a) \nabla S(a) = 0, \\ \{H, S\}_N(a) &= \nabla H(a)^\top \hat{J}(a) \nabla S(a) = 0. \end{aligned} $$

The results follows from the skew-symmetry of $\hat{J}(a)$ and the positive semi-definiteness of $\hat{G}(a)$, for any $a \in \mathbb{R}^N$.

## 4.4 Convergence results

Considering the convergence properties of the proposed semi-discretisation, we first notice that, in all the cases in which the mixed formulation is equivalent to the classical variational formulation, it enjoys the same convergence properties. In this section, we propose a proof of convergence in norm $\|\cdot\|_0$ which makes use of Lipschitz continuity hypotheses of the Fréchet derivatives and of the operators $\mathcal{J}$ and $\mathcal{G}$. In particular, let us consider the following assumptions:
(H1) The Fréchet derivatives $\delta\mathcal{H}$ and $\delta\mathcal{S}$ are Lipschitz continuous: there exist $c_H, c_S > 0$ such that

$$ \|\delta\mathcal{H}(u) - \delta\mathcal{H}(v)\|_{-s+d} \le c_H \|u - v\|_s, \quad \text{and} \quad \|\delta\mathcal{S}(u) - \delta\mathcal{S}(v)\|_{-s+e} \le c_S \|u - v\|_s. $$

for any $u, v \in X_s$.

(H2) The map $u \in X_s \mapsto \mathcal{J}(u) \in \mathcal{L}(X_{-s+d}, X_{-s})$ is Lipschitz continuous: there exists a constant $c_J > 0$ such that, for any $u, v \in X_s$,

$$ \|\mathcal{J}(u) - \mathcal{J}(v)\|_{\mathcal{L}(X_{-s+d}, X_{-s})} \le c_J \|u - v\|_s. $$

Analogously, concerning the dissipation, we assume that the map $u \in X_s \mapsto \mathcal{G}(u) \in \mathcal{L}(X_{-s+e}, X_{-s})$ is Lipschitz continuous: there exists a constant $c_G > 0$ such that, for any $u, v \in X_s$,

$$ \|\mathcal{G}(u) - \mathcal{G}(v)\|_{\mathcal{L}(X_{-s+e}, X_{-s})} \le c_G \|u - v\|_s. $$

**Proposition 4.9 (Convergence of the semi-discretisation scheme).** Let $s, d, e \in \mathbb{N} \setminus \{0\}$, $2s \ge d, e$. Under the hypotheses (H1) and (H2), the solution $u_N(t) \in V_N$ of the semi-discrete system (4.3) converges to the solution of the continuous mixed formulation (3.2): For any $t \in \mathcal{T}$,

$$ \lim_{N \to \infty} \|u(t) - u_N(t)\|_0 = 0. $$

*Proof.* In what follows, for the sake of compactness of the notation, we omit the time dependence on $u(t), z(t), y(t)$ and write simply $u, z, y$.

Let $P_N, Q_N$, and $R_N$ be three residual functions, orthogonal to $V_N$, that is, $\langle P_N, v_N \rangle_0 = 0$, $\langle Q_N, v_N \rangle_0 = 0$, and $\langle R_N, v_N \rangle_0 = 0$, for any $v_N \in V_N$, obtained by testing the numerical approximation of the mixed problem (4.1) against a generic element of the space, $v \in X_s$:

$$ \left\{ \begin{array}{l} \langle \partial_t u_N, v \rangle_0 = \langle \mathcal{J}(u_N) z_N + \mathcal{G}(u_N) y_N + Q_N, v \rangle_0, \\ \langle z_N - \delta\mathcal{H}(u_N) - P_N, v \rangle_0 = 0, \\ \langle y_N - \delta\mathcal{S}(u_N) - R_N, v \rangle_0 = 0. \end{array} \right. \qquad (4.8) $$

Note that, since, for $N \to \infty$, $V_N$ is dense in $X_s$, it holds $\|Q_N\|_s \to 0$ as $N \to \infty$ and similarly for $R_N$ and $P_N$ in the associated norms. Let the error be the function $e = u - u_N$ where $u$ is solution of the mixed formulation (3.2) and $u_N$ is solution of (4.8). We derive the mixed formulation for the error by taking the difference of the mixed formulation for the solution $u$ and the one for its numerical approximation $u_N$: for any $v \in X_s$,

$$ \left\{ \begin{array}{l} \langle \partial_t e, v \rangle_0 = \langle \mathcal{J}(u)z + \mathcal{G}(u)y - \mathcal{J}(u_N)z_N - \mathcal{G}(u_N)y_N - Q_N, v \rangle_0, \\ \langle z - z_N - \delta\mathcal{H}(u) + \delta\mathcal{H}(u_N) + P_N, v \rangle_0 = 0, \\ \langle y - y_N - \delta\mathcal{S}(u) + \delta\mathcal{S}(u_N) + R_N, v \rangle_0 = 0. \end{array} \right. $$

As a particular test function, we choose $v = e$ and get

$$ \left\{ \begin{array}{l} \langle \partial_t e, e \rangle_0 = \langle \mathcal{J}(u)z + \mathcal{G}(u)y - \mathcal{J}(u_N)z_N - \mathcal{G}(u_N)y_N - Q_N, e \rangle_0, \\ \langle z - z_N - \delta\mathcal{H}(u) + \delta\mathcal{H}(u_N) + P_N, e \rangle_0 = 0, \\ \langle y - y_N - \delta\mathcal{S}(u) + \delta\mathcal{S}(u_N) + R_N, e \rangle_0 = 0. \end{array} \right. $$
From the first equation, by making use of Hölder's inequality, we can write:

$$ \frac{1}{2} \partial_t \|e\|_0^2 \le \|\mathcal{J}(u)z + \mathcal{G}(u)y - \mathcal{J}(u_N)z_N - \mathcal{G}(u_N)y_N\|_{-s} \|e\|_s + \|Q_N\|_{-s} \|e\|_s. $$

We focus on the first term on the right-hand side. By means of the triangular inequality we can separate the Hamiltonian and dissipative parts,

$$ \|\mathcal{J}(u)z + \mathcal{G}(u)y - \mathcal{J}(u_N)z_N - \mathcal{G}(u_N)y_N\|_{-s} \le \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N\|_{-s} + \|\mathcal{G}(u)y - \mathcal{G}(u_N)y_N\|_{-s}. $$

We can treat the terms separately. Concerning the first term, we add and subtract $\mathcal{J}(u)z_N$ and make use of the triangular inequality:

$$ \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N\|_{-s} \le \|\mathcal{J}(u)(z - z_N)\|_{-s} + \|\mathcal{J}(u) - \mathcal{J}(u_N)\|_{-s} \|z_N\|_{-s}. $$

Using the Lipschitz continuity of $\mathcal{J}(u)$ with respect to $u$ (hypothesis $(H_2)$), we write:

$$ \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N\|_{-s} \le \hat{c}_J(u) \|z - z_N\|_{-s+d} + c_J \|e\|_s \|z_N\|_{-s+d}, \quad (4.9) $$

where $\hat{c}_J(u) := \|\mathcal{J}(u)\|_{\mathcal{L}(X_{-s+d}, X_{-s})}$. In order to bound $\|z - z_N\|_{-s+d}$ we consider the equation for the element of the cotangent bundle

$$ \langle z - z_N, v \rangle_0 = \langle \delta\mathcal{H}(u) - \delta\mathcal{H}(u_N) - P_N, v \rangle_0 \quad \forall v \in X_s. $$

Let us consider the norm in $X_{-s+d}$ as given by

$$ \|z\|_{-s+d} = \sup_{z_* \in X_{-s-d} \setminus \{0\}} \frac{|\langle z, z_* \rangle_0|}{\|z_*\|_{s-d}}. $$

Let $v_* \in X_{-s-d}$, $\|v_*\|_{-s-d} = 1$, be supremiser of the duality product $\langle z - z_N, v \rangle_0$. Since $X_s \subset X_{-s-d}$ with a dense inclusion, we consider the sequence in $X_s$ that converges to $v^*$. We get:

$$ \|z - z_N\|_{-s+d} = |\langle \delta\mathcal{H}(u) - \delta\mathcal{H}(u_N) - P_N, v_* \rangle_0| \le \|\delta\mathcal{H}(u) - \delta\mathcal{H}(u_N)\|_{-s+d} + \|P_N\|_{-s+d}. $$

The Hamiltonian is Fréchet differentiable, and, by virtue of the hypothesis $(H_1)$ we have Lipschitz continuity of the Fréchet derivative. This makes it possible to introduce a constant $c_H > 0$ such that:

$$ \|z - z_N\|_{-s+d} \le c_H \|e\|_s + \|P_N\|_{-s+d}. $$

We use this result into (4.9) and get:

$$ \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N\|_{-s} \le \hat{c}_J(u)(c_H \|e\|_s + \|P_N\|_{-s+d}) + c_J \|e\|_s \|z_N\|_{-s+d}, $$

which can be rewritten as

$$ \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N\|_{-s} \le (c_H \hat{c}_J(u) + c_J \|z_N\|_{-s+d}) \|e\|_s + \hat{c}_J(u) \|P_N\|_{-s+d}. $$

Concerning the dissipation mechanism, we can proceed in the same way, and get

$$ \|\mathcal{G}(u)y - \mathcal{G}(u_N)y_N\|_{-s} \le (c_S \hat{c}_G(u) + c_G \|y_N\|_{-s+e}) \|e\|_s + \hat{c}_G(u) \|R_N\|_{-s+e}, $$
where $\hat{c}_G(u) := \|G(u)\|_{L(X_{-s+e}, X_{-s})}$. Summing these two terms, we get that there exists $C_N(u) > 0$ such that

$$ \|\mathcal{J}(u)z - \mathcal{J}(u_N)z_N + \mathcal{G}(u)y - \mathcal{G}(u_N)y_N\|_{-s} \le C_N(u)\|e\|_s + \hat{c}_J(u)\|P_N\|_{-s+d} + \hat{c}_G(u)\|R_N\|_{-s+e}, $$

where $C_N(u) := c_H \hat{c}_J(u) + c_J \|z_N\|_{-s+d} + c_S \hat{c}_G(u) + c_G \|y_N\|_{-s+e}$. Henceforth, it holds

$$ \frac{1}{2} \partial_t \|e\|_0^2 \le C_N(u) \|e\|_s^2 + (\|Q_N\|_{-s} + \|\mathcal{J}(u)\|_{L(X_{-s+d}, X_{-s})} \|P_N\|_{-s+d} + \|\mathcal{G}(u)\|_{L(X_{-s+d}, X_{-s})} \|R_N\|_{-s+e}) \|e\|_s. $$

The equation is homogeneous in $\|e\|_s$ and since we can approximate the initial condition of the system, $e_0 = u(0) - u_N(0)$ is such that $\lim_{N \to \infty} \|e_0\|_s = 0$, we can conclude that the mixed semi-discrete formulation converges to the continuous mixed formulation of the system. Furthermore, the error with respect to the solution of the problem tends to zero in the norm of $X_0$. $\square$

# 5 Temporal discretisation of the semi-discrete problem

Concerning the numerical time integration of the semi-discrete problem (4.3) we consider two well-established methods: the Average Vector Field (AVF) introduced in [20] and the implicit midpoint rule.

First, we consider a partition of the temporal interval $\mathcal{T} = (0, T]$ into $N_t$ subintervals $(t^k, t^{k+1}]$, $0 \le k \le N_t - 1$. For simplicity of exposition, we assume this partition to be uniform so that $t^k = t^0 + k\Delta t$, for $0 \le k \le N_t$ with $\Delta t = N_t^{-1}T$. The extension to the non-uniform case is trivial. Given the ordinary differential equation

$$ d_t a = f(a), \quad a(0) = a^0, \qquad (5.1) $$

the AVF scheme approximates the solution at time $t^{k+1}$, $k \ge 0$, as

$$ \frac{a^{k+1} - a^k}{\Delta t} = \int_0^1 f((1-\xi)a^k + \xi a^{k+1}) d\xi. \qquad (5.2) $$

In our case, the velocity field of the flow is given by

$$ f(a) = \hat{J}(a)\nabla H(a) + \hat{G}(a)\nabla S(a). $$

**Lemma 5.1.** Assume that $J$ and $G$ in (4.3) are constant-valued. Then, the AVF scheme applied to (4.3)

* exactly preserves the Hamiltonian, when $G = 0$;

* satisfies the dissipative law, when $J = 0$;

* is a metriplectic time integrator.
Proof. First note that, if $J$ and $G$ are constant-valued then so are $\hat{J}$ and $\hat{G}$. Let $y(\xi) := (1-\xi)a^k + \xi a^{k+1}$ and let $\cdot$ denote the Euclidean inner product. Multiplying Equation (5.2) by the integral of $\nabla H(y(\xi))$ in $[0, 1]$, one gets

$$ \frac{1}{\Delta t} \int_{0}^{1} (a^{k+1} - a^k) \cdot \nabla H(y(\xi)) d\xi = \int_{0}^{1} \nabla H(y(\xi)) d\xi \hat{J} \int_{0}^{1} \nabla H(y(\xi)) d\xi + \int_{0}^{1} \nabla H(y(\xi)) d\xi \hat{G} \int_{0}^{1} \nabla S(y(\xi)) d\xi. $$

If $G=0$ or $G\nabla H(y) = 0$, for any $y \in \mathbb{R}^N$, then the above equality gives

$$ \frac{1}{\Delta t} (H(a^{k+1}) - H(a^k)) = 0 $$

owing to the skew-symmetry of $\hat{J}$. A similar reasoning can be carried out for the dissipation of $S$, leading to the conclusion. $\square$

Note that, the AVF scheme, in its classic version, is a second order temporal integrator. It is, however, not a symplectic integrator.

The implicit midpoint rule applied to the ODE (5.1) is given by

$$ \frac{a^{k+1} - a^k}{\Delta t} = f(a^{k+1/2}), \qquad (5.3) $$

where the midpoint state $a^{k+1/2}$ is defined as $a^{k+1/2} = \frac{a^{k+1} + a^k}{2}$ for any $0 \le k \le N_t$.

**Lemma 5.2.** Assume that $H$ is a quadratic function, i.e., $H(a) = a^\top La$ with $L$ symmetric positive definite. Then, the implicit midpoint rule applied to problem (4.3) ensures exact conservation of the Hamiltonian whenever $G=0$ or the system is metriplectic.

*Proof.* By the linearity of $\nabla H$ one has, for any $k \ge 0$,

$$ \begin{aligned} \frac{1}{\Delta t} (H(a^{k+1}) - H(a^k)) &= \frac{1}{\Delta t} \nabla H(a^{k+1/2}) \cdot (a^{k+1} - a^k) = \nabla H(a^{k+1/2}) \cdot f(a^{k+1/2}) \\ &= \nabla H(a^{k+1/2}) \cdot \hat{J}(a^{k+1/2}) \nabla H(a^{k+1/2}) + \nabla H(a^{k+1/2}) \cdot \hat{G}(a^{k+1/2}) \nabla S(a^{k+1/2}) \\ &= \nabla H(a^{k+1/2}) \cdot \hat{G}(a^{k+1/2}) \nabla S(a^{k+1/2}). \end{aligned} $$

The last term vanishes when either $G=0$ or the null condition $G\nabla H(y) = 0$ holds for any $y \in \mathbb{R}^N$. $\square$

An analogous reasoning holds for the dissipation of $S$.

**Lemma 5.3.** Assume that $S$ is a quadratic function, i.e., $S(a) = a^\top La$ with $L$ symmetric positive definite. Then, the implicit midpoint rule applied to problem (4.3) gives dissipation of $S$ whenever $J=0$ or the system is metriplectic.

The result of the previous Lemmas implies that, if $H$ and $S$ are quadratic functions of the state, then the implicit midpoint rule is a metriplectic integrator. Note that, for linear equations, the AVF scheme and the implicit midpoint rule coincides.
# 6 Numerical tests

In this section we present several numerical tests to illustrate the properties of the method. In all numerical tests we take $X_0 = L^2(\Omega)$ with the $L^2$ inner product denoted as $\langle \cdot, \cdot \rangle_0$.

The errors will be measured in relative $L^2$ norm, and we define:

$$e(t) = \frac{\|u(t) - u_N(t)\|_{L^2(\Omega)}}{\|u(t)\|_{L^2(\Omega)}}.$$

We will take the $L^\infty(\mathcal{T})$ norm in time:

$$\|e(t)\|_{L^\infty(\mathcal{T})} = \sup_{t \in \mathcal{T}} e(t),$$

and, in the following, we will refer to it as $L^\infty(\mathcal{T}, L^2(\Omega))$ relative error.

For the spatial approximation of the problem we consider a finite element partition $\Omega_N$ of the domain $\Omega$ in the sense of Ciarlet. In the numerical tests of the next sections $V_N$ will be the space of piecewise linear finite elements on the partition $\Omega_N$ and $N$ will correspond to the number of mesh nodes.

## 6.1 Korteveg-de Vries (KdV) equation

The KdV equation is an example of non-linear dispersive wave equation, which is an infinite dimensional integrable Hamiltonian system [7]. Let $\Omega = [0, 20\pi]$, we will consider here its standard form, with a dissipation in a form of a heat term, as presented in Eq. (2.5),

$$\partial_t u + \alpha u \partial_x u + \eta \partial_x^3 u = \nu \partial_x^2 u,$$

with periodic boundary conditions. In what follows, we consider $\alpha = 6$, $\eta = 1$, and we let $\nu$ vary. First we consider the non-dissipative case, which corresponds to the classical KdV equation, and then the dissipative one. Let us recall that the Hamiltonian is:

$$\mathcal{H}(u) = \int_{\Omega} \frac{\alpha}{6} u(x)^3 - \frac{\eta}{2} (\partial_x u(x))^2 dx,$$

and $\mathcal{J}(u) = \mathcal{J} = -\partial_x$, which does not depend on the solution. We consider the entropy

$$\mathcal{S}(u) = -\frac{1}{2} \int_{\Omega} u(x)^2 dx,$$

with $\mathcal{G}(u) = \mathcal{G} = -\nu \partial_x^2$.

In order to choose the semi-discretisation in space, let us write the mixed formulation by assuming the functions regular enough and check the minimal regularity required for the formulation to hold. Let $v \in H_{\text{per}}^1(\Omega)$ where the subscript refers to the periodic boundary conditions. We have

$$\left\{
\begin{array}{l}
\langle \partial_t u, v \rangle_0 = -\langle \partial_x z, v \rangle_0 + \nu \langle \partial_x y, \partial_x v \rangle_0, \\
\langle z, v \rangle_0 = \langle \frac{\alpha}{2} u^2, v \rangle_0 - \eta \langle \partial_x u, \partial_x v \rangle_0, \\
\langle y, v \rangle_0 = -\langle u, v \rangle_0.
\end{array}
\right.$$
Let us remark that, since the entropy is minus the energy of the solution $\delta S = -u$, the variable $y$ is somehow redundant. However, we keep it here in order to study the mixed formulation of the problem. We observe that choosing $V_N$ such that $V_N \subset H_{\text{per}}^1(\Omega)$ provides enough regularity to make the formulation wellposed. We use piecewise linear finite elements $V_N = \mathbb{P}_{\text{per}}^1(\Omega_N)$ on a uniform partition $\Omega_N$ with $N$ nodes.

The mass

$$M(u) = \int_{\Omega} u(x) \, dx$$

is a Casimir of the Hamiltonian part of the system; indeed $\delta M = 1 \in \ker(\mathcal{J})$. The dissipative term, given the periodic boundary conditions, does not affect the mass. We can therefore conclude that the mass is a conserved quantity of the system.

Given the results of Proposition 4.6, since $1 \in \mathbb{P}_{\text{per}}^1(\Omega_N)$, we can infer that the proposed semi-discretisation conserves the mass.

Let $M, K \in \mathbb{R}^{N \times N}$ be the mass and stiffness matrices respectively, whose components are defined as:

$$M_{i,j} := \langle v_j, v_i \rangle_0 \quad K_{i,j} := \langle \partial_x v_j, \partial_x v_i \rangle_0.$$

Using the notation from Section 4.1, the semi-discretisation of the system reads as follows

$$\left\{ \begin{array}{l} M d_t a = -Ab - \nu K a \\ M b = \frac{\alpha}{2} R a \otimes a - \eta K a \end{array} \right. \qquad (6.1)$$

where $A_{i,j} := \langle \partial_x v_j, v_i \rangle_0$ and $R \in \mathbb{R}^{N \times N^2}$ is the triple tensor defined as

$$R_{k,h} = \langle v_j v_k, v_i \rangle_0 \quad h = (i-1)N + j, \quad \forall 1 \le i, j \le N.$$

For the numerical time integration of the semi-discrete problem (6.1) we consider the implicit midpoint rule and the AVF method described in Section 5. In both discretisations the linear terms are evaluated at the midpoint state, while the nonlinear terms are discretised differently according to the chosen integrator. In particular, the fully discrete system corresponding to (6.1) reads

$$\left\{ \begin{array}{l} M \frac{a^{k+1} - a^k}{\Delta t} = -Ab^* - \nu K a^{k+1/2}, \\ Mb^* = \alpha R(d_2^{-1} a^{k+1} \otimes a^{k+1} + d_2^{-1} a^k \otimes a^k + d_1^{-1} a^{k+1} \otimes a^k) - \eta K a^{k+1/2}. \end{array} \right.$$

It can be easily verified that, for the implicit midpoint rule, $b^* = b^{k+1/2}$ which corresponds to $d_1 = 4, d_2 = 8$; whereas for the AVF, $d_1 = d_2 = 6$. We consider, as unknown, the vectors $a^{k+1}, b^*$, which are the discretisation of the elements of the manifold and its cotangent bundle, respectively. Moving the unknowns to the left hand side, we get

$$\left\{ \begin{array}{l} \left( M + \frac{\Delta t}{2} \nu K \right) a^{k+1} + \Delta t A b^* = \left( M - \frac{\Delta t}{2} \nu K \right) a^k \\ \left( \frac{\eta}{2} K - d_1^{-1} \alpha R a^k \otimes \right) a^{k+1} - d_2^{-1} \alpha R a^{k+1} \otimes a^{k+1} + Mb^* = \left( -\frac{\eta}{2} K + d_2^{-1} \alpha R a^k \otimes \right) a^k. \end{array} \right.$$

The nonlinearity in the latter equation is treated by means of a fixed point iteration.

In the numerical tests of this subsection we set the number $N$ of degrees of freedom of the spatial discretisation to $N = 512$ and consider $N_t = 800$ equispaced time iterations in the temporal interval $T = [0, 15]$.
Classical, non dissipative solutions

In this section, we consider $\nu = 0$, and get the classical KdV equation. We are in the periodic boundary conditions setting, but since the domain is large, the one soliton solution on the real line is close to the true solution of the system, and we will use it as a benchmark.

The analytical solution is:

$$u(x, t) = \mathrm{sech}^2\left(\frac{\sqrt{2}}{2}(x - 5\pi - 2t)\right), \qquad (6.2)$$

which corresponds to a solitary wave propagating at constant speed.

In Figure 1, on the left, we show a zoom of the solution, at final time, compared to the analytical one. We can see that the proposed method yields a numerical solution which does not exhibit spurious oscillations, and has a correct amplitude. The error is mainly due to a delay with respect to the analytical solution. On the same figure, on the right, we show the plot of the convergence of the solution towards the analytical one, in $L^\infty(\mathcal{T}, L^2(\Omega))$ relative error. In this figure, we compare the errors obtained with both the time integration schemes. We clearly see that they are equivalent, the errors being slightly smaller for the AVF method. We can distinguish a pre-asymptothic phase, in which the decay of the error is slower than $N^{-2}$, the expected rate, a phase in which the schemes are second-order in space, and, lastly, for larger space resolution, we can observe that the convergence slows down. This is due to the interaction between the space discretisation and the time integration errors.

<figure><img src="image_2.png" /><figcaption>Figure 1: On the left: zoom on the final time numerical solution compared to the analytical one. On the right, mesh convergence in log-log scale.</figcaption></figure>

In Figure 2 we report the conservation of mass (on the left) and of the Hamiltonian (on the right). The reference values are computed using the analytical solution (6.2). The proposed spatial discretisation ensures that the semi-discrete system preserves the mass and the Hamiltonian. The error that it is observed in the numerical tests is purely due to the time integration scheme. In terms of mass conservation, both the midpoint integration method and the AVF conserve the mass
up to the tolerance of the fixed-point iteration. The difference between the two time schemes can be seen on the Hamiltonian conservation, which is roughly 4 orders of magnitudes better for AVF, which conserves the Hamiltonian up to the tolerance of the nonlinear iteration. This is expected since the Hamiltonian for the KdV equation is not a quadratic function and, thus, the implicit midpoint scheme does not preserve it exactly.

<figure><img src="image_3.png" /><figcaption>Figure 2: Mass (on the left) and Hamiltonian (on the right) conservation in time, for the solution of the KdV equation with no dissipation.</figcaption></figure>

## Dissipative solutions

In this section, we assess the behaviour of the method when introducing dissipation, namely by setting $\nu = 1/4$.

In Figure 3, on the left, we show the time evolution of the numerical solution, obtained with the AVF time integration method. We plot the solution in space every 100 iteration, and highlight, in red, the initial condition. We can see that the dissipation mechanism induces a decrease in time of the amplitude of the soliton. Since the proposed approximation guarantees the separation between conservative and dissipative processes, the Hamiltonian part is the classical KdV equation, and hence a decrease in amplitude is related to a decrease in the soliton speed. This phenomenon can be seen visually when considering the soliton peak at multiple times. The mass of the solution is conserved, and for the sake of brevity we do not report the plot of the mass conservation, which is analogous to the one shown for the classical KdV solution in Figure 2.
<figure><img src="image_4.png" /><figcaption>Figure 3: On the left, time evolution of the solution of the dissipative KdV equation; the red curve is the initial condition evaluated at the mesh nodes. On the right, the discrete residual in time of the entropy production rate.</figcaption></figure>

It is interesting to consider, in this case, the evolution of the entropy $S(u) = -\frac{1}{2} \int_{\Omega} u(x)^2 dx$. We can see that, at continuous level, it holds:

$$d_t S(u) = \nu \int_{\Omega} \partial_x u(x)^2 dx.$$

When considering $\nu > 0$, we do not have an analytical solution for which we can exactly compute the entropy production rate. However, we can introduce the discrete residual of the entropy production

$$\mathcal{R}(u_N^{k+1}) := \frac{S(u_N^{k+1}) - S(u_N^k)}{\Delta t} - \nu \int_{\Omega} \left( \partial_x \left( \frac{u_N^{k+1} + u_N^k}{2} \right) \right)^2 dx.$$

The evolution of this quantity in time is shown, for both time integration schemes, in Figure 3, on the right. We can see that both methods have a good agreement, at discrete level, in terms of entropy production.

## 6.2 The 2D Navier-Stokes equations in vorticity formulation

We recall the formulation of the 2D Navier-Stokes equations as an infinite dimensional Hamiltonian flow, the incompressible Euler equations, to which we add a dissipation mechanism.

Let $\Omega \subset \mathbb{R}^3$ be a two-dimensional spatial domain. Let $\omega : \mathcal{T} \times \Omega \to \mathbb{R}$ be the vorticity defined as $\omega = \text{curl } u = \text{div}(\mathcal{J}u)$ where $u$ is the velocity field of the flow and

$$\mathcal{J} := \begin{pmatrix} 0 & \text{Id} \\ -\text{Id} & 0 \end{pmatrix}.$$

Introducing the stream function $\Psi : \mathcal{T} \times \Omega \to \mathbb{R}$ defined as $\text{curl} \Psi = -\mathcal{J} \nabla \Psi = u$, then the 2D
incompressible Navier-Stokes equations read

$$ \begin{cases} \partial_t \omega + \{\Psi, \omega\} = \nu \Delta \omega \\ \omega = \Delta \Psi \end{cases} $$

where $\nu$ is a positive constant, the stream function $\Psi$ is related to the vorticity via the Laplace-Beltrami operator $\Delta$, and $\{\cdot, \cdot\}$ is the Poisson bracket defined as $\{\Psi, \omega\} = \nabla\Psi \cdot \mathcal{J}\nabla\omega$.

The system can be written as a Hamiltonian system with dissipation; thereby

$$ \partial_t \omega = \mathcal{J}(\omega)\delta\mathcal{H}(\omega) + \nu\mathcal{G}\delta\mathcal{E}(\omega), \qquad (6.3) $$

where the function $\mathcal{H}$ is the kinetic energy

$$ \mathcal{H}(\omega) = -\frac{1}{2} \int_{\Omega} \Psi \omega \, dx, \quad \delta\mathcal{H} = -\Psi $$

and $\mathcal{E}$ is the entrophy

$$ \mathcal{E}(\omega) = \frac{1}{2} \int_{\Omega} |\omega|^2 \, dx, \quad \delta\mathcal{E} = \omega. $$

The metric tensor $\mathcal{G} = \Delta$ is the Laplacian, while the Poisson tensor is defined as $\mathcal{J}(\omega)z = \{z, \omega\} = \text{ad}_z^*\omega = -\nabla\omega \cdot \text{curl} z = -\nabla\omega \cdot \mathcal{J}\nabla z$. Observe that the entrophy $\mathcal{E}$ satisfies

$$ \mathcal{J}(\omega)\delta\mathcal{E}(\omega) = \mathcal{J}(\omega)\omega = \{\omega, \omega\} = 0, $$

and it is, thus, a Casimir of the $\mathcal{J}$ operator. Moreover, it dissipates since

$$ d_t \mathcal{E}(\omega) = \langle \delta\mathcal{E}(\omega), \mathcal{J}(\omega)\delta\mathcal{H}(\omega) + \nu\Delta\delta\mathcal{E}(\omega) \rangle_0 = -\nu \int_{\Omega} |\nabla\omega|^2 \, dx = -2\nu\mathcal{P}(\omega(t)), $$

where $\mathcal{P}$ is called palinstrophy. The kinetic energy $\mathcal{H}$ satisfies

$$ \mathcal{G}\delta\mathcal{H}(\omega) + \omega = -\Delta\Psi + \omega = 0 $$

and it dissipates since

$$ d_t \mathcal{H}(\omega) = \langle \delta\mathcal{H}(\omega), \nu\Delta\delta\mathcal{E}(\omega) \rangle_0 = -\nu \int_{\Omega} |\omega|^2 \, dx = -\nu 2\mathcal{E}(\omega(t)). $$

Let $X_s = X_1 = H^1(\Omega)$ and $X_{s-d} = X_{s-e} = X_1$. Following Section 3.1, the mixed formulation of problem (6.3) reads: for any $v \in X_1$,

$$ \begin{cases} \langle \partial_t \omega, v \rangle_0 = -\langle z, \mathcal{J}(\omega)v \rangle_0 - \nu \langle \nabla y, \nabla v \rangle_0, \\ \langle z, v \rangle_0 = -\langle \Psi, v \rangle_0, \\ \langle y, v \rangle_0 = \langle \omega, v \rangle_0, \end{cases} $$

where the stream function is solution of the Laplace equation $\langle \omega, v \rangle_0 = -\langle \nabla\Psi, \nabla v \rangle_0$ for any $v \in X_1$. From this formulation, we can check that the choice $V_N = \mathbb{P}^1(\Omega_N)$ is a simple pertinent choice to
discretise the system. We will use classical Lagrange piecewise linear finite element to discretise both the vorticity and the stream function.

Using the notation from Section 4.1, the semi-discrete system reads

$$ \begin{cases} M d_t a = J(a)b - \nu K a \\ K b = M a \end{cases} \qquad (6.4) $$

where $M, K \in \mathbb{R}^{N \times N}$ are the mass and stiffness matrix, respectively, while

$$ J(a)_{i,j} = \sum_{k=1}^{N} a_k \langle \{v_k, v_j\}, v_i \rangle_0 = \sum_{k=1}^{N} a_k \langle \nabla v_k \cdot J \nabla v_j, v_i \rangle_0 \quad 1 \le i, j \le N. $$

We verify that (6.4) has the Hamiltonian dissipative structure and preserves the conservation and dissipation laws at the semi-discrete level. First, observe that the discrete Hamiltonian and discrete entropy are

$$ \begin{aligned} H(a) &= -\frac{1}{2} \langle \omega_N, \Psi_N \rangle_0 = \frac{1}{2} a^\top M b = \frac{1}{2} a^\top M K^{-1} M a, & \nabla H(a) &= M K^{-1} M a, \\ E(a) &= \frac{1}{2} a^\top M a, & \nabla E(a) &= M a. \end{aligned} \qquad (6.5) $$

With these definitions, equation (6.4) can be written in the form

$$ d_t a = \hat{J}(a) \nabla H(a) + \nu \hat{G} \nabla E(a), $$

where $H$ and $E$ are as in (6.5), $\hat{G} = -M^{-1}KM^{-1}$ is symmetric negative definite, $\hat{J}(a)$ is defined as the matrix-valued function $\hat{J}: a \in \mathbb{R}^N \mapsto \hat{J}(a) \in \mathbb{R}^{N \times N}$ such that $\hat{J}(a) = M^{-1}J(a)M^{-1}$ or componentwise

$$ \hat{J}(a)_{i,j} := \sum_{k,l,s=1}^{N} M_{i,s}^{-1} M_{\ell,j}^{-1} a_k \langle \{v_k, v_\ell\}, v_s \rangle_0 \quad 1 \le i, j \le N. $$

The matrix $\hat{J}(a)$ is skew-symmetric for any $a \in \mathbb{R}^N$ and it satisfies the Jacobi identity, see Lemma 4.2.

We check the null condition, namely $\hat{J}(a)\nabla E(a) = 0$ for any $a \in \mathbb{R}^N$,

$$ \begin{aligned} (\hat{J}(a)Ma)_i &= \sum_{k,s=1}^{N} M_{i,s}^{-1} a_k J(a)_{s,k} = \sum_{k,l,s=1}^{N} M_{i,s}^{-1} a_k a_l \langle \{v_l, v_k\}, v_s \rangle_0 \\ &= \sum_{s=1}^{N} M_{i,s}^{-1} \langle \{\omega_N, \omega_N\}, v_s \rangle_0 = 0. \end{aligned} $$

The other condition is

$$ \hat{G} \nabla H(a) + a = -M^{-1} K M^{-1} (M K^{-1} M a) + a = 0. $$
For the numerical time integration of the semi-discrete problem (6.4), we consider the implicit midpoint rule and the AVF method described in Section 5; thereby

$$ \left(M + \frac{\Delta t}{2} \nu K\right) a^{k+1} - d_1^{-1} \Delta t J(a^{k+1}) b^{k+1} - d_2^{-1} \Delta t J(a^{k+1}) b^k - d_2^{-1} \Delta t J(a^k) b^{k+1} \\ = \left(M - \frac{\Delta t}{2} \nu K\right) a^k + d_1^{-1} \Delta t J(a^k) b^k. $$

The implicit midpoint rule corresponds to $d_1 = d_2 = 4$; whereas for the AVF scheme to $d_1 = 3$, $d_2 = 6$.

Since the Hamiltonian (6.5) of the Navier-Stokes problem is a quadratic function, in the absence of dissipation, it is conserved by the implicit midpoint rule as shown in Lemma 5.2.

## The Navier-Stokes equations on the torus

When using periodic boundary conditions, it is possible to have analytical solutions of the Navier-Stokes equations, similar to the classical Taylor vortex solutions. We refer to [26] for further details.

Let $\Omega = [0, 2\pi]^2$ and $V = H_{\text{per}}^1(\Omega)$ and let the temporal interval be $T = (0, 1]$. We take $\nu \in \mathbb{R}^+$ and $\lambda = 25$. We consider the stream function given by

$$ \Psi(t, x, y) = e^{-\nu\lambda t} \left[ \frac{1}{4} \cos(3x) \sin(4y) - \frac{1}{5} \cos(5y) - \frac{1}{5} \sin(5x) \right], \quad (t, x, y) \in \mathcal{T} \times \Omega. $$

We can verify that the vorticity $\omega$ satisfies

$$ \omega = -\lambda\Psi. $$

The evolution is such that the velocity field is always orthogonal to the gradient of the vorticity, and hence we can immediately verify that the above-stated solution verifies the Navier-Stokes equation in vorticity formulation.

In the present case, we consider several values of the diffusivity, $\nu \in \{10^{-2}, 10^{-4}, 0\}$. The last value corresponds to the steady solution of the incompressible Euler flow, for which $u \cdot \nabla\omega = 0$, and $\partial_t\omega = 0$. This limit case is challenging as there are no dissipation mechanisms acting. We consider $N_t = 200$ equispaced time steps.
<figure><img src="image_5.png" /><figcaption>Figure 4: Initial condition for the stream function $\Psi$ (on the left), and convergence in mesh (on the right), for the Navier-Stokes equation solution on the torus, when $\nu = 10^{-2}$.</figcaption></figure>

In Figure 4, on the left, we show the initial condition for the stream function $\Psi_0$. As we can see, it is periodic, but it does not have the symmetries of the classical Taylor vortex solution. On the right, the convergence in mesh is shown. In particular, we plot the error on the vorticity in $L^\infty(T, L^2(\Omega))$ when $\nu = 10^{-2}$. After a pre-asymptotic phase, we can see that we recover the expected rate $h^2 = N^{-1}$. The results are analogous when taking $\nu = 10^{-4}$.

<figure><img src="image_6.png" /><figcaption>Figure 5: Hamiltonian decay (left), and error in semi-logarithmic scale (right).</figcaption></figure>

In Figure 5, on the left, we show the decay of the Hamiltonian of the flow when $\nu = 10^{-2}$, and $N = 128^2$. We can see that both the AVF scheme and the implicit midpoint rule have the same behaviour, and follow the analytical trend. The error we have on the analytical value of the Hamiltonian depends solely on the fact that the initial condition does not belong to the space
$V_N$, and hence we have a discretisation error. On the right of Figure 5, we show the evolution of the error on the Hamiltonian decay, that is the difference between the values of the Hamiltonian $\mathcal{H}(\omega_N)$ obtained by the numerical scheme and the theoretical decay that the Hamiltonian should have, namely $\mathcal{H}(\omega_N(0))e^{-2\lambda t}$. The plot shows that both schemes have a good agreement (although not to machine precision) with respect to the theoretical decay.

<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the error in the enstrophy (left), and in the palinstrophy (right).</figcaption></figure>

In Figure 6 we show the same plot for the error in the decay of enstrophy and palinstrophy. For both these quantities, whose analytical decay in time is the same as for the Hamiltonian, we can see the same behaviour, for both time integration schemes.

Let us investigate the behaviour of the method when taking the limit case $\nu = 0$. In Figure 7, on the left, we show the evolution of the relative error in the vorticity approximation, as function of time, when taking $N = 128^2$. We can see that the maximum in time of the relative error is similar to what we got in the dissipative case. On the left, we plot the error with respect to the analytical evolution of the Hamiltonian, which is constant in time since there is no dissipation. We can see that this error is of the order of the machine precision, for both the time integration schemes. The method is able to correctly represent the fact that the Hamiltonian is a constant of motion. The same happens for the enstrophy. The error on the solution is mainly related to the fact that the initial condition, which is an equilibrium point for the infinite dimensional Hamiltonian system does not belong to $V_N$. The projection on $V_N$ of the initial condition is not an equilibrium point for the discrete Hamiltonian system, and this is why the solution, and hence the error, evolves in time.

In conclusion, the behaviour of the semi-discretisation in space is rather uniform with respect to the level of dissipation which we consider, including the case in which we do not have any dissipation.
<figure><img src="image_8.png" /><figcaption>Figure 7: On the left, evolution of the relative error in the approximation of the vorticity. On the right, error in the Hamiltonian evolution.</figcaption></figure>

## The Navier-Stokes equations on the sphere

In this section, we consider the Navier-Stokes equations on the sphere. The equations are the same, with the only difference that the gradient is replaced by the covariant gradient and the Laplace operator by the Laplace-Beltrami operator on the sphere.

Let the coordinates on the sphere be the angles $(\vartheta, \varphi) \in [0, \pi] \times [0, 2\pi]$. As an example of analytical solution, we consider the stream function at initial time:

$$ \Psi_0(\vartheta, \varphi) = \frac{1}{2} \sin(\vartheta) \cos(\varphi). $$

This function is related to spherical harmonics, and we can prove that we obtain a solution of the form:

$$ \Psi(t, \vartheta, \varphi) = \Psi_0 e^{-2\nu t}, $$

and $\omega(t, \vartheta, \varphi) = -2\Psi(t, \vartheta, \varphi)$.
<figure><img src="image_9.png" /><figcaption>Figure 8: On the left, initial condition for the stream function. On the right, convergence in mesh, for the two time integration schemes, when varying the space resolution.</figcaption></figure>

The numerical simulations are performed on a time interval $T = [0, 1]$, by taking $\Delta t = 10^{-2}$. The dissipation coefficient is taken, as before $\nu \in \{10^{-2}, 10^{-4}, 0\}$. Both the implicit midpoint rule and the AVF time integration methods are tested.

In Figure 8, on the left, we show the function $\Psi_0$ on the sphere. On the right, we report the plot of the relative $L^\infty(T, L^2(\Omega))$ error as function of the number of degrees of freedom used to discretise the mesh. As we can see, we obtain the same behaviour as for the planar case considered in the previous section. There are no substantial differences concerning the time integration scheme. The errors plotted in the figure are the ones obtained for $\nu = 10^{-2}$. The behaviour of the error obtained for $\nu = 10^{-4}$ and in the inviscid case are analogous and are not reported for the sake of brevity.

As a second test case, we are going to consider a scenario in which we do not have an analytical solution.

<figure><img src="image_10.png" /><figcaption>Figure 9: Evolution in time of the vorticity</figcaption></figure>
In particular, we initialise the vorticity as a sum of 512 point vortices, half of which with positive sign, the other half with negative sign, and intensity $\bar{\omega}_i = 400$.

In Figure 9 we show the solution obtained on the finest mesh $N \approx 1.2 \cdot 10^4$, when $\nu = 10^{-2}$, at $t = \{0.1, 0.3, 0.5\}$ using the midpoint rule as time integration scheme. In this figure we can see the dissipation of the enstrophy as well as the vortices coalescence phenomenon which is typical of two dimensional flows.

<figure><img src="image_11.png" /><figcaption>Figure 10: Evolution of the Hamiltonian (left), of the enstrophy (centre), and of the palinstrophy (right).</figcaption></figure>

In Figure 10 we plot the time evolution of the Hamiltonian (on the left), enstrophy (at the centre) and palinstrophy (on the right) for this simulation. As we can see the three quantities decay monotonically, as expected.

# 7 Conclusions and perspectives

We have addressed the problem of finding a geometric approximation of nonconservative evolution equations characterised by a Hamiltonian part and a gradient flow dissipation. We have proposed a conformal variational discretisation based on the mixed formulation of the problem that ensures preservation of the geometric structure of the problem and convergence of the approximate solution to the exact one. The focus of this work is on the spatial discretisation of problem (1.1), the investigation of numerical time integration scheme is left for future work.

To achieve sufficiently accurate solution the method might require a relatively large approximation space. The introduction of a dynamical adaptive discretisation in space could considerably speed up numerical simulations. An adaptive discretisation that does not hinder the geometric structure of the problem will be subject of future investigation.

# References

[1] R. Abraham and J. E. Marsden. *Foundations of mechanics*. Second edition. Addison-Wesley Publishing Company, Inc., Redwood City, CA., 1987.

[2] B. D. Andrews and P. E. Farrell. *High-order conservative and accurately dissipative numerical integrators via auxiliary variables*. 2024. arXiv: 2407.11904.
[3] S. Arnrich, A. Mielke, M. A. Peletier, G. Savaré, and M. Veneroni. "Passing to the limit in a Wasserstein gradient flow: from diffusion to reaction". In: *Calculus of Variations and Partial Differential Equations* 44.3 (2012), pp. 419–454.

[4] A. Bloch, P. Krishnaprasad, J. E. Marsden, and T. S. Ratiu. "The Euler-Poincaré equations and double bracket dissipation". In: *Communications in mathematical physics* 175.1 (1996), pp. 1–42.

[5] H. Egger. "Structure preserving approximation of dissipative evolution problems". In: *Numer. Math.* 143.1 (2019), pp. 85–106.

[6] C. Eldred and F. Gay-Balmaz. "Single and double generator bracket formulations of multicomponent fluids with irreversible processes". In: *Journal of Physics A: Mathematical and Theoretical* 53.39 (Aug. 2020), p. 395701.

[7] C. S. Gardner. "Korteweg-de Vries equation and generalizations. IV. The Korteweg-de Vries equation as a Hamiltonian system". In: *Journal of Mathematical Physics* 12.8 (Aug. 1971), pp. 1548–1551.

[8] G. Golo, V. Talasila, A. van der Schaft, and B. Maschke. "Hamiltonian discretization of boundary control systems". In: *Automatica* 40.5 (2004), pp. 757–771.

[9] M. Grmela and H. C. Öttinger. "Dynamics and thermodynamics of complex fluids. I. Development of a general formalism". In: *Phys. Rev. E* 56 (6 1997), pp. 6620–6632.

[10] R. Jordan, D. Kinderlehrer, and F. Otto. "The variational formulation of the Fokker-Planck equation". In: *SIAM journal on mathematical analysis* 29.1 (1998), pp. 1–17.

[11] A. N. Kaufman. "Dissipative Hamiltonian systems: A unifying principle". In: *Physics Letters A* 100.8 (1984), pp. 419–422.

[12] S. B. Kuksin. *Analysis of Hamiltonian PDEs*. Vol. 19. Oxford Lecture Series in Mathematics and its Applications. Oxford University Press, Oxford, 2000, pp. xii+212.

[13] B. Leimkuhler and S. Reich. "Hamiltonian PDEs". In: *Simulating Hamiltonian Dynamics*. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2005, pp. 316–356.

[14] J.-L. Lions. *Quelques méthodes de résolution des problèmes aux limites non linéaires*. Dunod, Paris; Gauthier-Villars, Paris, 1969, pp. xx+554.

[15] A. Mielke. *An introduction to the analysis of gradients systems*. 2023. arXiv: 2306.05026.

[16] A. Mielke and U. Stefanelli. "Weighted energy-dissipation functionals for gradient flows". In: *ESAIM: Control, Optimisation and Calculus of Variations* 17.1 (2011), pp. 52–85.

[17] P. J. Morrison. "A paradigm for joined Hamiltonian and dissipative systems". In: *Physica D: Nonlinear Phenomena* 18.1 (1986), pp. 410–419.

[18] H. C. Öttinger. "Relativistic and nonrelativistic description of fluids with anisotropic heat conduction". In: *Phys. A* 254.3 (1998), pp. 433–450.

[19] H. C. Öttinger and M. Grmela. "Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism". In: *Phys. Rev. E* 56 (6 1997), pp. 6633–6655.

[20] G. R. W. Quispel and D. I. McLaren. "A new class of energy-preserving numerical integration methods". In: *J. Phys. A* 41.4 (2008), pp. 045206, 7.

[21] L. G. Rebholz. "An energy- and helicity-conserving finite element scheme for the Navier-Stokes equations". In: *SIAM Journal on Numerical Analysis* 45.4 (2007), pp. 1622–1638.

[22] F. Santambrogio. "{Euclidean, metric, and Wasserstein} gradient flows: an overview". In: *Bulletin of Mathematical Sciences* 7 (2017), pp. 87–154.
[23] N. Sato and M. Yamada. "A reduced ideal MHD system for nonlinear magnetic field turbulence in plasmas with approximate flux surfaces". In: *Journal of Mathematical Physics* 65.9 (2024).

[24] A. van der Schaft. "Port-Hamiltonian systems: an introductory survey". In: *Proceedings of the International Congress of Mathematicians Vol. III*. suppl 2. European Mathematical Society Publishing House, 2006, pp. 1339–1365.

[25] G. Vissio and V. Lucarini. "Mechanics and thermodynamics of a new minimal model of the atmosphere". In: *The European Physical Journal Plus* 135.807 (2020).

[26] O. Walsh. "Eddy solutions of the Navier-Stokes equations". In: *The Navier-Stokes Equations II — Theory and Numerical Methods*. Ed. by J. G. Heywood, K. Masuda, R. Rautmann, and V. A. Solonnikov. Berlin, Heidelberg: Springer Berlin Heidelberg, 1992, pp. 306–309.

[27] Y. Zhang, A. Palha, M. Gerritsma, and Q. Yao. "A MEEVC discretization for two-dimensional incompressible Navier-Stokes equations with general boundary conditions". In: *Journal of Computational Physics* 510 (2024), p. 113080.