Prolongation-Collocation Variational Integrators

0.24 MB·English

by  Melvin Leok

PROLONGATION-COLLOCATION VARIATIONAL INTEGRATORS
MELVIN LEOK AND TATIANASHINGEL
Abstract. We introduce a novel technique for constructing higher-order variational integrators
forHamiltoniansystemsofODEs. Inparticular,weareconcernedwithgeneratinggloballysmooth
1 approximationsto solutions ofaHamiltonian system. Ourconstruction of thediscreteLagrangian
1 adopts Hermite interpolation polynomials and the Euler–Maclaurin quadrature formula, and in-
0
volves applying collocation to the Euler–Lagrange equation and its prolongation. Considerable
2
attentionisdevotedtotheorderanalysisoftheresultingvariationalintegratorsintermsofapprox-

n imation properties of the Hermite polynomials and quadrature errors. A performance comparison
a is presented on a selection of these integrators.
J

1
1

] 1. Introduction
A
N One of the major themes in geometric integration is symplectic methods for solving Hamiltonian
systems of ordinary differential equations. Viewed as maps, such methods preserve the symplectic
.
h
two-form underlying the dynamical evolution of the system. Variational integrators are an impor-
t
a tant class of symplectic integrators, which arise from discretizing Hamilton’s principle. We refer
m
the reader to [11] for a detailed discussion of the background theory of such methods. Variational
[
integrators are automatically symplectic and momentum preserving. Moreover, they exhibit good

1 energy behavior for exponentially long times.
v The construction of variational integrators combines the techniques from approximation theory
5
and numerical quadrature and is related to the Galerkin approach of converting a differential
9
9 operator equation into a discrete system. Galerkin methods are semi-analytic in the sense that the
1 discrete solution is described by an element of a finite-dimensional function space, which provides
.
1 an analytic expression for the numerical solution. Our goal is to develop variational integrators
0 that would lead to globally smooth approximations of the solution. To pursue this goal, we adopt
1
the space of piecewise Hermite polynomials in the Galerkin construction. It is a well-known result
1
: in approximation theory that Hermite polynomials allow for higher-order approximation of smooth
v
functions at relatively low cost. Recall that Lagrange polynomials are a common choice in the
i
X construction of higher-order variational integrators. However, the computation of an interpolating
r polynomial in Lagrange form becomes unstable when the degree of the polynomial is high. Unlike
a
Lagrange polynomials, higher-degree Hermite polynomials produce accurate and computationally
stable results.
Furthermore, the Galerkin construction based on piecewise Hermite polynomial interpolation
allows us to adequately address the question of variational order error analysis. The order error
analysis of variational integrators relies on determining the order with which a discrete Lagrangian
L :Q×Q → R approximates the exact discrete Lagrangian,
d
h
(1) LE(q ,q ;h) = L(q (t),q˙ (t))dt,
d 0 1 01 01
Z0
where q (t) is a solution curve of the Euler–Lagrange equation that satisfies the boundary condi-
01
tions q (0) = q , q (h) = q . The standard way of performing the error analysis is by comparing
01 0 01 1
the Taylor expansions of the exact discrete Lagrangian and the discrete Lagrangian. Consequently,
the result critically depends on the extent to which the discrete trajectory is able to approximate
the higher-derivatives of the exact Euler–Lagrange solution curve. To enable a robust variational
1
2 MELVINLEOKANDTATIANASHINGEL
order error analysis for the Galerkin variational integrators, we propose a novel application of the
collocation approach in the setting of discrete Lagrangian mechanics, which involves prolongations
of the Euler–Lagrange vector field. The proposed approach has the advantage that one is able to
prove optimal rates of convergence of the associated variational integrators as long as sufficiently
accurate quadrature formulas are used. The proof crucially depends on the preliminary estimates
on the higher-derivative approximation properties of prolongation-collocation curves.
1.1. Outline of the Paper. We present a brief review of discrete variational mechanics and
variational integrators in Section 2, and the Euler–Maclaurin quadrature formula in Section 3. In
Section 4, we introduce prolongation-collocation variational integrators, and perform variational
order error analysis in Section 5. In Section 6, we present a few numerical examples, and present
some conclusions and future directions in Section 7.
2. Variational Integrators
Let Q be the configuration manifold of a mechanical system, with generalized coordinates q.
Consider the Lagrangian L : TQ → R, where TQ is the tangent bundle of the configuration space
Q. The tangent bundle has local coordinates (q,v). Further, let C(Q) = C([0,T],Q) denote the
space of smooth trajectories q : [0,T] → Q in the configuration manifold Q. The action integral
S : C(Q) → R is defined as
T
S(q)= L(q(t),q˙(t))dt.
Z0
The variational principle known as Hamilton’s Principle states that
δS = 0, for δq(0) = δq(T) = 0,
which yields the Euler–Lagrange equations
∂L d ∂L
(2) (q,q˙)− (q,q˙) = 0.
∂q dt ∂q˙
(cid:18) (cid:19)
Itispossibletorewritetheequations(2)intermsofthegeneralized coordinates andmomenta(q,p)
on the cotangent bundle T∗Q (phase space). For this, we introduce the Legendre transformation
FL :TQ → T∗Q, defined by
∂L
FL: (q,q˙) 7→ q, .
∂q˙
(cid:18) (cid:19)
The Hamiltonian H : T∗Q → R is given by
H(q,p) = p·q˙−L(q,q˙)| .
p=∂L
∂q˙
One can show that equations (2) are equivalent to Hamilton’s equations (Theorem 1.3, p. 182 in
[6]),
∂H ∂H
(3) p˙ = − (p,q), q˙ = (p,q).
∂q ∂p
In the case of variational integrators, instead of discretizing the Euler–Lagrange equations (2),
one discretizes Hamilton’s principle. That is, one discretizes the action by introducing a discrete
Lagrangian and replaces the action integral by an action sum, and applies the discrete version of
the Hamilton’s variational principle. The major advantage of this approach is that the resulting
numerical algorithm automatically preserves the symplectic structure of the underlying dynamical
system.
PROLONGATION-COLLOCATION VARIATIONAL INTEGRATORS 3
The discrete Lagrangian L (q ,q ,h) is thought of as an approximation of the action integral
d 0 1
along the curve segment between the points q ≈ q(0) and q ≈ q(h). Formally, this can be
0 1
expressed as
h
(4) L (q ,q ,h) ≈ L(q(t),q˙(t))dt.
d 0 1
Z0
We will neglect the h-dependence and simply write L (q ,q ) when it is not essential for the
d 0 1
exposition of thematerial. Given thediscrete sequenceof times {t = hk|k =0,...,N}, h= T/N,
k
adiscrete curvein Q is denoted by {q }N , whereq ≈ q(t ). Thediscreteaction sumis afunction
k k=0 k k
that maps the discrete trajectories {q }N to R, and is given by
k k=0
N−1
S ({q }N ) = L (q ,q ).
d k k=0 d k k+1
k=0
X
The discrete Hamilton’s principle requires the discrete action to be stationary with respect to
variations vanishing at k = 0 and k = N. From this, one derives the discrete version of the
Euler–Lagrange equations, which are known as the discrete Euler–Lagrange equations,
(5) D L (q ,q )+D L (q ,q ) = 0,
2 d k−1 k 1 d k k+1
where k = 1,2,...,N −1. These equations implicitly define the one-step discrete Lagrangian map
F : Q×Q → Q×Q. The discrete Legendre transforms F±L : Q×Q → R are defined by
Ld d
F−L :(q ,q ) 7→ (q ,−D L (q ,q )),
d 0 1 0 1 d 0 1
F+L :(q ,q ) 7→ (q ,D L (q ,q )).
d 0 1 1 2 d 0 1
Pushing the discrete Lagrangian map F forward to T∗Q with the discrete Legendre transforms
Ld
gives the discrete Hamiltonian map F˜ : T∗Q → T∗Q by F˜ = F±L ◦ F ◦ (F±L )−1. The
Ld Ld d Ld d
fact that the definitions of F˜ are equivalent for the + and − case is implied by the following
Ld
commutative diagram (Theorem 1.5.2 in [11])
(q ,q ) (cid:31) FLd // (q ,q )
0 1 1 2
(cid:127)? (cid:127)?? (cid:127)? (cid:127)??
(6) F−Ld(cid:127)(cid:127)(cid:127)(cid:127)(cid:127) ?????F+Ld (cid:127)(cid:127)(cid:127)(cid:127)(cid:127) ?????F+Ld
(cid:127)(cid:127) (cid:127)(cid:127)(cid:127)(cid:127) ????(cid:31)(cid:31) (cid:127)(cid:127) (cid:127)(cid:127)(cid:127)(cid:127) F−Ld ????(cid:31)(cid:31)
(cid:127) (cid:127)
(q ,p ) (cid:31) // (q ,p ) (cid:31) // (q ,p )
0 0 1 1 2 2
F˜Ld F˜Ld
In coordinates, F˜ :(q ,p ) 7→ (q ,p ), where
Ld 0 0 1 1
(7) p = −D L (q ,q ), p = D L (q ,q ).
0 1 d 0 1 1 2 d 0 1
A numerical quadrature can be used to approximate the integral in (4). However, the functional
form of the solution curve q(t) is required when applying a quadrature rule and it is, in general,
unknown. Inpractice,onecanchooseaninterpolatingfunctionontheinterval[0,h]passingthrough
q , q . Then,aquadraturerulecan beapplied to theintegral of theLagrangian evaluated along the
0 1
interpolating function. This approach fits the general framework of Galerkin integration methods.
In more detail, for the construction of Galerkin Lagrangian variational integrators, one replaces the
path space C([0,T],Q), which is an infinite-dimensional function space, with a finite-dimensional
function space, Cs([0,T],Q). Commonly, one uses polynomial approximations to the trajectories,
letting
Cs([0,h],Q) = {q ∈ C([0,h],Q)|q is a polynomial of degree ≤ s}.
4 MELVINLEOKANDTATIANASHINGEL
An approximate action S(q) :Cs([0,h],Q) → R is
s
S(q) = h b L(q(c h),q˙(c h)),
i i i
i=1
X
where c ∈ [0,1] are quadrature points, b are quadrature weights, i = 1,...,s. The Galerkin
i i
discrete Lagrangian is
(8) L (q ,q ) = ext S(q).
d 0 1
q∈Cs([0,h],Q)
In particular, for higher-order methods one takes q ∈ Cs([0,h],Q) in the form
s
q(τh;qν,h) = qν˜l (τ),
0 0 ν,s
ν=0
X
where qν, ν = 1,...,s− 1, are the internal stages, ˜l (τ) are the Lagrange basis polynomials of
0 ν,s
degree s defined on the interval [0,1]. Then, the integration scheme (q ,p )7→ (q ,p ) is given by
0 0 1 1
s
−p = h b ∂L(c h)˜l (c )+ 1∂L(c h)˜l˙ (c ) ,
0 i i 0,s i i 0,s i
∂q h ∂q˙
i=1 (cid:20) (cid:21)
X
s
0= h b ∂L(c h)˜l (c )+ 1∂L(c h)˜l˙ (c ) ,ν = 1,s−1
i i ν,s i i ν,s i
∂q h ∂q˙
i=1 (cid:20) (cid:21)
X
s
p = h b ∂L(c h)˜l (c )+ 1∂L(c h)˜l˙ (c )
1 i i s,s i i s,s i
∂q h ∂q˙
i=1 (cid:20) (cid:21)
X
It has been established (see [6; 11]) that the above scheme is equivalent to a symplectic partitioned
Runge-Kutta method. Note that it yields a discrete solution that is, in general, only piecewise
regular. In what follows, we introduce a construction of higher-order variational integrators with
improved regularity across nodal times, which has as its natural variables the position, and its
derivatives at only the nodal times, without the use of internal stages.
3. Quadrature
Here,wepresentaquadratureformulathatwewilluseinourconstructionofahigh-orderdiscrete
Lagrangian. The advantage of this particular rule is that it only involves function evaluations at
the endpointsof the interval. When fast adaptive treecodes are used in conjunction with automatic
differentiation techniques [15], it is more efficient to obtain higher-order approximations using
higher-derivative information at the endpoints, rather than evaluating the integrand at a number
of internal stages.
Theorem 1. (Euler–Maclaurin quadrature formula)[1] If f is sufficiently differentiable on
(a,b), then for any m > 0
b θ N−1
f(x)dx = f(a)+2 f(a+kθ)+f(b)
2
Za " k=1 #
X
m
B B
− 2l θ2l f(2l−1)(b)−f(2l−1)(a) − 2m+2 Nθ2m+3f(2m+2)(ξ)
(2l)! (2m+2)!
Xl=1 (cid:16) (cid:17)
where B are the Bernoulli numbers, θ = (b−a)/N and ξ ∈ (a,b).
k
PROLONGATION-COLLOCATION VARIATIONAL INTEGRATORS 5
h
Let us apply Theorem 1 to approximate an integral f(x)dx in the simplest case when N = 1.
0
It is easy to see that we obtain the following quadrature rule
R
m
h B
(9) K(f)= [f(0)+f(h)]− 2l h2l f(2l−1)(h)−f(2l−1)(0) ,
2 (2l)!
Xl=1 (cid:16) (cid:17)
and the error of approximation is O(h2m+3).
4. The Prolongation-Collocation Method
In this section, we explain the construction of the discrete Lagrangian based on Hermite inter-
polation and the Euler–Maclaurin quadrature formula.
Motivation for Prolongation-Collocation Approach. The variational characterization of the
exact discrete Lagrangian (1) naturally leads to the variational Galerkin discrete Lagrangian (8),
where the infinite-dimensional function space C([0,h],Q) is replaced by a finite-dimensional sub-
space, and the integral is approximated by a quadrature formula. While this leads to a com-
putable discrete Lagrangian, one does not necessarily obtain an optimally accurate discrete La-
grangian whose variational order is related to the best approximation properties of the chosen
finite-dimensional function space. In particular, one finds that the variational Galerkin extremal
curves do not necessarily approximate the higher-derivatives of the Euler–Lagrange solution curves
with adequate accuracy.
Inretrospect,thefactthatthevariationalGalerkinapproachdoesnotreadilyleadtocomputable
discreteLagrangianswithprovableapproximationpropertiesisnottoosurprising. Byconstruction,
variational Galerkin discrete Lagrangians associated with a sequence of finite-dimensional function
spaces involve extremizers of a sequence of functionals. Since the sequence of finite-dimensional
function spaces converges to C([0,h],Q), the sequence of functionals converges to the functional
that appears in the variational characterization of the exact discrete Lagrangian. However, it is
unclearthatthesequenceofextremizersconverges totheextremizerofthelimitingfunctional, since
that corresponds to Γ-convergence [3] of the sequence of functionals. The issue of optimal rates of
convergence of the computable discrete Lagrangians involves establishing rates of convergence of
extremizers in terms of approximation rates of the finite-dimensional function spaces, which is an
even more complicated process.
As an alternative, we adopt the characterization of the exact discrete Lagrangian in terms of
the Euler–Lagrange solution curve, and construct a discrete curve which approximates higher-
derivatives of the Euler–Lagrange solution curve to an adequate level of accuracy. The latter is
explored in detail in Section 5.
Hermite Interpolation and Prolongation-Collocation. Wecommencebyreplacingq(t)in(4)
by its Hermite interpolant which is obtained by constructing a polynomial q (t) such that values of
d
q(t)andanynumberofitsderivativesatgivenpointsarefittedbythecorrespondingfunctionvalues
and derivatives of q (t). In this paper we are concerned with fitting function values of q(t) and
d
its derivatives at the end-points of the interval [0,h]. Consequently, a so-called two-point Hermite
interpolant q (t) of degree d= 2n−1 can be used, which has the form
d
n−1
(10) q (t) = q(j)(0)H (t)+(−1)jq(j)(h)H (h−t) ,
d n,j n,j
Xj=0(cid:16) (cid:17)
where
tj n−j−1 n+s−1
(11) H (t) = (1−t/h)n (t/h)s
n,j j! s
s=0 (cid:18) (cid:19)
X
6 MELVINLEOKANDTATIANASHINGEL
are the Hermite basis functions. Note that for n= 1, the interpolant is a straight line joining q(0)
and q(h). By choosing one of the simple quadrature rules to discretize the integral in (4) (e.g., the
midpoint rule or trapezoidal rule), one obtains a class of well-known integrators which are at most
second-order (see [11]). Therefore, the first nontrivial case of interest is n = 2, where we assume
that the position and velocity data at the end points are available. From now on, we only consider
n ≥ 2 when applying the Hermite interpolation formula. The detailed derivation of (10) can be
found, for example, in [4]. By construction,
q(r)(0) = q(r)(0), q(r)(h) = q(r)(h), r = 0,1,...n−1.
d d
Except for the step-size h, the discrete Lagrangian L (q ,q ,h) should only depend on q ≈ q(0),
d 0 1 0
q ≈ q(h). Therefore, letting q (0) = q and q (h) = q , we need to approximate the higher-order
1 d 0 d 1
derivatives of q(t) by expressions that only depend on q ,q . One natural approach, which is often
0 1
found in the literature, is to use finite differences. In this work, we propose to apply the idea of
collocation in conjunction with the Euler–Lagrange equations (2). The benefits of this approach
will be exemplified later when discussing the variational error analysis of the proposed class of
numerical integrators (see Section 5).
The collocation approach [5] is well-known in the theory of initial and boundary value problems
for ODEs [2; 7]. Roughly speaking, the technique consists of determining the unknown parameters
of a parameterized curve by requiring q (t) to satisfy the ODE at a given set of points (collocation
d
points). To define q (t) uniquely, one sets the number of collocation points to be equal to the
d
number of the available degrees of freedom. In our approach we use the method of collocation in a
slightly unusual manner. In particular, since the parameters in (10) correspond to the derivatives
of the solution curve q(t) at the end points of the interval [0,h], we are going to use t = 0 and t = h
as collocation points for the Euler–Lagrange equations (2) and consider the prolongation [12] of the
Euler–Lagrange equations in order to generate a sufficient number of conditions. In other words,
we increase the number of equations under consideration (not the number of collocation points) to
match the number of degrees of freedom.
For example, consider the case of the quintic Hermite interpolation, i.e., set n = 3 in (10). For
separable Lagrangians of the form L(q,q˙) = 1mq˙2 −V(q), where m is the mass and V(q) is the
2
potential energy term, the Euler–Lagrange equations (2) become a second-order ODE of the form
q¨(t)= f(q(t)),
and its first-order prolongation can be expressed as
q(3)(t)= f′(q(t))q˙(t).
We set the boundary conditions q (0) = q and q (h) = q and the collocation conditions
d 0 d 1
q¨ (0) = f(q (0)), q(3)(0) = f′(q (0))q˙ (0),
d d d d d
q¨ (h) = f(q (h)), q(3)(h) = f′(q (h))q˙ (h).
d d d d d
The above conditions constitute the system of six equations, which uniquely determines the fifth-
degree polynomial q (t) in the form of (10).
d
Ingeneral,fortheHermitepolynomialofdegree2n−1,onewouldneedtodifferentiatetheEuler–
(3) (n)
Lagrangeequation n−2times,thusderivingasystemofn−1equationsforq¨ (t),q (t),...,q (t).
d d d
Evaluated at 0 and h, these (together with the Euler–Lagrange equations) give 2n−2 collocation
equations, which together with boundary conditions (q (0) = q ,q (h) = q ) constitute a sufficient
d 0 d 1
number of conditions to determine the interpolant q (t) uniquely. Note that for large n, the
d
system of collocation conditions becomes nonlinear. Since the second and higher-order derivatives
(j) (j)
q (0),q (h) are given explicitly, the system can be recursively reduced to two implicit equations
d d
PROLONGATION-COLLOCATION VARIATIONAL INTEGRATORS 7
involving q˙ (0),q˙ (h). In this case, one would need to make use of a nonlinear root solver, such as
d d
the Newton–Raphson method, to determine q˙ (0),q˙ (h).
d d
Further, in order to discretize the integral in (4), we apply the Euler–Maclaurin quadrature
formula(9). Recallthattheformulainvolvesderivativesoftheintegrand,inourcasetheLagrangian
L, with respect to the independent variable evaluated at the end-points of the integration interval.
The latter, however, does not require the extensive computations typically associated with the
interpolating polynomial due to the use of the Hermite interpolation formula and the collocation
idea explained above. In more detail, we write
h h
L(q (t),q˙ (t))dt ≈ (L(q (0),q˙ (0))+L(q (h),q˙ (h)))
d d d d d d
2
Z0
m B d2l−1 d2l−1
− 2l h2l L(q (t),q˙ (t)) − L(q (t),q˙ (t)) .
(2l)! dt2l−1 d d dt2l−1 d d
l=1 (cid:18) (cid:12)t=h (cid:12)t=0(cid:19)
X (cid:12) (cid:12)
Provided that the degree of the interpolating polynomial is 2n(cid:12)−1, we choose m ≤ ⌊n/2(cid:12)⌋, where
(cid:12) (cid:12)
the brackets denote the greatest integer lower bound for n/2. So for even n, ⌊n/2⌋ = n/2 and for
odd n, ⌊n/2⌋ = (n−1)/2. This choice of m is justified by observing that the expressions for
d2l−1
L(q (t),q˙ (t)) , l = 1,2,...,m, τ = 0,h,
dt2l−1 d d
(cid:12)t=τ
(cid:12)
include the derivatives of qd(t) up to orde(cid:12)r 2⌊n/2⌋−1+1 ≤ n, which satisfy the corresponding
(cid:12)
collocation conditions.
Prolongation-Collocation Discrete Lagrangian. The Prolongation-Collocation discrete La-
grangian is defined as follows,
h
(12) L (q ,q ,h) = (L(q (0),q˙ (0))+L(q (h),q˙ (h)))
d 0 1 d d d d
2
⌊n/2⌋
B d2l−1 d2l−1
− 2l h2l L(q (t),q˙ (t)) − L(q (t),q˙ (t)) ,
(2l)! dt2l−1 d d dt2l−1 d d
l=1 (cid:18) (cid:12)t=h (cid:12)t=0(cid:19)
X (cid:12) (cid:12)
where q (t) ∈ Cs(Q) is determined by the boundary and prolon(cid:12)gation-collocation conditio(cid:12)ns,
d (cid:12) (cid:12)
q (0) = q q (h) = q ,
d 0 d 1
q¨ (0) = f(q ) q¨ (h) = f(q ),
d 0 d 1
(13) q(3)(0) = f′(q )q˙ (0) q(3)(h) = f′(q )q˙ (h),
d 0 d d 1 d
. .
. .
. .
dn dn
(n) (n)
q (0) = f(q (t)) q (h) = f(q (t))
d dtn d d dtn d
(cid:12)t=0 (cid:12)t=h
(cid:12) (cid:12)
One can include fewer than ⌊n/2⌋ terms i(cid:12)n the summation in (12). H(cid:12)owever, this will have an
(cid:12) (cid:12)
impact on the variational order of the corresponding integrator as will be further discussed in
Section 5. The system of equations (13) completely defines the discrete Lagrangian (12). Note
that we used the second-order ODE q¨(t)= f(q(t)) as a prototype of the Euler–Lagrange equations
for simplicity of notation only. The same idea applies for any smooth, not necessarily separable,
Lagrangian function and the corresponding Euler–Lagrange equations.
Given the initial conditions (q ,p ), the variational integrator has the form
0 0
p = −D L (q ,q ),
k 1 d k k+1
(14)
p = D L (q ,q ), k = 0,1,...,
k+1 2 d k k+1
8 MELVINLEOKANDTATIANASHINGEL
and defines a one-step map (q ,p ) 7→ (q ,p ). Generally, the equation p = −D L (q ,q )
k k k+1 k+1 k 1 d k k+1
together with thesystem of collocation conditions (13)can bereducedto a system of implicit equa-
tions withrespectto q ,q˙ (t ),q˙ (t ). As soonas thesolution isobtained via someappropriate
k+1 d k d k+1
nonlinear root-finding method, it is inserted into the equation p = D L (q ,q ).
k+1 2 d k k+1
When the Lagrangian has a relatively simple form, it makes sense to compute the expression for
thediscreteLagrangian(12)symbolically,whichcanbedoneusingthesymbolicmoduleinMatlab
or symbolic software such as Mathematica or Maple. Having computed a priori closed-form
expressions for the right-hand side in (14), makes the implementation of the integrator particularly
simple and fast. The reader is referred to Section 6 for some numerical examples.
5. Variational Order Calculation
The construction of variational integrators in the Galerkin framework naturally leads to the
question of how it can be reconciled with the results from approximation theory of function spaces
and numerical analysis of quadrature schemes. In particular, our goal is to explore the way the
quantitative characteristics of the approximation errors enter the calculation of the convergence
order of the respective integrators. Variational error analysis provides the right framework to
pursue this goal.
The variational error analysis introduced in [11], and refined in [14], is based on the idea that
rather than considering how closely the numerical trajectory matches the exact flow, one can con-
sider how the discrete Lagrangian approximates the exact discrete Lagrangian (1) which generates
the exact flow map of the Euler–Lagrange equations. In other words, we are looking at the approx-
imation error in
h
L (q(0),q(h)) ≈ ext L(q(t),q˙(t))dt = LE(q(0),q(h),h).
d d
q∈C([0,h],Q) Z0
q(0)=q0,q(h)=q1
We say that a given discrete Lagrangian is of order r if there exist an open subset U ⊂ TQ with
v
compact closure and constants C and h > 0 so that
v v
(15) kL (q(0),q(h),h)−LE(q(0),q(h),h)k ≤ C hr+1
d d v
for all solutions q(t) of the Euler–Lagrange equations with initial condition (q(0),q˙(0)) ∈ U and
v
for all h ≤ h . In [11], the authors prove the equivalence of (15) (cf. Theorem 2.3.1 in [11]) to:
v
(i) the discrete Hamiltonian map F˜ being of order r;
Ld
(ii) the discrete Legendre transforms F±L being of order r.
d
In particular, the discrete Hamiltonian map is of order r if
(16) kF˜ (q(0),p(0),h)−F˜ (q(0),p(0),h)k ≤ C˜ hr+1,
Ld LEd v
for all solutions (q(t),p(t)) of the Hamilton’s equations with initial condition (q(0),p(0)) ∈ U ⊂
w
T∗Q and for all h ≤ h . The order of the discrete Legendre transforms is defined analogously.
w
Recall from the diagram (6) that
F˜ : (q ,p )7→ (q ,p ).
Ld 0 0 1 1
By construction, F˜ (q(0),q(h),h) produces the values (q(h),p(h)) corresponding to the exact so-
LE
d
lution of the Hamiltonian system (3), whereas F˜ (q(0),q(h),h) produces the approximate values
Ld
(q ,p ), q ≈ q(h), p ≈ p(h). To summarize, the estimate (16) provides the local order of con-
1 1 1 1
vergence of the discrete trajectory (q ,p ) to the exact flow (q(t),p(t)) of the Hamiltonian vector
k k
field. This order is the same as the order to which the discrete Lagrangian approximates the exact
discrete Lagrangian, which we focus on.
Before we explore the inequality (15) for the Prolongation-Collocation discrete Lagrangian dis-
(j)
cussed in Section 4, we would like to establish the approximation error of q (t) in comparison to
d
PROLONGATION-COLLOCATION VARIATIONAL INTEGRATORS 9
q(j)(t) for j = 1,2...,n, where q(t) is the exact solution of the Euler–Lagrange equation, and q (t)
d
is the Hermite interpolating polynomial (10) of degree 2n−1, constructed by letting q (0) = q(0)
d
and q (h) = q(h), and imposing the prolongation-collocation conditions discussed in Section 4 at
d
the endpoints. Note that this can be a difficult task in general, since the complexity of the col-
location procedure escalates with the degree of the Hermite polynomial. However, we are only
interested in the approximation order at the end-points of the interval [0,h], in which case the
analysis is straightforward.
Lemma 1. For q(t), q (t) as above, if q˙ (τ) = q˙(τ)+O(hp) for some p > 0 and τ = 0,h, then
d d
q(j)(τ) = q(j)(τ)+O(hp), j = 3,...,n.
d
Proof. Note that q¨(τ) coincides with q¨ (τ) by construction. Therefore, we consider j starting
d
from 3. As before, we restrict the proof to the case of a separable Lagrangian. Indeed, since the
Euler–Lagrange equation is equivalent to the second-order ODE
(17) q¨(t) =f(q(t)),
it follows that for τ = 0,h,
q(3)(τ)−q(3)(τ) = f′(q (τ))q˙ (τ)−f′(q(τ))q˙(τ)
d d d
= f′(q(τ))q˙ (τ)−f′(q(τ))q˙(τ) = f′(q(τ))(q˙ (τ)−q˙(τ)) = O(hp),
d d
provided there exists a uniform bound for f′. Consecutively differentiating (17) and substituting
the corresponding expressions for lower order derivatives, one can see that q(j)(τ) (and q(j)(τ)) can
d
be represented as a polynomial in powers of q˙(τ) (resp. q˙ (τ)) with coefficients which only depend
d
on q(τ) = q (τ),
d
q(4)(τ) = f′′(q(τ))q˙(τ)2+f′(q(τ))f(q(τ))
q(5)(τ) = f(3)(q(τ))q˙(τ)3 +q˙(τ)(3f′′(q(τ))f(q(τ))+f′(q(τ))2)
.
.
.
As soon as f has bounded higher-order derivatives, the above formulas imply that the order of the
approximation of q(j)(τ) by q(j)(τ) as a function of h is equal to the order of the approximation of
d
q˙(τ) by q˙ (τ). We consider these expressions up to the j = n case, where n is determined by the
d
number of collocation equations that were used to compute q . (cid:3)
d
In the next lemma, we determine the value of p in the relation q˙ (τ) = q˙(τ)+O(hp) in terms of
d
the degree of the polynomial q (t).
d
Lemma 2. Consider a polynomial q (t) of degree d= 2n−1 given by the formula
d
n−1
q (t) = a H (t)+(−1)ja H (h−t) ,
d 0j n,j 1j n,j
j=0
X(cid:0) (cid:1)
where H (t) are the basis polynomial functions (11). By construction,
n,j
(j) (j)
q (0) = a , q (h) = a , j =0,...n−1.
d 0j d 1j
We let a = q(0), a = q(h). The coefficients a , a , j = 1,2,...,n − 1 are obtained from
00 10 0j 1j
the system of equations consisting of the Euler–Lagrange equation (2) and its prolongations. In
particular, these are n−1 differential equations evaluated on q (t) at t = 0 and t = h. Then for
d
τ = 0,h,
q˙ (τ) = q˙(τ)+O(h2n−1).
d
10 MELVINLEOKANDTATIANASHINGEL
Proof. Let q(t) ∈ C2n([0,h],Q) be the solution of the Euler–Lagrange equation (2) with boundary
conditions q(0) and q(h). Then, q¨(t) can be written in the form,
q¨(t) = P [q¨](t)+R [q¨](t),
2n−3 n−1
where
n−2
P [q¨](t) = q(j+2)(0)H (t)+(−1)jq(j+2)(h)H (h−t) ,
2n−3 n−1,j n−1,j
Xj=0(cid:16) (cid:17)
q2n+2(ξ)
R [q¨](t) = tn−1(h−t)n−1, ξ ∈ (0,h),
n−1
(2n−2)!
and H (t) are the basis polynomials defined by (11). See [4] for the proof of the above formula
n−1
for sufficiently smooth functions. Note that q¨ (t) is a polynomial of degree 2n−3 to which the
d
same formula can be applied. In the latter case, the remainder term is identically zero. Hence,
n−2
q¨ (t) = q(j+2)(0)H (t)+(−1)jq(j+2)(h)H (h−t) .
d d n−1,j d n−1,j
Xj=0(cid:16) (cid:17)
Subtracting q¨ (t) from q¨(t) gives
d
(18) q¨(t)−q¨ (t)
d
n−2
= (q(j+2)(0)−q(j+2)(0))H (t)+(−1)j(q(j+2)(h)−q(j+2)(h))H (h−t)
d n−1,j d n−1,j
Xj=1(cid:16) (cid:17)
+R [q¨](t).
n−1
Next, we integrate the above expression from 0 to h. The left-hand side of the equation becomes
h
[q¨(t)−q¨ (t)]dt = [q˙(h)−q˙ (h)]−[q˙(0)−q˙ (0)].
d d d
Z0
Further, observe that
h h
H (t)dt = H (h−t)dt = C h2, j = 1,...n−2
n−1,j n−1,j j
Z0 Z0
h
tn−1(h−t)n−1 = Ch2n−1,
Z0
where C ,C > 0 are constants which do not depend on h. Now, let q˙ (0) = q˙(0) + O(hp) and
j d
q˙ (h) = q˙(h) + O(hp), where we wish to determine the value of p. It is easy to see that after
d
integrating both sides of (18) and rewriting it in terms of the order conditions we arrive at
O(hp)= O(hp+2)+O(h2n−1).
Note that weused Lemma1to estimate thehigher-order derivatives in (18). Itfollows immediately
that p = 2n−1 and the proof is finished. (cid:3)
We are now ready to prove the main result of this section.
Theorem 2. Assume that a Lagrangian function L : TQ → R is sufficiently smooth and its
partial derivatives are uniformly bounded. Then, the discrete Lagrangian L (q ,q ,h) constructed
d 0 1
according to (12), (13) with q = q(0) and q = q(h), approximates the exact discrete Lagrangian
0 1
LE(q(0),q(h),h) with order 2⌊n/2⌋ +3, for n ≥ 3. In particular, when n is even, ⌊n/2⌋ = n/2,
d
and the order is n+3. For odd n, ⌊n/2⌋ = (n−1)/2, and the order is n+2. For n= 2, the order
is equal to 3.