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 portHamiltonian 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 semidefinite 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 infinitedimensional systems, including ordinary and partial differential equations and differentialalgebraic 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 skewadjoint 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 timestepping 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 portHamiltonian 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 structurepreserving 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 righthand 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 semidiscretization 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 structurepreserving 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 finitedimensional 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 structurepreserving discretization.
As a second step of our analysis, we will show that the energy identities in the integral form (1.3) remain valid after timediscretization 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 nonsmooth 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 semidiscretization in space, equation (2.1) amounts to a finitedimensional system of ordinary or differentialalgebraic 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 wellposedness of the discretized problem but rather derive properties of semidiscrete 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 nonconforming 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 semidiscretized 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 wellposedness, but rather prove properties of discrete solutions.
Remark 3.
It suffices to consider scalarvalued 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 righthand 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 semidiscretization 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 antisymmetric, 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 antisymmetry 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 lefthand 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 scalarvalued 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 lefthand 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 LobattoIIIA 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 wellknown timediscretization 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 MagnetoQuasistatics
As a first test problem, we consider the equations of magnetoquasistatics, 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 magnetoquasistatic 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 semidiscretization which automatically inherits the energy identity derived above. After choice of a basis for the space
𝕍
h
, the semidiscrete 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 energydissipation 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 righthand 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 semidiscrete 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 finitedimensional 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 realvalued 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 timediscretization 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 wellknown driftoff phenomenon; see, e.g., [14, Chapter 5] or [29, Chapter 14] for details. Let us note that some artificial driftoff 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 driftoff in index2 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 MagnetoQuasistatics
We consider a geometry
Ω
=
Ω
x
y
×
ℝ
with infinite extent in the zdirection 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 magnetoquasistatics 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.
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 LobattoIIIA 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 energydissipation 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).
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 differentialalgebraic 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 index1 reformulation of the original problem. In compact notation, the above system can be written as
𝒞
(
u
)
∂
t
u
=

ℋ
′
(
u
)
with antisymmetric matrix
𝒞
(
u
)
and energy derivative vector
ℋ
′
(
u
)
. The above index1 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 timediscretization, 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.
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 driftoff 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 timediscretization. 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 nonquadratic 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: Eg331/21
References
[1] B. M. Afkham and J. S. Hesthaven, Structurepreserving modelreduction 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 LargeScale 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 structurepreserving model reduction, Systems Control Lett. 58 (2009), no. 3, 225–232. Search in Google Scholar
[6] C. Beattie and S. Gugercin, Structurepreserving model reduction for nonlinear portHamiltonian systems, 50th IEEE Conference on Decision and Control and European Control Conference (CDCECC), IEEE Press, Piscataway (2011), 6564–6569. Search in Google Scholar
[7] P. Benner, V. Mehrmann and D. C. Sorensen, Dimension Reduction of LargeScale 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, Energypreserving methods on Riemannian manifolds, Math. Comp. 89 (2020), no. 322, 699–716. Search in Google Scholar
[10] E. Celledoni and E. H. Hoiseth, Energypreserving and passivityconsistent numerical discretization of portHamiltonian systems, preprint (2017), https://arxiv.org/abs/1706.08621. Search in Google Scholar
[11] S. Chaturantabut, C. Beattie and S. Gugercin, Structurepreserving model reduction for nonlinear portHamiltonian systems, SIAM J. Sci. Comput. 38 (2016), no. 5, B837–B865. Search in Google Scholar
[12] D. Cohen and E. Hairer, Linear energypreserving 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. EichSoellner 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, Energypreserving 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, Energydiminishing 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 DifferentialAlgebraic 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, DifferentialAlgebraic 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 portHamiltonian 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
Received: 20200229
Revised: 20200827
Accepted: 20201116
Published Online: 20201204
Published in Print: 20210401
© 2021 Walter de Gruyter GmbH, Berlin/Boston