On the Energy Stable Approximation of Hamiltonian and Gradient Systems

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

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

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,

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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 .