# On the Energy Stable Approximation of Hamiltonian and Gradient Systems

Herbert Egger, Oliver Habrich and Vsevolod Shashkov

# Abstract

A general framework for the numerical approximation of evolution problems is presented that allows to preserve an underlying dissipative Hamiltonian or gradient structure exactly. The approach relies on rewriting the evolution problem in a particular form that complies with the underlying geometric structure. The Galerkin approximation of a corresponding variational formulation in space then automatically preserves this structure which allows to deduce important properties for appropriate discretization schemes including projection based model order reduction. We further show that the underlying structure is preserved also under time discretization by a Petrov–Galerkin approach. The presented framework is rather general and allows the numerical approximation of a wide range of applications, including nonlinear partial differential equations and port-Hamiltonian systems. Some examples will be discussed for illustration of our theoretical results, and connections to other discretization approaches will be highlighted.

## 1 Introduction

The modeling of dynamical systems often leads to problems with a generalized Hamiltonian or gradient structure and various examples have been studied intensively in the literature. In this paper, we consider abstract evolution problems of the form

(1.1) 𝒞 ( u ) t u = - ( u ) + f ( u ) ,

which include Hamiltonian and gradient systems as special cases. Here : 𝕍 represents a continuously differentiable energy or storage functional with derivative denoted by : 𝕍 𝕍 and 𝕍 . The functions f : 𝕍 𝕍 and 𝒞 : 𝕍 ( 𝕍 , 𝕍 ) are assumed to be sufficiently smooth, defined on some Banach space 𝕍 with dual 𝕍 , and ( 𝕍 , 𝕍 ) denoting the space of linear bounded operators mapping from 𝕍 to 𝕍 . We further assume that 𝒞 ( u ) : 𝕍 𝕍 is positive semi-definite for any u 𝕍 . For every point in time, equation (1.1) is thus understood in the sense of linear functionals in 𝕍 . As we will see below, the particular form (1.1), with the derivative of the Hamiltonian appearing explicitly in the equation, is canonical for the approximation by Galerkin methods. One might also replace 𝒞 ( u ) t u and f ( u ) by more general expressions 𝒞 ( u , t u ) and f ( u , t u ) .

Equation (1.1) is rather general and may represent finite- or infinite-dimensional systems, including ordinary and partial differential equations and differential-algebraic problems. Questions of existence and uniqueness of solutions can hardly be addressed on this general level. However, properties of sufficiently regular solutions can be derived, if they exist. Any classical solution u C 1 ( 0 , T ; 𝕍 ) of system (1.1), for instance, can be shown to satisfy

(1.2) d d t ( u ( t ) ) = ( u ( t ) ) , t u ( t ) = f ( u ( t ) ) , t u ( t ) - 𝒞 ( u ( t ) ) t u ( t ) , t u ( t )

for all 0 t T ; here , denotes the duality product on 𝕍 × 𝕍 . The physical meaning of this energy–dissipation identity is that the total energy ( u ) of the system changes in time only due to the work done by external forces and internal dissipation. Through integration in time, one can obtain a corresponding identity in integral form

(1.3) ( u ( t ) ) = ( u ( s ) ) + s t f ( u ( r ) ) , t u ( r ) - 𝒞 ( u ( r ) ) t u ( r ) , t u ( r ) d r ,

which holds for all 0 s t T and which remains valid for less regular solutions, e.g. u W 1 , p ( 0 , T ; 𝕍 ) ; in this case, equation (1.2) is valid only for a.e. 0 t T . Let us note that, in the derivation of (1.2) and (1.3), only the variational form

(1.4) 𝒞 ( u ( t ) ) t u ( t ) , v = - ( u ( t ) ) , v + f ( u ( t ) ) , v

of the problem with v = t u ( t ) was used and that the validity of (1.4) for all v 𝕍 yields an equivalent variational or weak form of the system (1.1) under consideration.

We briefly mention two important particular cases of problem (1.1) that will be covered automatically by our considerations. If 𝒞 ( u ) = - 𝒞 ( u ) * is skew-adjoint and f ( u ) 0 , then equation (1.1) models a Hamiltonian system and the energy ( u ) is preserved for all time. On the other hand, if 𝒞 ( u ) = 𝒞 ( u ) * is positive definite and f ( u ) 0 , then the problem describes a gradient system and energy decays until a steady state is reached.

Energy identities (1.2) and (1.3) play a fundamental role in the analysis of evolution problem (1.1), and they encode important physical principles. Therefore, much research has been devoted to the construction and analysis of numerical methods that still satisfy similar identities even after discretization. In [19], the discrete derivative methods, also often called discrete gradient methods, were introduced, which, applied to (1.1), lead to time-stepping schemes of the form

𝒞 ¯ ( u n , u n - 1 ) u n - u n - 1 τ = - ¯ ( u n , u n - 1 ) + f ¯ ( u n , u n - 1 ) , n 1 .

Particular approximations 𝒞 ¯ ( u n , u n - 1 ) , ¯ ( u n , u n - 1 ) , f ¯ ( u n , u n - 1 ) result in the discrete derivative and the average vector field methods studied in [19, 26]; second order convergence in time was established for these methods. The average vector field collocation methods, which were proposed in [12, 20] for the numerical integration of Hamiltonian systems, yield corresponding methods of higher order. Their application to port-Hamiltonian systems was investigated in [10], and the generalization to gradient and Hamiltonian systems on Riemannian manifolds was studied in [21, 8, 9]. The general ideas underlying the discrete gradient approach can also be utilized for the space discretization of nonlinear evolution problems; we refer to the monograph [17] for details. Many further results concerning the structure-preserving discretization of Hamiltonian systems can be found in [25, 23] as well as the references given there.

In a recent paper [13], we studied the systematic approximation of dissipative nonlinear dynamical systems of another canonical form, namely

𝒬 ( u ) t u = 𝒜 ( u ) ,

by means of Galerkin approximation in space and discontinuous Galerkin discretization in time. In this system, the operator 𝒬 ( u ) is related to the derivative of the energy functional and 𝒜 ( u ) encodes the dissipation terms, which essentially drive the evolution. In the current manuscript, we follow a similar route but address problems of a different form, namely (1.1), in which dissipative terms are encoded in the operator 𝒞 ( u ) in front of the time derivative, while the derivative ( u ) of the energy functional appears as a driving force on the right-hand side.

As a first step of our analysis, we will show that the particular structure of the problem is preserved exactly by Galerkin projection of the corresponding variational formulation (1.4) in space. Such discretization strategies are frequently used for the semi-discretization of evolutionary partial differential equations, see e.g. [16, 30], and also for model order reduction of dynamical systems [4, 7]. Let us mention [1, 5, 6, 11, 28] as particular examples addressing the structure-preserving model reduction of dissipative Hamiltonian systems. Note that the structure preservation under Galerkin projection in space does in general not hold for other, seemingly equivalent formulations of the evolution problem, like

(1.5) t u = ( u ) u ( u ) + g ( u ) ,

which are often considered in the literature; see [10, 19, 20, 26] for instance. In a finite-dimensional setting, the gradient u ( u ) may be identified with the derivative ( u ) . If ( u ) is invertible, we may transform this equation into

(1.6) - ( u ) - 1 t u = - u ( u ) - ( u ) - 1 g ( u ) ,

which, apart from the identification of u ( u ) with ( u ) is of the particular form (1.1) with 𝒞 ( u ) = - ( u ) - 1 and f ( u ) = - ( u ) - 1 g ( u ) . The equivalence of (1.5) and (1.6), however, is usually lost after approximation. To see this, assume g ( u ) = 0 for simplicity, and consider the Galerkin approximation in space realized by some projection 𝒫 . Then projection of (1.5) leads to 𝒫 ( u ) - 1 t u = 𝒫 u ( u ) , while that of (1.6) yields 𝒫 t u = 𝒫 ( u ) u ( u ) . Equivalence of the two formulations would require that ( u ) commutes with 𝒫 , which is usually not the case. As we will see below, the particular form (1.1) of the problem allows for a natural approximation by Galerkin projection and therefore seems canonical for a structure-preserving discretization.

As a second step of our analysis, we will show that the energy identities in the integral form (1.3) remain valid after time-discretization of (1.1) or its variational form (1.4) by certain Petrov–Galerkin methods. This is again strongly linked to the particular structure of the evolution problem under investigation, which should therefore be taken into account already in the modeling stage. We will demonstrate by examples that this can be achieved without difficulties for typical problems of rather different origin. We will also discuss in detail the connection of our approach to the discrete derivative and the average vector field collocation methods mentioned earlier, which can in fact be viewed as special instances or inexact realizations of the methods proposed in this paper. Our considerations may therefore provide further insight also into the analysis and generalization of these methods and the construction of new schemes.

The remainder of the manuscript is organized as follows. In Section 2, we discuss the space discretization of problem (1.1) by Galerkin approximation of the corresponding variational form (1.4). In addition, we discuss inexact variants of the resulting methods, which may be more convenient for a practical realization. In Section 3, we then study the time discretization by a Petrov–Galerkin approach, and we highlight the connection to other methods that have been discussed in the literature. Again, some level of inexactness is allowed in the realization which may facilitate the actual numerical implementation. The wide applicability of our methods will be demonstrated in Section 4, where we discuss some typical test problems in finite and infinite dimensions. Corresponding numerical results are presented in Section 5 for illustration of our theoretical results.

## 2 Space Discretization

Let 𝕍 h 𝕍 be a closed subspace of the Hilbert space 𝕍 . We consider the Galerkin approximation of the variational principle (1.4) in 𝕍 h given by

(2.1) 𝒞 ( u h ( t ) ) t u h ( t ) , v h = - ( u h ( t ) ) , v h + f ( u h ( t ) ) , v h .

This equation is assumed to hold for all v h 𝕍 h and for all 0 t T . As a direct consequence of the particular problem structure, we obtain the following structure result.

## Theorem 1.

Let u h C 1 ( [ 0 , T ] ; V h ) satisfy (2.1) for all 0 t T and v h V h . Then the discrete energy dissipation relation

(2.2) d d t ( u h ( t ) ) = f ( u h ( t ) ) , t u h ( t ) - 𝒞 ( u h ( t ) ) t u h ( t ) , t u h ( t )

holds for all 0 t T , and as a direct consequence, one also has the integral identity

( u h ( t ) ) = ( u h ( s ) ) + s t f ( u h ( r ) ) , t u h ( r ) - 𝒞 ( u h ( r ) ) t u h ( r ) , t u h ( r ) d r

for all 0 s t T . The second identity remains valid for non-smooth solutions, e.g., if u h H 1 ( 0 , T ; V h ) , while the first identity holds only for a.e. 0 t T in that case.

## Proof.

The first identity follows by formal differentiation of ( u ( t ) ) with respect to time and use of identity (2.1) with v h = t u h ( t ) ; this is possible since t u h ( t ) 𝕍 h is an admissible test function. The second identity then simply follows by integration of the first relation with respect to time. ∎

After semi-discretization in space, equation (2.1) amounts to a finite-dimensional system of ordinary or differential-algebraic equations. Sufficient conditions on 𝒞 ( u ) and ( u ) guaranteeing the existence and uniqueness of solutions for the corresponding initial value problems could be stated, but for sake of generality, we do not address the well-posedness of the discretized problem but rather derive properties of semi-discrete solutions, assuming that they exist.

## Remark 1.

Let us emphasize that the structure preservation under Galerkin projection is a direct consequence of the particular canonical form of the evolution problem under consideration. This would not be the case if the discretization was based on a different formulation; compare with (1.5) and the corresponding remarks in the introduction.

## Remark 2.

A similar result can be obtained if certain approximations for the individual terms in the discrete variational problem (2.1) are used. As an example, we may consider an inexact Galerkin approximation of the form

𝒞 h ( u h ( t ) ) t u h ( t ) , v h = - h ( u h ( t ) ) , v h + f h ( u h ( t ) ) , v h .

Energy identities similar to (2.2) or (1.3) then hold with ( u ) , 𝒞 ( u ) and f ( u ) replaced appropriately by h ( u h ) , 𝒞 h ( u h ) and f h ( u h ) . One could also utilize a discretization dependent duality product , h on 𝕍 h × 𝕍 h and consider non-conforming Galerkin approximations with spaces 𝕍 h 𝕍 ; we refer to [17] for further considerations in this direction and to Sections 4 and 5 for particular examples.

## 3 Time Discretization

We now turn to the time discretization. Since the variational form (1.4) of the problem is inherited by Galerkin approximation in space, see the previous section, the following arguments also cover problems that have already been semi-discretized in space.

Let T N = { 0 = t 0 < < t N = T } be a partition of the interval [ 0 , T ] ; set τ n = t n - t n - 1 , and denote by P k ( T N ; 𝕍 ) , k 0 , the space of piecewise polynomial functions over the partition T N with values in 𝕍 . For the time discretization of (1.4), we consider the following Petrov–Galerkin approach: find

u N P k + 1 ( T N ; 𝕍 ) H 1 ( 0 , T ; 𝕍 )

satisfying

(3.1) t n - 1 t n 𝒞 ( u N ( t ) ) t u N ( t ) , v ~ N ( t ) d t = t n - 1 t n f ( u N ( t ) ) , v ~ N ( t ) - ( u N ( t ) ) , v ~ N ( t ) d t

for all v ~ N P k ( T N ; 𝕍 ) and all time slabs 1 n N . We note that the ansatz functions u N are piecewise polynomials of degree k + 1 and globally continuous in time, while the test functions v ~ N are piecewise polynomials of degree k and may be discontinuous at time points t n , n = 1 , , N - 1 ; therefore, (3.1) is a Petrov–Galerkin approximation [2, 30]. Solvability of system (3.1) can again be certainly obtained under appropriate structural assumptions on the operators 𝒞 ( u ) and the functional ( u ) , which however depend on the particular examples to be studied. We again stay on the general level and do not consider issues of well-posedness, but rather prove properties of discrete solutions.

## Remark 3.

It suffices to consider scalar-valued test functions v ~ N P k ( T h ; ) in the discrete variational problem (3.1), in which case one obtains the equivalent formulation

(3.2) t n - 1 t n f ( u N ( t ) ) v ~ N ( t ) - 𝒞 ( u N ( t ) ) t u N ( t ) v ~ N ( t ) d t = t n - 1 t n ( u N ( t ) ) v ~ N ( t ) d t ,

which now has to be understood as an equation in 𝕍 ; compare with system (1.1) and its equivalent variational formulation (1.4) on the continuous level. We thus obtain a time discretization of (1.1) as well. This observation will be useful for our discussions later on.

With arguments similar to the previous sections, we now obtain the following result.

## Theorem 2.

Let u N P k + 1 ( T N ; V ) H 1 ( 0 , T ; V ) solve (3.1) for all v ~ N P k ( T N ; V ) and all 1 n N (or equivalently (3.2) for all v ~ N P k ( T N ; R ) ). Then

( u N ( t n ) ) = ( u N ( t m ) ) + t m t n f ( u N ( t ) ) , t u N ( t ) - 𝒞 ( u N ( t ) ) t u N ( t ) , t u N ( t ) d t

for all time instances 0 t m t n T of the respective time grid T N .

## Proof.

From the fundamental theorem of calculus, we can deduce

( u N ( t n ) ) = ( u N ( t n - 1 ) ) + t n - 1 t n ( u N ( t ) ) , t u N ( t ) d t = (i) .

By choice of the trial and test spaces, the time derivative v ~ N = t u N is an admissible test function for problem (3.1), and we can replace the right-hand side (i) by

(i) = t n - 1 t n f ( u N ( t ) ) , t u N ( t ) - 𝒞 ( u N ( t ) ) t u N ( t ) , t u N ( t ) d t .

This proves the result for m = n - 1 ; the general case then follows by induction over n. ∎

## Remark 4.

Similar to the corresponding result for the semi-discretization in space, the proof and validity of the discrete energy identity above strongly relies on the particular structure of problem (1.1) and its variational formulation (1.4). Further, note that the energy identity here only holds at specific points in time and, in general, we do not have a pointwise energy identity like (1.2) or (2.2) for the time discrete approximation.

Let us next comment on the connection to other approximation schemes that have been proposed for the time discretization of Hamiltonian and gradient systems.

## Remark 5.

Consider the case k = 0 in the Petrov–Galerkin method (3.1). Then u N is piecewise linear in time, and v ~ N is piecewise constant. Using the abbreviation u n = u N ( t n ) , scheme (3.1) can be written equivalently as

𝒞 ¯ ( u n , u n - 1 ) u n - u n - 1 τ n = - ¯ ( u n , u n - 1 ) + f ¯ ( u n , u n - 1 ) ,

with averages

𝒞 ¯ ( u n , u n - 1 ) = 1 τ n t n - 1 t n 𝒞 ( u n ( s ) ) d s , ¯ ( u n , u n - 1 ) = 1 τ n t n - 1 t n ( u n ( s ) ) d s , f ¯ ( u n , u n - 1 ) = 1 τ n t n - 1 t n f ( u n ( s ) ) d s ;

here we used u n ( s ) = u n - 1 + t - t n - 1 τ n ( u n - u n - 1 ) for abbreviation. For k = 0 , method (3.1) thus coincides with a discrete gradient method as outlined in the introduction; see [19, 26, 10, 21] for further details.

## Remark 6.

If f ( u ) = 0 and 𝒞 ( u ) is anti-symmetric, then exact energy conservation is implied by Theorem 2. In this case, the proof essentially relies on the exact realization of the basic identity

( u N ( t n ) ) - ( u N ( t n - 1 ) ) = t n - 1 t n ( u N ( t ) ) , t u N ( t ) d t

and the anti-symmetry of the operator 𝒞 ( u ) , which allows for some degree of inexactness in the realization. Using a quadrature rule ( γ j , ω j ) , 0 j J , to approximate the integral on the left-hand side of (3.2), for instance, leads to

t n - 1 t n 𝒞 ( u N ( t ) ) t u N ( t ) v ~ N ( t ) d t τ n j = 0 J 𝒞 ( u N ( t j n ) ) t u N ( t j n ) v ~ N ( t j n ) ω j ,

with intermediate time points t j n = t n - 1 + γ j ( t n - t n - 1 ) . For convenience, we here employed formulation (3.2) with scalar-valued test functions v ~ N P k ( T N ; ) . A similar approximation can be used for the term t n - 1 t n f ( u N ( t ) ) v ~ N ( t ) d t if f ( u ) 0 . Note, however, that an exact integration of the nonlinear term involving ( u N ) is always required to guarantee exact energy conservation. When choosing v ~ N ( t ) as Lagrange polynomials for the quadrature points t j n , j = 0 , , J = k , this immediately leads to an average vector field collocation method discussed in [10, 20, 21]. If 𝒞 ( u ) = 𝒞 is independent of u and t, and f ( u ) = f ( t ) is a polynomial of degree k + 1 in time, then exact integration of the left-hand side of (3.2) can be achieved by the Gauss quadrature rule. In that case, the Petrov–Galerkin approximations and the average vector field collocation methods coincide. The Petrov–Galerkin scheme of order k is actually closely related to the Lobatto-IIIA method with s = k + 1 stages, which can be obtained by appropriate inexact numerical integration of all terms in (3.1). The proposed Petrov–Galerkin approach can thus be used to derive well-known time-discretization methods in a systematic manner, and it provides a framework to generalize these methods to a wider class of problems.

## 4 Examples

We now illustrate the viability of the proposed discretization strategies in space and time by discussing their application for some typical problems we have in mind.

### 4.1 Magneto-Quasistatics

As a first test problem, we consider the equations of magneto-quasistatics, which arise in the eddy current approximation of Maxwell’s equations at low frequency [3]. The magnetic flux density 𝐛 = curl 𝐚 is then usually represented by a magnetic vector potential 𝐚 , which is assumed to satisfy

σ t 𝐚 + curl ( ν ( curl 𝐚 ) curl 𝐚 ) = 𝐣 in Ω .

Here σ denotes the electric conductivity, ν is the incremental magnetic reluctivity tensor, and 𝐣 a given source current density, and Ω 3 some bounded Lipschitz domain.

For ease of presentation, we further assume that σ is bounded and uniformly positive. Then 𝐚 is determined uniquely without gauging, and 𝐞 = - t 𝐚 amounts to the electric field. Moreover, we consider homogeneous boundary conditions 𝐚 × 𝐧 = 0 on Ω , where 𝐧 is the outward pointing unit normal vector on Ω . The variational formulation for the above problem then reads

Ω σ t 𝐚 ( t ) 𝐰 d x + Ω ( ν ( curl 𝐚 ( t ) ) curl 𝐚 ( t ) ) curl 𝐰 d x = Ω 𝐣 ( t ) 𝐰 d x ,

which is supposed to hold for all t of relevance and for all functions 𝐰 H 0 ( curl ; Ω ) having a square integrable weak curl and satisfying zero boundary conditions 𝐰 × 𝐧 = 0 on Ω ; see [3, 27] for details. We assume that there exists a magnetic energy density η : 3 such that 𝐛 η ( 𝐛 ) = ν ( 𝐛 ) 𝐛 ; in other words,

η ( 𝐛 ) = 0 𝐛 ( ν ( 𝐛 ) 𝐛 ) 𝑑 𝐛 ;

for ν ( 𝐛 ) = ν 0 , with scalar constant ν 0 and denoting the identity matrix, we would obtain the usual term η ( 𝐛 ) = ν 0 2 | 𝐛 | 2 . The expression ( 𝐚 ) = Ω η ( curl 𝐚 ) d x now defines the magnetic energy of the system. By differentiation in time, we obtain

d d t ( 𝐚 ( t ) ) = Ω ( ν ( curl 𝐚 ( t ) ) curl 𝐚 ( t ) ) curl t 𝐚 ( t ) d x = Ω 𝐣 ( t ) t 𝐚 ( t ) d x - Ω σ t 𝐚 ( t ) t 𝐚 ( t ) d x ,

which expresses the fact that the magnetic energy of the system is only altered by the power supplied through excitation currents and the dissipation caused by eddy currents.

By setting 𝕍 = H 0 ( curl ; Ω ) and u = 𝐚 , the magneto-quasistatic problem can be seen to fit into the abstract framework of the previous sections, with ( u ) = Ω η ( 𝐚 ) d x , derivative ( u ) = curl ( ν ( curl 𝐚 ) curl 𝐚 ) , which has to be understood in a distributional sense, generalized time derivative 𝒞 ( u ) t u = σ t 𝐚 and forcing term f ( u ) = 𝐣 .

We then consider the Galerkin approximation in space using an appropriate finite element space 𝕍 h H 0 ( curl ; Ω ) ; see [27] for details. Following our previous considerations, this leads to a semi-discretization which automatically inherits the energy identity derived above. After choice of a basis for the space 𝕍 h , the semi-discrete problem can be cast into a system of ordinary differential equations

M σ t a ( t ) + K ν ( a ( t ) ) a ( t ) = j ( t ) ,

and the time discretization of this system by a Petrov–Galerkin approximation, as proposed in Section 3, leads to a fully discrete scheme which satisfies a corresponding energy identity in integral form (1.3) at discrete points t n of the time mesh. We thus obtain a discretization scheme which automatically inherits the physical principle of energy conservation.

### Remark 7.

Before closing this section, let us briefly comment on some obvious generalizations. Without any complications, one can consider other types of boundary conditions and a field dependent conductivity σ = σ ( t 𝐚 ) . If σ vanishes identically on a subdomain Ω nc Ω , then one has to restrict 𝐚 in Ω nc by some gauging conditions [3, 27]; these can be treated, e.g., as additional constraints with arguments similar to those in Section 4.2 below.

### 4.2 Cahn–Hilliard Equation

A simple mathematical model for the phase separation in binary fluids is given by the Cahn–Hilliard equation; see e.g. [15]. After elimination of the chemical potential, the simplest form of the problem can be phrased as

t u = Δ ( - γ Δ u + ψ ( u ) ) in Ω , t > 0 .

Here u represents the difference of the phase fractions of the two components, ψ is a double well potential with two minima, and γ > 0 is a positive constant. When complementing the system with homogeneous Neumann boundary conditions

0 = n u = n Δ u on Ω , t > 0 ,

which implies that n ψ ( u ) = ψ ′′ ( u ) n u = 0 , one can see that d d t Ω u ( t ) d x = Ω t u ( t ) d x = 0 , i.e., the total mass in the system does not change over time. For ease of notation, we will assume that Ω u ( t ) d x = 0 in the sequel. We further denote by ( - Δ ) - 1 the solution operator for the Neumann problem

- Δ w = f in Ω , n w = 0 on Ω .

For any sufficiently regular f with zero average, this equation has a unique weak solution w H 1 ( Ω ) with Ω w d x = 0 . By formally applying this operator to the Cahn–Hilliard equation stated above, we obtain the equivalent system

( - Δ ) - 1 t u = γ Δ u - ψ ( u ) in Ω ,
0 = n u on Ω ,

which will be the basis for further considerations. The weak form of this problem reads

( - Δ ) - 1 t u ( t ) , v = - γ u ( t ) , v - ψ ( u ( t ) ) , v ,

which is assumed to hold for any sufficiently regular test function v and for any time t under consideration. Here a , b = Ω a b d x amounts to the L 2 -scalar product on Ω.

We next recall that, besides the conservation of mass, the solutions of the Cahn–Hilliard problem also are governed by decay of the free energy ( u ) = Ω γ 2 | u | 2 + ψ ( u ) d x . By formal differentiation in time, we obtain

d d t ( u ( t ) ) = γ u ( t ) , t u ( t ) + ψ ( u ( t ) ) , t u ( t ) = - ( - Δ ) - 1 t u ( t ) , t u ( t ) .

In the second step, we used the weak form of the problem with test function v = t u ( t ) . Due to the Poincaré inequality, the operator ( - Δ ) - 1 is positive definite on the space of functions with vanishing average, and hence the free energy ( u ( t ) ) is monotonically decreasing until the system reaches a steady state. The Cahn–Hilliard problem in its second equivalent form thus has exactly the canonical structure (1.1), with 𝒞 ( u ) = ( - Δ ) - 1 independent of u, f ( u ) 0 , and 𝕍 = { v H 1 ( Ω ) : Ω v d x = 0 } . Discretizations which inherit the underlying energy-dissipation structure can thus be formally obtained with the arguments outlined in Sections 2 and 3.

Let us note that a standard Galerkin approximation of the weak form of the Cahn–Hilliard system leads to an approximation u h with values in 𝕍 h 𝕍 characterized by a discrete variational principle of the form

( - Δ - 1 t u h ( t ) , v h = - γ u h ( t ) , v h - ψ ( u h ( t ) ) , v h .

While theoretically sound, such a method is not very useful in practice since the application of the inverse Laplacian ( - Δ ) - 1 can in general not be realized exactly. To overcome this problem, one may utilize a discrete approximation ( - Δ h ) - 1 , which is defined via the solution w h = ( - Δ h ) - 1 f 𝕍 h of the discretized problem

w h , η h = f , η h for all η h 𝕍 h .

Again, a zero average condition for the right-hand side f and the solution w h allow to guarantee existence of a unique solution. The approximation u h for the Cahn–Hilliard system is then characterized by the discrete variational principle

( - Δ h ) - 1 t u h ( t ) , v h = - γ u h ( t ) , v h - ψ ( u h ( t ) ) , v h ,

which is again required to hold for all v h 𝕍 h and all relevant t. With similar reasoning as on the continuous level, one can verify the energy identity

d d t ( u h ( t ) ) = - ( - Δ h ) - 1 t u h ( t ) , t u h ( t ) ,

which shows that the energy of semi-discrete solutions again decreases until a steady state is reached. This approach corresponds to an inexact Galerkin approximation in space, as discussed in Remark 2. Further approximations, e.g., numerical quadrature for the nonlinear terms, may be used to facilitate the numerical implementation; see Section 5.

For the subsequent time discretization, we can employ a Petrov–Galerkin approximation as outlined in Section 3. Since the dissipation term here is linear and the force term vanishes, this coincides with a particular average vector field collocation method based on Gauss quadrature; see [21], and refer to Remark 6 for discussion in this direction.

### 4.3 Constrained Hamiltonian Systems

As another relevant class of applications, we consider finite-dimensional Hamiltonian systems with holonomic constraints given by

(4.1) t q = H p ( p , q ) ,
(4.2) t p = - H q ( p , q ) + f ( p , q ) + g q ( p , q ) λ ,
(4.3) 0 = g ( q ) .

Such systems arise, for instance, in the modeling of multibody dynamics: see, e.g., [14, Chapter 1], [25, Chapter 7] or [29, Chapter 13]. Here p and q are, respectively, the vectors of generalized coordinates and momenta, H ( p , q ) is the Hamiltonian, i.e., the energy or storage functional, f ( p , q ) denotes external forces, λ is the vector of Lagrange multipliers for the constraints, and g q ( p , q ) λ are the corresponding constraint forces. We assume in the following that p, q, λ are real-valued vectors and use subscripts to denote partial derivatives with respect to these variables. The symbol , then stands for the Euclidean scalar product on n .

By differentiating (4.3) with respect to time, we obtain the linearized constraint

(4.4) 0 = g q ( q ) t q ,

which is equivalent to the nonlinear constraint (4.3) up to a constant factor that can be fixed by an appropriate initial condition. We further introduce the principal function of the Lagrange multiplier Λ ( t ) = 0 t λ ( s ) d s . The variational form of system (4.1), (4.2) complemented by the linearized constraint (4.4) can then be written as

(4.5) t q ( t ) , v = H p ( p ( t ) , q ( t ) ) , v ,
(4.6) t p ( t ) , w = - H q ( p ( t ) , q ( t ) , w + f ( p ( t ) , q ( t ) ) , w + g q ( q ( t ) ) t Λ ( t ) , w ,
(4.7) 0 = g q ( q ( t ) ) t q ( t ) , η .

These identities are assumed to hold for all test vectors v, w, η of appropriate size, and any time t > 0 . For every sufficiently smooth solution ( p , q , Λ ) of (4.5)–(4.7), we then get

d d t H ( p ( t ) , q ( t ) ) = H p ( p ( t ) , q ( t ) ) , t p ( t ) + H q ( p ( t ) , q ( t ) ) , t q ( t ) = t q ( t ) , t p ( t ) - t p ( t ) , t q ( t ) + f ( p ( t ) , q ( t ) ) , t q ( t ) + g q ( q ( t ) ) t Λ ( t ) , t q ( t ) = f ( p ( t ) , q ( t ) ) , t q ( t ) .

In the last step, we used that g q ( q ) t Λ , t q = g q ( q ) t q , t Λ = 0 , which follows from testing equation (4.7) with η = t Λ . This identity describes that the total energy of the system can only change due to the work done by the external forces.

After rearrangement of terms, one can see that system (4.1), (4.2) with linearized constraint (4.4) is again of the abstract form 𝒞 ( u ) t u = - ( u ) + f ( u ) with state vector u = ( p , q , Λ ) , energy functional ( u ) = H ( p , q ) and operators

𝒞 ( u ) = ( 0 - I 0 I 0 - g q ( q ) 0 g q ( q ) 0 ) and f ( u ) = ( 0 f ( p , q ) 0 ) .

Hence all results about the approximation by Galerkin methods in space and the Petrov–Galerkin method in time obtained in Sections 2 and 3 apply immediately. In particular, the Galerkin projection of (4.5)–(4.7) into a conforming subspace 𝕍 h 𝕍 automatically inherits the energy identity stated above. Our results in particular also cover model order reduction approaches based on Galerkin projection [7]. Furthermore, the energy identity in its integral form also remains valid after time-discretization if the appropriate Petrov–Galerkin approximation is used; see Theorem 2 and Remark 3.

### Remark 8.

In the Petrov–Galerkin scheme (3.1), we may use test functions of the form v ~ ( t ) = 1 | [ t n - 1 , t n ] v ¯ with v ¯ 𝕍 arbitrary. For the system under investigation, this yields

0 = t n - 1 t n g q ( q ( t ) ) t q ( t ) d t = t n - 1 t n d d t g ( q ( t ) ) d t = g ( q ( t n ) ) - g ( q ( t n - 1 ) ) .

The original constraint (4.3) thus remains valid for all time points t n of the time grid T N if it was valid at initial time t 0 ; we only used that the piecewise constant functions in time are elements of the test space V ~ N in (3.1). Hence, in contrast to other time discretization schemes, the proposed Petrov–Galerkin approximation of (4.1), (4.2) and (4.4) does not suffer from the well-known drift-off phenomenon; see, e.g., [14, Chapter 5] or [29, Chapter 14] for details. Let us note that some artificial drift-off may be observed in practice due to inexact integration of (3.1), which can be removed by standard stabilization or projection steps; see [14, Chapter 5] and [29, Chapter 14]. An alternative approach to avoid drift-off in index-2 formulations of multibody dynamics is the GGL approach [18], also see [29, Chapter 14], which however requires to extend the system by additional unknowns. Another approach by “minimal extension” is discussed in [24, Section 6.4]. We are not aware, however, of a provably energy conserving time discretization strategy for these approaches.

## 5 Numerical Tests

For illustration of our theoretical results and in order to demonstrate the practicability of the proposed methods, we now briefly discuss their application and numerical realization for the test problems mentioned in the previous section.

### 5.1 Magneto-Quasistatics

We consider a geometry Ω = Ω x y × with infinite extent in the z-direction and assume that 𝐣 = ( 0 , 0 , j z ) and that j z and σ are independent of z. Then solutions of the eddy current problem can be found in the form 𝐚 = ( 0 , 0 , a z ) with a z independent of z, and the equations of magneto-quasistatics reduce to

σ t a z = div ( ν ( a z ) a z ) - j z in Ω x y .

We assume that Ω x y 2 is bounded and complement the system by homogeneous boundary conditions; see Figure 1 for an illustration of the geometric arrangement.

As magnetic reluctivity function for our tests, we take ν ( 𝐛 ) = ν 0 ( 1 - α β + | 𝐛 | 2 ) with constants 0 < α < β leading to the energy density u ( 𝐛 ) = ν 0 2 ( | 𝐛 | 2 - α log ( β + | 𝐛 | 2 ) ) . The excitation currents are defined as

j z ( x , y , t ) = I ( χ Ω + - χ Ω - ) sin ( ω t )

with I, ω prescribed, and Ω ± denoting the support of positive and negative part of the impressed current; see Figure 1. In our simulations, we use the parameters α = 0.9 , β = 1 , σ = 1 , ω = 2 π and I = 1 . For discretization of the problem, we chose a Galerkin approximation in space with piecewise linear finite elements and a Petrov–Galerkin approximation of orders k = 1 , 2 in time with constant time step size τ n = τ > 0 .

Figure 1

(a) Geometric setup of the model problem. (b) Potential a z and magnetic flux 𝐛 (arrows) for t = 1 / 2 . (c) Energy H ( t ) (red, dashed) at time t and total supply H supp ( t ) (black, dotted) as well as total loss H loss ( t ) (blue, solid) up to time t.

### (c)

In the middle of Figure 1, we show the potential a z and the magnetic field vectors 𝐛 at time t = 1 2 . In the right plot, we display the evolution of the total energy H n = ( a h ( t n ) ) as well as the time accumulated supplied and dissipated energy of the system. Note that according to the discrete energy balance

H n = H 0 + H supp n - H loss n ,

with symbols H supp n = τ k = 1 n P supp k and H loss n = τ k = 1 n P loss k denoting the total accumulated supply and loss energy, respectively, where

P supp k = 1 τ t k - 1 t k Ω x y j z ( t ) t a h ( t ) d x d t and P loss k = 1 τ t k - 1 t k Ω x y σ | t a h ( t ) | 2 d x d t

are the average supply and loss power at time step k. For a better evaluation, we list in Table 1 the maximal discrepancy in the above energy identity obtained with the proposed Petrov–Galerkin time discretization and the Crank–Nicolson scheme.

Table 1

Energy error max n | H n - H supp n + H loss n | obtained for the Lobatto-IIIA method with s = 2 and s = 3 stages and the corresponding Petrov–Galerkin methods of order k = 1 and k = 2 .

 τ Lobatto ( s = 2 ) PG ( k = 1 ) Lobatto ( s = 3 ) PG ( k = 2 ) 0.1 1.1000 ⋅ 10 - 3 8.7137 ⋅ 10 - 14 9.2094 ⋅ 10 - 6 1.7073 ⋅ 10 - 13 0.05 2.7303 ⋅ 10 - 4 5.9014 ⋅ 10 - 13 1.6664 ⋅ 10 - 6 1.4951 ⋅ 10 - 12 0.01 1.0740 ⋅ 10 - 5 1.4769 ⋅ 10 - 12 2.0862 ⋅ 10 - 9 2.2324 ⋅ 10 - 12

As predicted, the Petrov–Galerkin method preserves the energy identity on the discrete level up to roundoff errors which may accumulate with the number of time steps. In contrast to that, the Crank–Nicolson scheme is not able to exactly reproduce the underlying energy-dissipation structure on the discrete level.

### 5.2 Cahn–Hilliard Equation

As a second test problem, we consider the Cahn–Hilliard equation with a logarithmic potential

ψ ( x ) = θ 2 [ ( 1 + x ) ln ( 1 + x 2 ) + ( 1 - x ) ln ( 1 - x 2 ) ] + θ c 2 ( 1 - x 2 )

and Neumann boundary conditions; see Section 4.2. For our computations, we use Ω = ( 0 , 1 ) 2 as a computational domain and choose γ = 0.002 , θ = 0.3 , θ c = 0.8 . For discretization, we use the proposed Galerkin Petrov–Galerkin approximation with piecewise linear finite elements in space and piecewise linear approximation in time. The fully discrete solution is characterized by

( - Δ h ) - 1 u h n - u h n - 1 τ n , v h h = - γ u h n + u h n - 1 2 , v h - ψ ( u h n , u h n - 1 ) , v h h .

Here we used

u h n = u h , N ( t n ) and ψ ¯ ( u h n , u h n - 1 ) = 1 τ t n - 1 t n ψ ( u h ( t n - 1 + s τ ) ) d s

for abbreviation. Note that this average vector field approximation automatically results from the Petrov–Galerkin time integration of lowest order. In the terms characterized by , h , we use numerical integration by the vertex rule; see Remark 2. The inverse of the discrete Laplacian is defined by w h = ( - Δ h ) - 1 f with w h , η h = f , η h h for all η h and zero average conditions. The assembly of the inverse Laplacian is circumvented by introducing a discrete potential μ h in the numerical realization; cf. Section 4.2.

For our computations, we utilize a uniform triangulation with mesh size h = 1 64 and a time step τ = 4 10 - 4 . The nonlinear system in every time step is solved by a simplified Newton iteration with tolerance 10 - 12 . The nodal values of the initial datum u h ( 0 ) are chosen randomly around 0 with standard deviation of 10 - 3 .

Figure 2

(a)–(c) Snapshot of the phase fraction u at time t = 0 , 0.3 , 2 . (d) Evolution of the free energy (blue, solid), dissipated energy (black, dotted), total energy (red, dashed).

### (d)

In Figure 2, we display some snapshots of the numerical solution, and we depict the evolution of the discrete free energy, which is defined by

H h n = Ω γ 2 | u h ( t n ) | 2 + I h ψ ( u h ( t n ) ) d x .

Here I h : H 1 ( Ω ) V h denotes the piecewise linear interpolation for which we observe that Ω I h u d x = u , 1 h exactly amounts to numerical integration by the vertex rule. As a consequence, the discrete energy identities derived in Sections 2 and 3 carry over almost verbatim, which yields to the conclusion that the decay of the discrete free energy should be compensated exactly by the total accumulated dissipation. In our computational tests, this can be observed up to machine accuracy.

### 5.3 Constrained Hamiltonian Dynamics

As a simple example for a constrained dynamical system, we consider a pendulum with a mass m = 1 at position ( x , y ) swinging at a fixed distance around the origin. The equations of motion can then be written as

x ˙ = v , v ˙ = - 2 x λ , y ˙ = w , w ˙ = - g - 2 y λ , 0 = x 2 + y 2 - 2 ,

which comprises a differential-algebraic system of index 3; see [22, 24] for details. This example fits into the class of problems considered in Section 4.3 with q = ( x , y ) , p = ( v , w ) , constraint g ( q ) = x 2 + y 2 - 2 , and energy H ( x , y , v , w ) = 1 2 ( v 2 + w 2 ) + g y , consisting of a kinetic and a potential contribution. Following the considerations presented in Section 4.3, we obtain the canonical form of the system, which here has the form

( 0 0 - 1 0 0 0 0 0 - 1 0 1 0 0 0 - 2 x 0 1 0 0 - 2 y 0 0 2 x 2 y 0 ) ( v ˙ w ˙ x ˙ y ˙ Λ ˙ ) = - ( v w 0 g 0 ) .

In the derivation, the algebraic constraint was differentiated once, and the Lagrange multiplier was integrated. The resulting system turns out to be an equivalent index-1 reformulation of the original problem. In compact notation, the above system can be written as 𝒞 ( u ) t u = - ( u ) with anti-symmetric matrix 𝒞 ( u ) and energy derivative vector ( u ) . The above index-1 formulation therefore has the canonical form (1.1) with force vector f ( u ) 0 , and we obtain exact conservation of the energy for all time.

For the time-discretization, we use a Petrov–Galerkin method as outlined in Section 4.3. Since 𝒞 ( u ) and ( u ) both depend affine linearly on u, all time integrals can be evaluated exactly. The approximation of lowest order with piecewise linear functions in time and constant time step τ then leads to a recursion of the form

(5.1) 𝒞 ¯ ( u n , u n + 1 ) u n + 1 - u n τ = - ¯ ( u n , u n + 1 ) , n 0 ,

with 𝒞 ¯ ( u n , u n + 1 ) = 1 2 ( 𝒞 ( u n ) + 𝒞 ( u n + 1 ) ) and ¯ ( u n , u n + 1 ) = 1 2 ( ( u n ) + ( u n + 1 ) ) . As before, u n = u ( t n ) denote the values of the discrete solution at the time points t n .

Figure 3

(a) Snapshot of the position ( x , y ) at time t = 0 , 0.2 , , 3 . (b) Evolution of kinetic (blue, dashed), potential (red, dotted) and total energy (black, solid) of the system. (c) Evolution of the position x (blue, dashed), y (red, dotted) and length (black, solid) of the pendulum.

### (c)

For our numerical tests, we chose the initial conditions as x 0 = - 1 , y 0 = 1 , v 0 = 0 , w 0 = 0 and Λ 0 = 0 . We further set g = 1 and choose a time step size of τ = 0.2 . The nonlinear systems at every time step are solved with the Newton iteration. Three to four iterations were usually required to achieve the prescribed tolerance of 10 - 12 .

In Figure 3, we display some snapshots of the numerical solution and the evolution of the energy and the position as well as length over time. As predicted by our theoretical results, scheme (5.1) is exactly energy preserving and, despite the fact that the constraint was differentiated, does not suffer from the drift-off phenomenon; see Remark 8. In our computations, the length of the pendulum as well as the total energy of the system are conserved up to machine precision even over hundreds of periods.

## 6 Discussion

We proposed an abstract framework for the numerical approximation of evolution problems with an underlying Hamiltonian or gradient structure, and we showed that this underlying geometric structure can be preserved exactly under discretization with Galerkin methods in space and Petrov–Galerkin approximation in time if the evolution problem was written in a canonical form already on the continuous level. We also illustrated that some inexactness in the numerical realization of the Galerkin approximations is possible which may facilitate the numerical realization. The discrete derivative and average vector field collocation methods could be interpreted as such inexact realizations of the Petrov–Galerkin time-discretization. The viability of the proposed approach and the practicability of the resulting discretization schemes was demonstrated for some typical test problems. We note that, in case of complex non-quadratic Hamiltonians, the efficient numerical handling of the term ( u ) , v arising in (1.4) and the corresponding discrete equations is somewhat subtle and may deserve further considerations.

Funding source: Deutsche Forschungsgemeinschaft

Award Identifier / Grant number: TRR 146 project C3

Award Identifier / Grant number: TRR 154 project C04

Award Identifier / Grant number: Eg-331/2-1

### References

[1] B. M. Afkham and J. S. Hesthaven, Structure-preserving model-reduction of dissipative Hamiltonian systems, J. Sci. Comput. 81 (2019), no. 1, 3–21. Search in Google Scholar

[2] G. Akrivis, C. Makridakis and R. H. Nochetto, Galerkin and Runge–Kutta methods: Unified formulation, a posteriori error estimates and nodal superconvergence, Numer. Math. 118 (2011), no. 3, 429–456. Search in Google Scholar

[3] A. Alonso Rodríguez and A. Valli, Eddy Current Approximation of Maxwell Equations. Theory, Algorithms and Applications, Model. Simul. Appl. 4, Springer, Milan, 2010. Search in Google Scholar

[4] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems, Adv. Des. Control 6, Society for Industrial and Applied Mathematics, Philadelphia, 2005. Search in Google Scholar

[5] C. Beattie and S. Gugercin, Interpolatory projection methods for structure-preserving model reduction, Systems Control Lett. 58 (2009), no. 3, 225–232. Search in Google Scholar

[6] C. Beattie and S. Gugercin, Structure-preserving model reduction for nonlinear port-Hamiltonian systems, 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC), IEEE Press, Piscataway (2011), 6564–6569. Search in Google Scholar

[7] P. Benner, V. Mehrmann and D. C. Sorensen, Dimension Reduction of Large-Scale Systems, Lect. Notes Comput. Sci. Eng. 45, Springer, Berlin, 2005. Search in Google Scholar

[8] E. Celledoni, S. L. Eidnes, B. Owren and T. R. Ringholm, Dissipative numerical schemes on Riemannian manifolds with applications to gradient flows, SIAM J. Sci. Comput. 40 (2018), no. 6, A3789–A3806. Search in Google Scholar

[9] E. Celledoni, S. L. Eidnes, B. Owren and T. R. Ringholm, Energy-preserving methods on Riemannian manifolds, Math. Comp. 89 (2020), no. 322, 699–716. Search in Google Scholar

[10] E. Celledoni and E. H. Hoiseth, Energy-preserving and passivity-consistent numerical discretization of port-Hamiltonian systems, preprint (2017), https://arxiv.org/abs/1706.08621. Search in Google Scholar

[11] S. Chaturantabut, C. Beattie and S. Gugercin, Structure-preserving model reduction for nonlinear port-Hamiltonian systems, SIAM J. Sci. Comput. 38 (2016), no. 5, B837–B865. Search in Google Scholar

[12] D. Cohen and E. Hairer, Linear energy-preserving integrators for Poisson systems, BIT 51 (2011), no. 1, 91–101. Search in Google Scholar

[13] H. Egger, Structure preserving approximation of dissipative evolution problems, Numer. Math. 143 (2019), no. 1, 85–106. Search in Google Scholar

[14] E. Eich-Soellner and C. Führer, Numerical Methods in Multibody Dynamics, Eur. Consort. Math. Ind., B. G. Teubner, Stuttgart, 1998. Search in Google Scholar

[15] C. M. Elliott, The Cahn–Hilliard model for the kinetics of phase separation, Mathematical Models for Phase Change Problems (Óbidos 1988), Internat. Ser. Numer. Math. 88, Birkhäuser, Basel (1989), 35–73. Search in Google Scholar

[16] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements, Appl. Math. Sci. 159, Springer, New York, 2004. Search in Google Scholar

[17] D. Furihata and T. Matsuo, Discrete Variational Derivative Method, Chapman & Hall/CRC Numer. Anal. Sci. Comput., CRC Press, Boca Raton, 2011. Search in Google Scholar

[18] C. W. Gear, B. Leimkuhler and G. K. Gupta, Automatic integration of Euler–Lagrange equations with constraints, J. Comput. Appl. Math. 12–13 (1985), 77–90. Search in Google Scholar

[19] O. Gonzalez, Time integration and discrete Hamiltonian systems, J. Nonlinear Sci. 6 (1996), no. 5, 449–467. Search in Google Scholar

[20] E. Hairer, Energy-preserving variant of collocation methods, JNAIAM. J. Numer. Anal. Ind. Appl. Math. 5 (2010), no. 1–2, 73–84. Search in Google Scholar

[21] E. Hairer and C. Lubich, Energy-diminishing integration of gradient systems, IMA J. Numer. Anal. 34 (2014), no. 2, 452–461. Search in Google Scholar

[22] E. Hairer, C. Lubich and M. Roche, The Numerical Solution of Differential-Algebraic Systems by Runge–Kutta methods, Lecture Notes in Math. 1409, Springer, Berlin, 1989. Search in Google Scholar

[23] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration, Springer Ser. Comput. Math. 31, Springer, Heidelberg, 2010. Search in Google Scholar

[24] P. Kunkel and V. Mehrmann, Differential-Algebraic Equations. Analysis and Numerical Solution, EMS Textbk. Math., European Mathematical Society, Zürich, 2006. Search in Google Scholar

[25] B. Leimkuhler and S. Reich, Simulating Hamiltonian Dynamics, Cambridge Monogr. Appl. Comput. Math. 14, Cambridge University, Cambridge, 2004. Search in Google Scholar

[26] R. I. McLachlan, G. R. W. Quispel and N. Robidoux, Geometric integration using discrete gradients, R. Soc. Lond. Philos. Trans. Ser. A Math. Phys. Eng. Sci. 357 (1999), no. 1754, 1021–1045. Search in Google Scholar

[27] P. Monk, Finite Element Methods for Maxwell’s Equations, Numer. Math. Sci. Comput., Oxford University, New York, 2003. Search in Google Scholar

[28] R. V. Polyuga, Discussion on: “Passivity and structure preserving order reduction of linear port-Hamiltonian systems using Krylov subspaces” [mr2731511], Eur. J. Control 16 (2010), no. 4, 407–409. Search in Google Scholar

[29] K. Strehmel and R. Weiner, Numerik gewöhnlicher Differentialgleichungen, Teubner Math. Textb., B. G. Teubner, Stuttgart, 1995. Search in Google Scholar

[30] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, 2nd ed., Springer Ser. Comput. Math. 25, Springer, Berlin, 2006. Search in Google Scholar

Revised: 2020-08-27
Accepted: 2020-11-16
Published Online: 2020-12-04
Published in Print: 2021-04-01

© 2021 Walter de Gruyter GmbH, Berlin/Boston