Constructing equivalence-preserving Dirac variational integrators with forces

TL;DR

This paper develops a class of Dirac variational integrators with forces that preserve geometric structures in interconnected mechanical systems, enhancing long-term simulation accuracy.

Contribution

It introduces a novel construction of equivalence-preserving Dirac variational integrators with forces for interconnected systems, extending variational integrator theory.

Findings

Derived a new class of Dirac variational integrators with force preservation properties

Demonstrated improved structure preservation in interconnected systems

Discussed potential applications and future research directions

Abstract

The dynamical motion of mechanical systems possesses underlying geometric structures, and preserving these structures in numerical integration improves the qualitative accuracy and reduces the long-time error of the simulation. For a single mechanical system, structure preservation can be achieved by adopting the variational integrator construction. This construction has been generalized to more complex systems involving forces or constraints as well as to the setting of Dirac mechanics. Variational integrators have recently been applied to interconnected systems in Parks and Leok (2017), which are an important class of practically useful mechanical systems whose description in terms of Dirac structures and Dirac mechanical systems was elucidated in Jacobs and Yoshimura (2014). Since these interconnected systems are modeled as a collection of subsystems with forces of interconnection, we revisit some of the properties of forced variational integrators. In particular, we derive a class of Dirac variational integrators with forces that exhibit preservation properties that are critical when applying variational integrators to the discretization of interconnected Dirac systems. We close with a discussion of ongoing and future research based on these findings.

Tables

Table 1. Table 1. Fourth and second-order forced integrators show the expected convergence rates. The second-order integrator uses Trapezoidal quadrature and second-order Taylor shooting. The fourth-order integrator uses Simpson’s Rule quadrature and fourth-order Taylor shooting. Both methods simulate a damped harmonic oscillator with m = k = 1 𝑚 𝑘 1 m=k=1 and damping coefficient c = 0.01 𝑐 0.01 c=0.01 .

Steps per period	20	40	80	160
$2^{nd}$ order error	0.2275	0.0557	0.0138	0.0035
at $t=5T$
error ratio	–	4.0844	4.0362	3.9429
$4^{th}$ order error	0.6551$\times {10}^{-4}$	0.0516$\times {10}^{-4}$	0.0034$\times {10}^{-4}$	0.0002$\times {10}^{-4}$
at $t=5T$
error ratio	–	12.696	15.176	17

Equations

$$\delta \sum _{k=0}^{N-1}{L}_d(q_k,q_{k+1})=0,$$

$$D_2{L}_d(q_{k-1},q_k)+D_1{L}_d(q_k,q_{k+1})=0.$$

$${\mathbb{F}}^{+}{L}_d(q_k,q_{k+1})$$

$${\mathbb{F}}^{-}{L}_d(q_k,q_{k+1})$$

$$p_{k,k+1}^{-}$$

$$p_{k-1,k}^{+}$$

$$p_{k-1,k}^{+}=p_{k,k+1}^{-}.$$

$$p_{k+1}$$

$$p_k$$

$$L_d^E(q_0,q_1;h)={\int }_0^{h}L(q_{01}(t),{\dot{q}}_{01}(t))dt.$$

$$L_d^E(q_0,q_1;h)=\underset{\begin{array}{c}q\in C^2([0,h],Q)\\ q(0)=q_0,q(h)=q_1\end{array}}{\mathrm{ext}}{\int }_0^{h}L(q(t),\dot{q}(t))dt.$$

$$\delta \sum _{k=0}^{N-1}{L}_d(q_k,q_{k+1})+\sum _{k=0}^{N-1}[f_d^{-}(q_k,q_{k+1})\cdot \delta q_k+f_d^{+}(q_k,q_{k+1})\cdot \delta q_{k+1}]=0,$$

$$D_2{L}_d(q_{k-1},q_k)+f_d^{+}(q_{k-1},q_k)+D_1{L}_d(q_k,q_{k+1})+f_d^{-}(q_k,q_{k+1})=0.$$

$$p_{k+1}={\mathbb{F}}^{f+}{L}_d(q_k,q_{k+1})$$

$$p_k={\mathbb{F}}^{f-}{L}_d(q_k,q_{k+1})$$

$$\frac{\partial L}{\partial q}(q,\dot{q})-\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}(q,\dot{q})\right)+f_L(q,\dot{q})=0,$$

$$L_d^E(q_0,q_1,h)$$

$$f_d^{E+}(q_0,q_1,h)$$

$$f_d^{E-}(q_0,q_1,h)$$

$${\omega }_{d+}^a(q_k,q_{k+1})={\omega }^a(q_k,{\mathcal{R}}_{q_k}^{-1}(q_{k+1})).$$

$${\Delta }_Q^{d+}=\{(q_k,q_{k+1})\in Q\times Q|{\omega }_{d+}^a(q_k,q_{k+1})=0,a=1,\dots ,m\}.$$

$$\delta \sum _{k=0}^{N-1}[L_d(q_k,q_k^{+})+p_{k+1}(q_{k+1}-q_k^{+})]=0,$$

$$0$$

$$q_{k+1}$$

$$p_{k+1}$$

$$p_k$$

$$(X_d^k,{\mathfrak{D}}^{+}{L}_d(q_k,q_k^{+}))\in D_{{\Delta }_Q}^{d+}.$$

$${\omega }_{d-}^a(q_k,q_{k+1})={\omega }^a(q_{k+1},{\mathcal{R}}_{q_{k+1}}^{-1}(q_k))$$

$$\delta \sum _{k=0}^{N-1}[L_d(q_{k+1}^{-},q_{k+1})-p_k(q_k-q_{k+1}^{-})]=0.$$

$$0$$

$$q_k$$

$$p_k$$

$$p_{k+1}$$

$$(X_d^k,{\mathfrak{D}}^{-}{L}_d(q_{k+1}^{-},q_{k+1}))\in D_{{\Delta }_Q}^{d-}.$$

$$||L_d(q(0),q(h),h)-L_d^E(q(0),q(h),h)||\le C_L{h}^{r+1},$$

$$||f_d^{\pm }(q(0),q(h),h)-f_d^{E\pm }(q(0),q(h),h)||\le C_{f\pm }{h}_{f\pm }^{r+1},$$

Full text

1. Introduction

Reframing analytic mechanics as geometric mechanics by considering mechanical systems as evolving on symplectic or Poisson manifolds reveals the depth of geometric structure in the physical world. Examining the range of geometric structures and their associated invariants—momenta, energy, symplectic forms, etc.—relevant to a particular mechanical motion provides a qualitative picture of the motion. Computational geometric mechanics allows these geometric structures to guide the development of numerical methods, leading to improved qualitative behavior of the numerical solutions.

Analytical mechanics and its geometric properties can be derived from variational or almost-variational principles. Variational integrators mimic this process by deriving integrators as the discrete evolution equations implied by a discrete variational principle. This process was first developed for discretizing Hamilton’s principle in the case of conservative, nondegenerate, unconstrained Lagrangian systems. The discrete Hamilton’s principle produces a set of discrete evolution equations known as the Discrete Euler–Lagrange (DEL) equations. Considered as a one-step map, the DEL equations are symplectic, approximately energy conserving, and satisfy a discrete Noether’s theorem. This mirrors the symplecticity, energy conservation, and momentum conservation at the continuous level.

The variational construction has been extended to accommodate forcing, holonomic constraints [7], and nonholonomic constraints [2]. Most recently, the variational construction was applied to Lagrange–Dirac mechanical systems, which is the Lagrangian view of Dirac mechanics [5]. Dirac mechanics generalizes both Lagrangian and Hamiltonian mechanics to a formulation that accommodates forces and constraints as well as degeneracy. Forces, constraints, and degeneracy each change the geometric structure of the system away from straightforward conservation laws, and Dirac geometry and mechanics elucidates the specifics of these changes [9; 10].

Dirac mechanics was also recently shown to be a useful setting for the study of interconnected systems [4]. Interconnected systems consist of multiple mechanical components joined together to operate as a whole. Many practical engineering systems can be conceptualized this way, and often it is more natural to model each component individually than to attempt to model the entire large-scale system monolithically. In [8], we developed the foundations for incorporating the ease of subsystem-level modeling to geometric integration through an interconnected discrete Dirac simulation framework. That framework is based on incorporating the ideas of [4] into the discrete foundation developed in [5].

This paper takes the first step by extending Dirac variational integrators to a formulation that includes discrete forcing. Connected subsystems exert interaction forces on one another, so accommodating forces may prove useful in interconnected applications. We can also view the interconnected modeling process as breaking a large (possibly conservative) system into smaller, more manageable parts. From this point of view, interaction forces are artificial modeling constructs, and we should ensure that they cancel properly in the final integration, so as to ensure that the large monolithic model is equivalent to the model based on interconnecting the component subsystems. This paper addresses that concern by examining notions of equivalence at both the continuous and discrete levels and providing a criterion for defining forced variational or Dirac integrators which are well defined with respect to changes in representation of equivalent systems.

The paper begins with a review of classic, forced, and Dirac variational integrators in Section 2. Section 3 examines order of accuracy and equivalence preservation in forced variational integrators, including several numerical examples. This study informs the later design and implementation of forced Dirac variational integrators. Section 4 describes the extension of Dirac variational integrators to forced Dirac variational integrators, closing with a numerical example. Section 5 summarizes the implications of this work and suggests future research directions.

2. Background: A review of variational integrators

2.1. Classic variational integrators

Consider a conservative, unconstrained, nondegenerate Lagrangian system. Given the Lagrangian $L:TQ\rightarrow\mathbb{R}$, we can construct a discrete Lagrangian $L_{d}:Q\times Q\rightarrow\mathbb{R}$, which can be viewed as approximating the Type I generating function for the symplectic time-$h$ flow map of the Euler–Lagrange vector field. With this, variational integrators are constructed from the discrete Hamilton’s principle, which states that the discrete action sum is stationary,

[TABLE]

for fixed endpoints, $\delta q_{0}=\delta q_{N}=0$. Computing the infinitesimal variation of the discrete action sum, and collecting the terms involving $\delta q_{k}$ yields the discrete Euler–Lagrange equations,

[TABLE]

Here, $D_{1}L_{d}$ and $D_{2}L_{d}$ denote the derivative of the discrete Lagrangian $L_{d}$ with respect its first and second argument, respectively. We view $Q\times Q$ as the discrete analog of the tangent bundle $TQ$ and the discrete Euler–Lagrange equations as defining a one-step map from $Q\times Q$ to itself that determines $(q_{k},q_{k+1})$ from $(q_{k-1},q_{k})$.

One can also define a discrete flow along $T^{*}Q$ by introducing two discrete Legendre transforms, $\mathbb{F}^{\pm}L_{d}:Q\times Q\rightarrow T^{*}Q$,

[TABLE]

These give two possible definitions for the momentum at the point $q_{k}$, which are defined in terms of the data $(q_{k},q_{k+1})$ and $(q_{k-1},q_{k})$,

[TABLE]

Then, the DEL equations enforce a momentum matching condition along the discrete flow, i.e.,

[TABLE]

Thus, the momentum $p_{k}=p_{k-1,k}^{+}=p_{k,k+1}^{-}$ is well-defined along solutions of the DEL equations, and we can define a one-step map in phase space by

[TABLE]

This map is called the discrete Hamiltonian map corresponding to the DEL equations.

The discrete Euler–Lagrange equations and discrete Hamiltonian map each preserve the appropriate symplectic form and satisfy a discrete Noether’s theorem, so that symmetries in the discrete Lagrangian $L_{d}$ result in conservation of the component of the discrete momenta in the direction of the infinitesimal generators of the symmetry. As symplectic integrators, they also approximately conserve energy over exponentially long time scales. These conservation properties hold and ensure that the qualitative behavior of the Lagrangian system is well-approximated for any choice of $L_{d}$.

The choice of $L_{d}$ is guided by the existence of an exact discrete Lagrangian, $L_{d}^{E}$, which yields a discrete Hamiltonian map that samples the exact Hamiltonian flow. For points $q_{0},q_{1}\in Q$, and $h$ sufficiently small, let $q_{01}(t)$ be the unique curve such that $q_{01}(0)=q_{0},\ q_{01}(h)=q_{1}$, and $q_{01}(t)$ satisfies the Euler–Lagrange equations for $L$ on $[0,h]$. Then $L_{d}^{E}$ is given by

[TABLE]

This is related to Jacobi’s solution of the Hamilton–Jacobi equation, and can viewed as evaluating the action integral on the solution of an Euler–Lagrange boundary-value problem, or it can equivalently be described variationally,

[TABLE]

One can approximate the exact discrete Lagrangian $L_{d}^{E}$ by incorporating a wide variety of standard numerical techniques in order to obtain a computable discrete Lagrangian $L_{d}$ [6; 3]. By the variational error analysis described in [7], it can be shown that if the discrete Lagrangian $L_{d}(q_{0},q_{1};h)=L_{d}^{E}(q_{0},q_{1};h)+\mathcal{O}(h^{r+1})$, then the associated discrete Hamiltonian map is an order $r$ approximation of the exact flow map.

2.2. Forced variational integrators

Given a continuous force $f_{L}:TQ\rightarrow T^{*}Q$ that is a fiber-preserving map over the identity, we introduce discrete forces $f_{d}^{\pm}:Q\times Q\rightarrow T^{*}Q$ that approximate the virtual work of the forces. Then, the forced discrete Euler–Lagrange equations are derived from the discrete Lagrange–d’Alembert principle,

[TABLE]

where as before, the endpoints are fixed, $\delta q_{0}=\delta q_{N}=0$. The forced discrete Euler–Lagrange equations are given by

[TABLE]

The discrete Legendre transforms and discrete Hamiltonian map now incorporate the contribution of the discrete forces,

[TABLE]

With these momentum definitions, the forced DEL equations can again be viewed as a momentum matching condition.

These equations reduce to the DEL equations in the absence of forcing, and they satisfy a forced discrete Noether’s theorem for symmetries of the discrete Lagrangian where the discrete forces do no work in the direction of the infinitesimal generators of the symmetry. Forced variational integrators exhibit better energy behavior than non-geometric integrators in practice, in the sense that the discrete energy evolution better reflects the exact energy evolution of the system. This lacks rigorous explanation since the forced equations are no longer symplectic nor energy preserving at either the continuous or the discrete level.

Once again, exact discrete quantities guide both the error analysis and the practical choice of $L_{d}$ and $f_{d}^{\pm}$. For points $q_{0},q_{1}\in Q$, define $q_{01}(t)$ as above, except that we now require $q_{01}(t)$ to satisfy the forced Euler–Lagrange equations,

[TABLE]

for $L$ and $f_{L}$ on $[0,h]$. Then, $L_{d}^{E}$ and $f_{d}^{E\pm}$ are given by

[TABLE]

We obtain computable $L_{d}$ and $f_{d}^{\pm}$ by approximating the curve $q_{01}(t)$ that satisfies the forced Euler–Lagrange boundary-value problem and the integrals that arise in the definition of $L_{d}^{E},f_{d}^{E\pm}$. Again, the orders of these approximations determine the order of accuracy of the discrete Hamiltonian map. This fact is stated without proof in [7], and we provide an explicit proof in the following section.

2.3. Dirac variational integrators

Dirac variational integrators were developed in [5] and generalize the variational integrator construction to the case of Dirac mechanics. In [5], the authors develop both the variational theory of discrete Dirac mechanics and explicit discrete analogues of Dirac structures. The authors arrive at two formulations, (+) and ($-$)-discrete Dirac mechanics, stemming from a choice of generating function used to define the discrete analogue of the symplectic flat map, $\Omega_{d\pm}^{\flat}$.

Discretization begins from a continuous problem described by a Lagrangian function $L:TQ\rightarrow\mathbb{R}$ and a continuous constraint distribution $\Delta_{Q}\subset TQ$. The constraint distribution is determined by its annihilator, which we write as $\Delta_{Q}^{\circ}=\text{span}\{\omega^{a}\}_{a=1}^{m}\subset T^{*}Q$. In each formulation of discrete Dirac mechanics, we construct a discrete Lagrangian $L_{d}$ and discrete annihilating one-forms $\omega_{d\pm}^{a}$ based on the continuous problem and a choice of retraction $\mathcal{R}:TQ\to Q$.

2.3.1. (+)-discrete Dirac mechanics

In (+)-discrete Dirac mechanics, we define the discrete one-forms as

[TABLE]

These then define a discrete distribution as follows

[TABLE]

The (+)-discrete Lagrange–d’Alembert–Pontryagin principle is

[TABLE]

with variations that vanish at the endpoints, $\delta q_{0}=\delta q_{N}=0$, and discrete constraints $(q_{k},q_{k+1})\in\Delta_{Q}^{d+}$. We also impose a constraint on the variations $\delta q_{k}\in\Delta_{Q}(q_{k})$ after computing variations inside the sum. The variable $q_{k}^{+}$ serves as the discrete analog to the introduction of $v$ in the continuous principle. This process produces the (+)-discrete Lagrange–Dirac equations of motion,

[TABLE]

where $a=1,\dots,m$, and the last equation uses the Einstein summation convention. These equations simplify to the DEL equations in the unconstrained case, and they recover the nonholonomic integrators of [2]. They can also be expressed in terms of discrete Dirac structures as

[TABLE]

where $X_{d}^{k}=((q_{k},p_{k}),(q_{k+1},p_{k+1}))$ is the discrete vector field, $\mathfrak{D}^{+}L_{d}$ is the (+)-discrete Dirac differential, and $D_{\Delta_{Q}}^{d+}$ is the (+)-discrete Dirac structure induced by the continuous distribution $\Delta_{Q}$. The discrete symplectic flat map, $\Omega_{d+}^{\flat}$, contributes to the definition of $D_{\Delta_{Q}}^{d+}$. We show (3) now to highlight its similarity to the continuous expression for constrained Dirac mechanics, $(X,\mathfrak{D}L)\in D_{\Delta_{Q}}$. We provide explicit descriptions of the discrete objects as needed in the development of forced Dirac integrators below.

2.3.2. ($-$)-discrete Dirac mechanics

This formulation defines the discrete one-forms as

[TABLE]

and uses the ($-$)-discrete Lagrange-d’Alembert-Pontryagin principle

[TABLE]

We impose the same constraints on the variations of the discrete principle as in the (+) case, using $\Delta_{Q}^{d-}$ defined from $\{\omega_{d-}^{a}\}_{a=1}^{m}$ instead of $\Delta_{Q}^{d+}$. The variable $q_{k}^{-}$ now plays the role of the discrete velocity. This yields the ($-$)-discrete Lagrange–Dirac equations

[TABLE]

Equivalently,

[TABLE]

Again, we avoid defining each discrete object until necessary, and the interested reader is referred to [5] for the details.

3. Forced Integrators revisited

In [7], the authors assert without proof the existence of a theorem relating the order of a forced variational integrator to the order of approximation of $L_{d},f_{d}^{\pm}$ to $L_{d}^{E},f_{d}^{E\pm}$. Here, we provide an explicit proof of the most practical direction of that theorem and a recipe for constructing forced integrators of known order. We then analyze the conditions under which this process yields a well-defined, equivalence-preserving integrator and discuss the implications of equivalence-preservation in interconnected applications. The section closes with several numerical examples illustrating both our construction process and the consequences of equivalence-preservation.

3.1. Determining the order of a forced variational integrator

Theorem 1 below explicitly proves that the orders of approximation of $L_{d}$ and $f_{d}^{\pm}$ to $L_{d}^{E}$ and $f_{d}^{E\pm}$ determine the order of accuracy of the discrete Hamiltonian map they define. This gives the most useful direction of the order theorem mentioned in [7], since the order of $L_{d}$ and $f_{d}^{\pm}$ are easier to calculate than the order of the integrator.

We need a few preliminary definitions before we state the theorem. First, we explicitly define what we mean by the order of $L_{d},f_{d}^{\pm}$ and $\mathbb{F}^{f\pm}L_{d}$, using the same definitions as in [7]. Thus, a given $L_{d}$ is of order $r$ if there exist constants $C_{L}>0$, $h_{L}>0$ and an open subset $U_{L}\subset TQ$ with compact closure such that

[TABLE]

for all solutions $q(t)$ of the forced Euler–Lagrange equations with initial conditions satisfying $(q(0),\dot{q}(0))\in U_{L}$, and for all $h\leq h_{L}$. Similarly, $f_{d}^{\pm}$ are of order $r$ if there exist constants $C_{f\pm}>0$, $h_{f\pm}>0$ and open subsets $U_{f\pm}\subset TQ$ with compact closure such that

[TABLE]

for all solutions $q(t)$ of the forced Euler–Lagrange equations with initial conditions satisfying $(q(0),\dot{q}(0))\in U_{f}^{\pm}$, and for all $h\leq h_{f}^{\pm}$. We will say that $L_{d},f_{d}^{\pm}$ are simultaneously of order $r$ if there exists $U\subset U_{L}\cap U_{f}^{+}\cap U_{f}^{-}$ such that $U$ is a nontrivial open set with compact closure.

For the discrete Legendre transforms to be of order $r$, we require constants $C^{\pm}>0,h^{\pm}>0$ and an open set $U^{\pm}\subset TQ$ such that

[TABLE]

for all solutions $q(t)$ of the forced Euler–Lagrange equations with initial conditions satisfying $(q(0),\dot{q}(0))\in U^{\pm}$, and for all $h\leq h^{\pm}$.

Theorem 1**.**

Consider a hyperregular Lagrangian $L$ with force $f_{L}$, corresponding Hamiltonian $H$, and corresponding Hamiltonian force $f_{H}$. If $L_{d}$ and $f_{d}^{\pm}$ are simultaneously of order $r$, then

1. a.

the forced discrete Legendre transforms $\mathbb{F}^{f\pm}L_{d}$ are of order $r$. 2. b.

the discrete Hamiltonian map is of order $r$.

Proof.

a. Assume that $L_{d}$ and $f_{d}^{\pm}$ are simultaneously of order $r$. Then $L_{d}$ of order $r$ implies existences of a function $e$ such that

[TABLE]

and $||e(q(0),q(h),h)||\leq C_{L}$ on $U$. Similarly, $f_{d}^{+}$ of order $r$ implies the existence of a function $e^{+}$ such that

[TABLE]

and $||e^{+}(q(0),q(h),h)||\leq C_{f+}$ on $U$.

Taking derivatives of (4) with respect to $q(h)$ gives

[TABLE]

We assume both $L_{d}$ and $L_{d}^{E}$ have continuous derivatives in $q(h)$, so, by the definition of $e$ in (4), $D_{2}e$ is continuous and bounded on the compact set $\operatorname*{cl}(U)$. Combining this with (5), we have

[TABLE]

From the above arguments, $D_{2}e(q(0),q(h),h)+e^{+}(q(0),q(h),h)$ is bounded on the set $\operatorname*{cl}(U)$. Thus, $\mathbb{F}^{f+}L_{d}$ is of order $r$. Taking derivatives of $L_{d}$ with respect to $q(0)$ and using a similar calculation with $f_{d}^{-}$ shows that $\mathbb{F}^{f-}L_{d}$ is of order $r$.

b. Let $F^{f}_{L_{d}}$ denote the integrator defined by the forced discrete Euler Lagrange equations and $\tilde{F}^{f}_{L_{d}}$ denote the discrete Hamiltonian integrator defined by the forced discrete Legendre transforms in (2). Then, from the definitions of $F_{L_{d}},\tilde{F}_{L_{d}}$, and $\mathbb{F}^{f\pm}L_{d}$, we have the commutative diagram (6), which is the forced equivalent of diagram (1.5.3) in [7].

[TABLE]

Thus, we can express the forced discrete Hamiltonian map as $\tilde{F}_{L_{d}}^{f}=\mathbb{F}^{f+}L_{d}\circ(\mathbb{F}^{f-}L_{d})^{-1}$, directly analogous to the unforced discrete Hamiltonian map, $\tilde{F}_{L_{d}}=\mathbb{F}^{+}L_{d}\circ(\mathbb{F}^{-}L_{d})^{-1}$. As such, the proof of b. from a. given in [7] carries over exactly, and we choose not to reproduce it here. ∎

The next theorem provides a recipe for constructing forced triples $(L_{d},f_{d}^{\pm})$ of known order.

Theorem 2**.**

Consider a $p^{th}$ order accurate boundary-value solution method and a $q^{th}$ order accurate quadrature formula. Assume both the Lagrangian and the forcing function are Lipschitz continuous in both variables. Use the boundary-value solution method to obtain approximations $(q^{i},v^{i})\approx(q(c_{i}h),v(c_{i}h))$, at the quadrature nodes $c_{i}$ of the solution $(q(t),v(t))$ of the forced Euler–Lagrange boundary-value problem. Then the associated discrete Lagrangian given by

[TABLE]

and the discrete forces given by

[TABLE]

all have order of accuracy $\min(p+1,q)$.

Proof.

The result follows from a few straightforward calculations. We begin with $L_{d}$.

[TABLE]

For $f_{d}^{+}$,

[TABLE]

A similar calculation shows that

[TABLE]

∎

Taken together, Theorems 1 and 2 show that we can construct forced variational integrators of known order from a choice of a quadrature rule and a boundary-value solution method. In fact, we could choose up to three different quadrature rules and boundary-value solution methods to define $L_{d}$, $f_{d}^{+}$, and $f_{d}^{-}$. The proofs above would still apply, and the resulting integrator would have order min($p_{1}+1,p_{2}+1,p_{3}+1,q_{1},q_{2},q_{3})$ for $p_{i}$ the orders of the boundary-value solutions and $q_{i}$ the orders of the quadrature rules. However, this produces integrators with unpredictable results, as discussed in the next section.

3.2. Notions of equivalence and equivalence-preservation

A forced Lagrangian system has equations of motion

[TABLE]

Since the equations of motion are defined by the combination of $L$ and $f_{L}$, it is possible for pairs $(L^{1},f_{L}^{1})$ and $(L^{2},f_{L}^{2})$ with $L^{1}\neq L^{2}$ and $f_{L}^{1}\neq f_{L}^{2}$ to define the same equations of motion. We refer to this as equivalence of the pairs $(L^{1},f_{L}^{1})$ and $(L^{2},f_{L}^{2})$. Most numerical methods simulate mechanics by numerically integrating the differential equations of motion, so they produce the same numerical approximation whether the motion was originally described using $(L^{1},f_{L}^{1})$ or $(L^{2},f_{L}^{2})$.

Variational integrators simulate mechanics using the discrete equations of motion defined by a particular choice of $L_{d}$ and $f_{d}^{\pm}$. Thus, a variational integrator is defined by the choice of discretizations $L\mapsto L_{d}$ and $f_{L}\mapsto f_{d}^{\pm}$. It is therefore natural to ask which discretization schemes for the discrete Lagrangian and discrete forces when applied to equivalent representations of forced Euler–Lagrange systems lead to equivalent discrete equations of motion. This reflects whether or not the resulting variational integrator is well-defined on the equivalence class of representations of forced Euler–Lagrange systems.

More explicitly, given two equivalent representations $(L^{1},f_{L}^{1})$ and $(L^{2},f_{L}^{2})$ of a forced Euler–Lagrange system, the resulting solution trajectories are the same. Since a variational integrator for a forced Euler–Lagrange system aims to approximate that solution, it is desirable to consider well-defined variational integrators that produce the same approximation irrespective of which of the two equivalent representations $(L^{1},f_{L}^{1})$ or $(L^{2},f_{L}^{2})$ it is applied to. Well-definedness in this sense will be important in interconnected applications where we intentionally alter the forced representation of a system to view it as a collection of interacting subsystems. In particular, if we considered a system with many components in terms of the constituent free-body diagrams, it is essential when combining the free-body diagrams for the internal forces to cancel out in order to recover the original system.

We discuss the implications of equivalence at the continuous and discrete levels. Then, we provide a simple criterion for constructing well-defined forced variational integrators.

3.2.1. Continuous equivalence

Suppose we have two canonical Lagrangians, $L^{i}=\frac{1}{2}\dot{q}^{T}M^{i}\dot{q}-V^{i}(q)$. Then, equality of the equations of motion implies

[TABLE]

so that

[TABLE]

Comparing variables on each side, we conclude that $M^{2}=M^{1}$ and $f_{L}^{1}(q,v)-f_{L}^{2}(q,v)=\nabla V^{1}(q)-\nabla V^{2}(q)$.

More generally, assume two Lagrangians of the form $L^{i}(q,v)=K(q,v)-V^{i}(q)=\frac{1}{2}g(v,v)-V^{i}(q)$ for some metric tensor $g$. Then, continuous equivalence again implies $f_{L}^{1}(q,v)-f_{L}^{2}(q,v)=\nabla V^{2}(q)-\nabla V^{1}(q)$.

3.2.2. Notions of discrete equivalence.

In [7], the authors define strongly equivalent discrete Lagrangians to be those that generate equivalent discrete Hamiltonian maps and weakly equivalent discrete Lagrangians to be those that generate equivalent discrete Lagrangian maps. That is to say that strongly equivalent discrete Lagrangians will generate the same discrete sequence $\{(q_{k},p_{k})\}_{k=0}^{N}\subset T^{*}Q$, whereas weakly equivalent discrete Lagrangians will only generate the same sequence $\{q_{k}\}_{k=0}^{N}\subset Q$.

We can define strong and weak equivalence of discrete triples $(L_{d}^{i},f_{d}^{i\pm})$ in the same way. From diagram (6), we see that we have strong equivalence if and only if we have equivalence of the forced discrete Legendre transforms, $\mathbb{F}^{f\pm}L_{d}$, just as in the unforced case. Then, Theorem 1 still holds for triples which are strongly equivalent to a triple meeting the assumptions. Moreover, Dirac variational integrators directly generalize the discrete Hamiltonian implementation of forced variational integrators. Thus, we focus on strong equivalence.

Strongly equivalent discrete triples $(L_{d}^{1},f_{d}^{1\pm})$ and $(L_{d}^{2},f_{d}^{2\pm})$ have equivalent forced discrete Legendre transforms, i.e. $\mathbb{F}^{f\pm}L_{d}^{1}=\mathbb{F}^{f\pm}L_{d}^{2}$. This implies that

[TABLE]

and

[TABLE]

So each difference in discrete forcing in two strongly equivalent triples must be integrable in at least one of the position variables. If we think of $f_{d}^{\pm}$ as right and left discrete forces, then each discrete force difference must be integrable in its base-point variable, mirroring the continuous conclusion that $f_{L}^{1}(q,v)-f_{L}^{2}(q,v)=\nabla V^{2}(q)-\nabla V^{1}(q)$.

3.2.3. Preserving continuous equivalence

Theorems 1 and 2 provide a means of constructing forced variational integrators of known order by choosing quadrature rules and boundary-value solution methods. It is tempting to try to optimize the overall discretization by tailoring the discretization of $L_{d}$, $f_{d}^{+}$, and $f_{d}^{-}$ individually, but the resulting integrator is no longer well-defined with respect to the equivalence relation defined above. The simplest way to generate an integrator that preserves continuous equivalence by our method is to choose the same quadrature rule and boundary-value solution method for all three discrete quantities.

In this case, we have

[TABLE]

and

[TABLE]

Suppose $(L^{1},f_{L}^{1})$ and $(L^{2},f_{L}^{2})$ are equivalent, so

[TABLE]

Constructing $(L_{d}^{1},f_{d}^{1\pm})$ and $(L_{d}^{2},f_{d}^{2\pm})$ from $(L^{1},f_{L}^{1})$ and $(L^{2},f_{L}^{2})$, we then have

[TABLE]

Similarly, $f_{d}^{2-}(q_{0},q_{1};h)-f_{d}^{1-}(q_{0},q_{1};h)=D_{1}L_{d}^{1}(q_{0},q_{1};h)-D_{1}L_{d}^{2}(q_{0},q_{1};h)$. Thus, $(L_{d}^{1},f_{d}^{1\pm})$ and $(L_{d}^{2},f_{d}^{2\pm})$ are strongly equivalent.

3.3. Numerical examples

We construct our implementations by choosing a quadrature rule and a one-step map, using shooting to solve the boundary-value problem, which is analogous to the shooting-based variational integrator proposed in [6]. Quadrature of order $q$ and a one-step map of order $p$ produce an integrator of order min($p,q$) since the shooting solution is order $p-1$.

Table 1 shows this construction converges as predicted. The second-order integrator uses trapezoidal quadrature and a second-order Taylor’s method shooting solution. The fourth-order integrator uses Simpson’s Rule quadrature and a fourth-order Taylor’s method shooting solution. Both were run on a damped harmonic oscillator with mass $m=1$, spring constant $k=1$, and damping coefficient $c=0.01$. The table shows the effect of doubling the number of time steps per period $T$ on the error after five periods. In both cases, the error decreases at the predicted rate.

Figure 1 illustrates the unpredictability of the results of non-equivalence-preserving discretizations. Two implementations use second-order Taylor’s method shooting as the boundary-value solution method for both $L_{d}$ and $f_{d}^{\pm}$ with a Trapezoidal and/or Midpoint quadrature. The equivalence-preserving discretization uses Trapezoidal quadrature on both $L_{d}$ and $f_{d}^{\pm}$. The non-equivalence-preserving implementation uses a Trapezoidal quadrature on $L_{d}$ and Midpoint quadrature on $f_{d}^{\pm}$. Both integrators were run on a conservative harmonic oscillator eleven times, with the potential force increasingly represented as an external force, i.e.

[TABLE]

for $\alpha=0,0.1,0.2\dots,1$. The step size was $h=0.05$ for all runs. All values of $\alpha$ produce the same continuous equations of motion. The equivalence-preserving discretization produces the same solution regardless of $\alpha$. The non-equivalence-preserving discretization produces $\alpha$-dependent solutions which veer away from the unforced representation’s solution as $\alpha$ increase. We also compared a Midpoint-Midpoint construction to a Midpoint-Trazezoidal construction in the same way. The results lie directly atop the Trap-Trap vs. Trap-Mid results shown. This unpredictability of results makes non-equivalence-preserving integrators unsuited for applications like interconnected systems in which a system’s forced representation will purposely be modified by modeling the component subsystems individually, with internal forces to account for the interaction with other subsystems.

Figure 2 shows a different example where failure to preserve equivalence produces wildly unpredictable results. We constructed a shooting based variational integrator with fourth-order Runge–Kutta as the underlying one-step method, and used the trapezoidal rule on $L_{d}$ and $f_{d}^{\pm}$ to construct a second-order equivalence-preserving integrator, and Simpson’s rule on $L_{d}$ and trapezoidal rule on $f_{d}^{\pm}$ to construct a second-order non-equivalence-preserving integrator. Both of these methods were applied to simulate a system whose equations of motion simplify to an unforced harmonic oscillator. We inserted an artificial potential of $100q^{5}$, then canceled that potential force with a forcing function. As can be seen from the figure, the non-equivalence-preserving formulation fails to capture the motion and exhibits a qualitatively different behavior. The example is artificial, but it highlights the need to preserve equivalence in developing a well-defined integrator.

4. Dirac variational integrators with forces

We now lay out the generalization of Dirac variational integrators to include forces. We begin from the variational perspective, where we can directly apply the ideas of forced variational integrators. We then formulate the resulting integrators from the perspective of discrete Dirac structures and finish with a basic numerical implementation.

4.1. Variational formulation

4.1.1. (+)-discrete Dirac mechanics

In [5], the discrete variational principle for a Dirac system without external forces is the (+)-discrete Lagrange-d’Alembert-Pontryagin principle, given by

[TABLE]

where the variations vanish at the endpoints and are constrained to lie on a constraint distribution, and pairs of consecutive points in the discrete solution lie in a discrete constraint distribution. We define the forced principle in direct analogy to (1) as

[TABLE]

where the second sum represents the virtual work associated with the external forces. We have chosen the arguments of $f_{d}^{\pm}$ to match those of $L_{d}$, and we constrain the variations in the same way as [5]. Thus, $\delta p_{k}$ and $\delta q_{k}^{+}$ are arbitrary and $\delta q_{0}=\delta q_{n}=0$. We impose $\delta q_{k}\in\Delta_{Q}(q_{k})$ on the remaining $\delta q_{k}$ after taking variations inside the sum, and we insist that $(q_{k},q_{k+1})\in\Delta_{Q}^{d+}$. For the forced principle (7), these constraints yield

[TABLE]

These equations recover the (+)-discrete Lagrange–Dirac equations of [5] when $f_{d}^{\pm}=0$ and the forced discrete Euler–Lagrange equations [7] in the absence of constraints. When both the forces and constraints are zero, we recover the classic discrete Euler–Lagrange equations.

4.1.2. ($-$)-discrete Dirac mechanics

In [5], ($-$)-discrete Dirac mechanics derives from the ($-$)-discrete Lagrange-d’Alembert principle,

[TABLE]

with constraints on the variations, and the pairs of sequential points as before. We add forces to this principle as

[TABLE]

again choosing the arguments of $f_{d}^{\pm}$ to match those of $L_{d}$ and constraining the variations as in [5]. In this case, then, $\delta p_{k}$ and $\delta q_{k}^{-}$ are arbitrary and $\delta q_{0}=\delta q_{N}=0$. We impose $\delta q_{k}\in\Delta_{Q}(q_{k})$ on the remaining $\delta q_{k}$ after taking variations inside the sums, and we insist that $(q_{k},q_{k+1})\in\Delta_{Q}^{d-}$. With these conditions, computing variations of the principle (8) yields

[TABLE]

In the unforced case, these equations recover the ($-$)-discrete Dirac equations of [5]. In the unconstrained case, we recover the forced discrete Euler–Lagrange equations of [7]. In the absence of both we again recover the classic discrete Euler–Lagrange equations.

4.2. Discrete Dirac structure formulation

The dynamics of an unforced, continuous Lagrange–Dirac dynamical system can be expressed in terms of Dirac structures as $(X,\mathfrak{D}L)\in D$, where $X$ is the partial vector field of the motion, $\mathfrak{D}L$ the Dirac differential of the Lagrangian $L$ and $D$ a Dirac structure on $T^{*}Q$. To accomodate forces in the continuous case, the force $F:TQ\to T^{*}Q$ is lifted to a map $\tilde{F}:TQ\to T^{*}T^{*}Q$ such that

[TABLE]

Locally, $\tilde{F}$ is given by $\tilde{F}(q,v)=(q,p,F(q,v),0)$. The forced Lagrange–Dirac dynamical system is then given by

[TABLE]

Without forces, we can express (+) and ($-$)-discrete Dirac mechanics in terms of discrete Dirac structures as

[TABLE]

and

[TABLE]

which amount to the conditions

[TABLE]

To introduce forces into these expressions in a manner that is consistent with the discrete equations of motion that are derived variationally, we require

[TABLE]

To this end, define $\tilde{F}_{d+}$ and $\tilde{F}_{d-}$ as

[TABLE]

Then, we can write conditions (10) as

[TABLE]

mimicking the continuous expression (9). In (11) we regard $\mathfrak{D}^{\pm}L_{d}-\tilde{F}_{d\pm}$ as an abstract expression and impose base-point matching after the subtraction.

4.3. Numerical Example

We close with a simple implementation of the (+) forced discrete Dirac equations. We simulate a basic RLC resonator with a single resistor, inductor, and capacitor. We assume some existing voltage in the system but no replenishing voltage source.

[TABLE]

The system is described as a Lagrange–Dirac system in terms of the charges and currents in each component, $(q^{C},q^{L},q^{R},i^{C},i^{L},i^{R})$. Let $L,R,$ and $C$ denote the inductance, resistance, and capacitance of the components. Then the system has Lagrangian function

[TABLE]

and force

[TABLE]

Kirchhoff’s current law imposes the following constraints,

[TABLE]

Thus, our constraints can be described in terms of the annihilator distribution $\Delta_{Q}^{\circ}=span\{\omega^{1},\omega^{2}\}$ for $\omega^{1}=dq^{L}-dq^{R}$ and $\omega^{2}=dq^{R}-dq^{C}$. We follow the discretization in [5] to construct $\omega_{d+}^{1,2}$ as $\omega_{d+}^{a}(q_{k},q_{k+1})=\omega^{a}(q_{k},\mathcal{R}_{q_{k}}^{-1}(q_{k+1}))$ for $\mathcal{R}_{q}(v)=q+vh$.

Existence of exact discrete quantities in the Dirac setting has not yet been established, so we cannot turn to them to guide the discretization for $L_{d}$ and $f_{d}^{\pm}$. Establishing such quantities is obviously desirable and a topic for future work. For now, we choose a simple $L_{d}$ and $f_{d}^{\pm}$ shown to work well in both the forced, non-Dirac setting [7] as well as the nonholonomic integrator setting [2]. Namely, we choose $L_{d}$ and $f_{d}^{\pm}$ to be

[TABLE]

and

[TABLE]

[7] introduces these forces as the natural complement to $L_{d}^{1/2}$. The discretizations $L_{d}^{1/2}$ and $f_{d}^{1/2\pm}$ both correspond to choosing midpoint quadrature and a linear boundary-value solution to approximate $L_{d}^{E}$ and $f_{d}^{E\pm}$ in the forced, non-Dirac setting. Thus, the set $(L_{d}^{1/2},f_{d}^{1/2\pm})$ is natural and equivalence-preserving by our earlier analysis as well.

Figure 3 shows that the simulated charge in the capacitor over time (dotted line) is almost identical to the exact solution (solid line) using a time-step of $h=0.05$ for 1000 time-steps.

5. Conclusions and future work

In this work we have successfully integrated the presence of external forcing into the discrete Dirac framework originally presented in [5] so that our forced discrete Dirac integrators recover the forced DEL integrators of [7]. We have also studied the effects of shifting forced representations on the accuracy of numerical simulations and presented a straightforward method for constructing integrators that avoid these spurious effects.

The work is motivated by the applications to the interconnection of discrete Dirac mechanics. Continuous interconnections can be equivalently represented as either constraints or interaction forces acting between subsystems. One would hope to obtain a similar equivalence at the discrete level. In recent work [8], we have derived discrete interconnections from the point of view of constraints. However, we do not yet have an equivalent representation in terms of discrete interaction forces. There is a lack of equivalence between forced and constrained representations throughout the discrete literature, though recent work has rectified this issue on vectors spaces through Hamel’s formalism [1].

6. Acknowledgements

HP has been supported by the NSF Graduate Research Fellowship grant number DGE-1144086. ML has been supported in part by NSF under grants DMS-1010687, CMMI-1029445, DMS-1065972, CMMI-1334759, DMS-1411792, DMS-1345013.