Quadratic choreographies

Applied Numerical Mathematics 75 (2014) 108–122

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

Quadratic choreographies P. Ryckelynck a,b,c , L. Smoch a,b,c,∗ a b c

ULCO, LMPA, F-62100 Calais, France Université Lille Nord de France, F-59000 Lille, France CNRS, FR 2956, France

a r t i c l e

i n f o

Article history: Available online 10 July 2013 Keywords: Calculus of variations Functional equations Discretization Quadratic eigenvalue problems Periodic and almost-periodic solutions

a b s t r a c t This paper addresses the classical and discrete Euler–Lagrange equations for systems of n particles interacting quadratically in Rd . By highlighting the role played by the center of mass of the particles, we solve the previous systems via the classical quadratic eigenvalue problem (QEP) and its discrete transcendental generalization. Next, we state a conditional convergence result, in the Hausdorff sense, for the roots of the discrete QEP to the roots of the classical one. At last, we focus especially on periodic and choreographic solutions and we provide some numerical experiments which confirm the convergence. © 2013 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction This paper seeks to continue the development of the theory for the discrete calculus of variations which was initiated by Cresson et al., see [4,5]. It consists originally in replacing the derivative x˙ (t ) of the dynamic variable x(t ) defined on [t 0 , t f ] with a 2N + 1 terms scale derivative

2x(t ) =

N  γj j =− N

ε

x(t + j ε )χ− j (t ),

t ∈ [t 0 , t f ],

γj ∈ C

(1)

where χ j (t ) denotes the characteristic function of [t 0 , t f ] ∩ [t 0 + j ε , t f + j ε ], for some time delay ε . We consider a Lagrangian L of n particles in Rd , where d denotes the “physical” dimension. The principle of least action may be extended to the case of non-differentiable dynamic variables, see [4,5] and [12,13] for various statements of this principle. The equations of motion may be returned as the following two dynamic sets of equations

x¨ j (t ) = F j (x1 , . . . , xn , x˙ 1 , . . . , x˙ n ),

(2)

−2 2x j (t ) = F˜ j (x1 , . . . , xn , 2x1 , . . . , 2xn )

(3)



where 2 denotes the analogous operator to 2 where the coefficients γ are replaced with γ− , and the functions F j , F˜ j are built from the specific interaction between the particles. While (2) is a set of ODE, (3) consists in functional difference equations. Note by the way that for conservative systems, the forces in (2) and (3) do not depend explicitly on the variables x˙ k and 2xk . In contrast, the use of variational integrators (see [8,11] for instance) approximates the action integral by a finite sum defined in some finite-dimensional vector space of which the elements are the trajectories samples. The minimization of the

*

Corresponding author at: ULCO, LMPA, F-62100 Calais, France. E-mail addresses: [email protected] (P. Ryckelynck), [email protected] (L. Smoch).

0168-9274/$36.00 © 2013 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2013.03.010

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

109

action leads to recurrence equations and not to ODE nor to functional equations. Nevertheless, these recurrence equations of motion and the numerical scheme that we devise for solving (3) are quite similar, provided t ∈ [t 0 + 2N ε , t f − 2N ε ] and L is autonomous. We investigate the existence of pseudo-periodic solutions of (2) and (3) of the shape

u(t ) = u0 +

K 

e λ t u ,

(4)

=1

where u0 and u constitute a family of K + 1 vectors of Cd and (λ ) ∈ (C ) K is a sequence of K distinct complex numbers. When the Lagrangian is quadratic and autonomous so that the forces F j and F˜ j are linear, it is well known that quadratic eigenvalue problem (QEP) may be applied to solve (2) (see [14]). In order to solve (3), we introduce a transcendental eigenvalue problem which seems to be new and allows to investigate qualitative features of the solution. The rest of this paper is organized as follows. Section 2 is devoted to the derivation of the classical and discrete  Euler– Lagrange equations (respectively abbreviated as C.E.L. and D.E.L.) and highlights the role played by xs (t ) = j x j (t ), or equivalently, by the center of mass n1 xs (t ). In Section 3, we present a method for solving the equations of motion for generic Lagrangians for C.E.L. as well as D.E.L. The first step of this method determines xs (t ) from some generalized (quadratic or transcendental) eigenvalue problem. The second step seeks x j (t ) from xs (t ) by solving another eigenvalue problem. Section 4 is devoted to the convergence of the generalized eigenvalue problem linked to the D.E.L. as ε tends to 0. Section 5 deals with the existence and the features of periodic and choreographic solutions. Section 6 presents some numerical experiments illustrating the phenomenon of convergence as ε tends to 0, for some various operators 2. Finally, Section 7 gives an overall view and conclusions. 2. Equations of motion for symmetric quadratic Lagrangians of n particles systems in R d 2.1. Actions in the continuous and discrete settings We denote by C p w the space of the functions x : [t 0 , t f ] → Rd continuous on each interval ]t 0 + j ε , t 0 + ( j + 1)ε [∩[t 0 , t f ] for all j  0 and small enough, i.e. j 

t f −t 0

ε

, and continuous at t 0 and t f . The subspace of C p w the elements of which

are the functions x : [t 0 , t f ] → R twice-continuously is denoted by C 2 . The space C p w is endowed with the ordinary scalar  tf product f, g = t0 t f(t ).g(t ) dt. If X = (x1 , . . . , xn ) denotes a system of n functions in C p w , we may think of X as the set of d

dynamic variables describing the state of a system of n interacting particles in Rd . We consider the actions Acont and Adisc of the shape

t f Acont (X) =



 ˙ (t ) dt , L X(t ), X

t f Adisc (X) =

t0

  L X(t ), 2X(t ) dt .

(5)

t0

These actions are defined respectively on the subspaces of (C 2 )n and C np w the elements X of which satisfy Dirichlet conditions at t 0 and t f . In contrast, variational integrators [8,11] generate discrete action sums of the following shape. Let Y j ∈ (Rd )n be the sample Y j = X( j + t ε ) and M be the number of time samples. If Y = (Y0 , . . . , Y M ) is the approximation of the whole trajectory over the interval [t 0 , t f ], then the discrete action is AVI (Y) =

 M −1 j =0

Ld (Y j , Y j+1 , ε ).

The study of convergence of the discrete actions is two-fold oriented: first, to prove that Adisc (X) and AVI (Y) converge to Acont (X) as the time step tends to 0, and second, to prove that accordingly, the solutions of the equations of motion for AVI converge to those for Acont . The case of variational integrators is discussed in details in [8, Chap. IX]. 2.2. n-Body problem associated to quadratic Lagrangians We introduce a general quadratic Lagrangian of n particles in Rd , compatible with discrete symmetries of the system. Let J 1 , . . . , J 5 ∈ C 1 ([t 0 , t f ], Rd×d ), J 1 , . . . , J 4 be symmetric matrices and J 6 , J 7 ∈ C 1 ([t 0 , t f ], Rd ). From now on, we drop t from the formulas when it is clear enough. For an isolated particle with position x and velocity y we may set

L1 (x, y) =

1t 2

y J 1y +

1t 2

x J 2 x + t x J 5 y + t J 6 y + t J 7 x.

Next, two particles with positions x j , xk , and velocities y j , yk are interacting for pairs in conformity with the following Lagrangian

L2 (x j , y j , xk , yk ) = t y j J 3 yk + t x j J 4 xk . Therefore, the Lagrangian of the whole system is

110

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

L=

n 

L1 (x j , x˙ j ) +



j =1

L2 (x j , x˙ j , xk , x˙ k ).

(6)

j =k

Eq. (6) is a paradigm for many models arising in physics such as electromagnetism, mechanics of continua, lattices but also in quantum mechanics and relativity. When all the matrices are constant, J 6 and the symmetric part of J 5 contribute to null-Lagrangians and do not receive any interpretation. The diagonal part of J 1 contains masses while the nondiagonal part is linked to anisotropy of the frame. For each particle, J 2 and J 7 may be thought respectively as a stiffness matrix and a constant force exerted on it. In classical mechanics (see e.g. [1]), when one considers centrifugal forces, the opposite of the skew-symmetric part of J 5 is equal to the matrix of the cross-product with the angular momentum in the rotating frame. Considering at last the interaction of two particles, the matrix J 4 occurs in harmonic oscillator lattices or Hamiltonian networks of weakly coupled oscillators. In contrast, J 3 is rather rare and appears for example in the no-core shell model. The previous model of Lagrangian (6) takes into account the fact that the system is symmetric w.r.t. the permutations of the particles, i.e. under the symmetric group Sn . 2.3. Hamilton’s principle Lemma 2.1. The adjoint operator of 2 w.r.t. the product ·,· is the operator 2 obtained from (1) by substituting γ− for γ . Proof. Without loss of generality, by using bilinearity, we may suppose d = 1. In this case, for all f , g ∈ C p w , we have

 f , 2g =

t f N  γk k=− N

=

ε

f (s − kε ) g (s)χ0 (s)χk (s) ds

t0

t f N  γk k=− N

f (t ) g (t + kε )χ−k (t ) dt

t0

t f N  γk k=− N

=

ε

ε





g (t ) f (t − kε )χk (t ) dt = 2 f , g ,

t0

by using some simple changes of variables and indices. Theorem 2.1. Let X = (x1 , . . . , xn ) in C np w and xs =

2

n

j =1 x j .

A necessary and sufficient condition for X to be a critical point of Acont in (C 2 )n is that X satisfies the dynamic system

  ( J 1 − 2 J 3 )¨x j + (−2 ˙J 3 + t J 5 + ˙J 1 − J 5 )˙x j + t ˙J 5 − J 2 + 2 J 4 x j = −2 J 3 x¨ s − 2 ˙J 3 x˙ s + 2 J 4 xs + ( J 7 − ˙J 6 ), ∀ j .

A necessary and sufficient condition for X to be a critical point of Adisc in of equations

C np w

(7)

is that X satisfies the linear functional recurrence system

    2 ( J 1 − 2 J 3 )2x j + 2 t J 5 x j + J 5 2x j + ( J 2 − 2 J 4 )x j   = −22 ( J 3 2xs ) − 2 J 4 xs − 2 J 6 + J 7 , ∀ j .

(8)

Proof. The differential equations (7) are straightforward consequences of the classical Euler–Lagrange equations (see for instance [1]). Let X, H ∈ C np w be given, then the mapping η → Adisc (X + ηH) − Adisc (X) is a quadratic polynomial. By using direct expansion, we find that the coefficient of η is equal to

D Adisc (X)(H) =

t f  n t t0



σ j (t )2h j (t ) + t τ j (t )h j (t ) dt

j =1

where

σ j = J 1 2x j + t J 5 x j + 2



J 3 2xk + J 6 = ( J 1 − 2 J 3 )2x j + 2 J 3 2xs + t J 5 x j + J 6 ,

k= j

τ j = J 5 2x j + J 2 x j + 2

 k= j

J 4 xk + J 7 = ( J 2 − 2 J 4 )x j + 2 J 4 xs + J 5 2x j + J 7 ,

(9)

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

since xs =

n

j =1 x j .

111

Using the preceding lemma, we get

D Adisc (X)(H) =

n  

2 σ j + τ j , h j .

(10)

j =1

Since X is a critical point, we have D Adisc (X)(H) = 0, ∀H. Thus, we have for all j the relation 2 σ j + τ j = 0 which amounts to (8) and ends the proof. 2 Remark 2.1. We may easily deduce Eqs. (8) from the general discrete Euler–Lagrange equations stated in [12, Theorem 4.1], that is

2

 ∂L   ∂L  t , X(t ), 2X(t ) + t , X(t ), 2X(t ) = 0, ∂2x j ∂x j

∀ j,

(11)

which hold for arbitrary Lagrangians and arbitrary operators 2. 2.4. Relative motion around the center of mass We notice that Eqs. (7) and (8) are quite uncoupled since the coupling is realized only through the vector xs . We mention two simple consequences of Theorem 2.1. The first one arises from summing all equations in (7) or in (8) over j, and the second one deals with time-independent Lagrangians such that J 5 is skew-symmetric. Corollary 2.1. If X is a solution to (7), then the sum xs satisfies the dynamic system









J 1 + 2(n − 1) J 3 x¨ s + ˙J 1 + 2(n − 1) ˙J 3 + t J 5 − J 5 x˙ s

  − J 2 + 2(n − 1) J 4 − t ˙J 5 xs + n( ˙J 6 − J 7 ) = 0.

(12)

Similarly, if X is a solution to (8), then xs satisfies the functional equation

  2 ( J 1 2xs ) + 2(n − 1)2 ( J 3 2xs ) + 2 t J 5 xs   + J 5 2xs + J 2 + 2(n − 1) J 4 xs + n2 J 6 + n J 7 = 0.

(13)

Corollary 2.2. Suppose that the functions J k (t ), k ∈ {1, . . . , 7}, are constant w.r.t. time and J 5 is skew-symmetric. The systems of equations (12) and (7) simplify respectively into









J 1 + 2(n − 1) J 3 x¨ s − 2 J 5 x˙ s − J 2 + 2(n − 1) J 4 xs = n J 7

(14)

and

( J 1 − 2 J 3 )¨x j − 2 J 5 x˙ j − ( J 2 − 2 J 4 )x j = −2 J 3 x¨ s + 2 J 4 xs + J 7 .

(15)

Similarly, the systems of functional recurrence equations (13) and (8) simplify respectively into













J 1 + 2(n − 1) J 3 2 2xs + J 5 2xs − 2 xs + J 2 + 2(n − 1) J 4 xs

= −n2 (1) J 6 − n J 7

(16)

and

  ( J 1 − 2 J 3 )2 2x j + J 5 2x j − 2 x j + ( J 2 − 2 J 4 )x j = −2 J 3 2 2xs − 2 J 4 xs + 2 J 6 − J 7 .

(17)

˙ ) − U (X) and consequently, it is conservative, i.e. the Remark 2.2. If J 5 = 0, the system has a Lagrangian of the shape T (X energy n

 1 j =1

t x˙ J x˙ − 1 t x J x + t J x˙ − t J x 6 j 7 j j 1 j j 2 j 2 2

is a constant of motion.

 +

 t j =k

x˙ j J 3 x˙ k − t x j J 4 xk



(18)

112

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

3. Solutions of equations of motion in the general case 3.1. Preliminaries on Quadratic Eigenvalue Problems We provide in this section the solutions to problems presented in Corollaries 2.1 and 2.2. From now on we suppose that the vectors and matrices J k , k = 1, . . . , 7, are time-independent and that J 5 is skew-symmetric. By general case, we mean that the set of coefficients ( J k )k=1,...,7 , satisfying the conditions (21) and (25) below is everywhere dense in (Rd×d )5 × (Rd )2 . According to [9,10,14] we define the Quadratic Eigenvalue Problem associated to ( A , B , C ) as the search of the complex roots of the discriminantal equation





det A λ2 + B λ + C = 0

(19)

where the l.h.s. is a polynomial of λ of degree 2d, together with the description of the various kernels ker( A λ2 + B λ + C ). The reader is referred to [14] for a survey of theory applications and algorithms of the QEP. It is a classical old fact [7,9,14] that, if all the roots λ ,  ∈ {1, . . . , 2d}, of A λ2 + B λ + C are distinct and C is invertible, the general solution to the dynamic system A x¨ + B x˙ + C x = k0 has the shape (4), with K = 2d, u0 = C −1 k0 and u ∈ ker(λ2 A + λ B + C ), ∀. When the number of roots is less than 2d, slight more complicated expressions for the solutions may be found in [9,10,14]. Because of Eqs. (12) to (17) and the previous discussion, we may provide, under specific assumptions, the shape of the solutions to C.E.L. and D.E.L. 3.2. The case of C.E.L. We introduce the matrix-valued function

    Pν (λ) := J 1 + 2(ν − 1) J 3 λ2 − 2 J 5 λ − J 2 + 2(ν − 1) J 4

(20)

and the following subsets of C



 Qν := λ ∈ C/ det Pν (λ) = 0 ,

∀ν ∈ R.

Proposition 3.1. We assume that

|Qn | = |Q0 | = 2d, Qn ∩ Q0 = ∅,     det Pn (0) = 0 and det P0 (0) = 0

(21)

where n denotes the number of interacting particles. Then all the solutions xs and x j to (14) and (15) respectively are of the shape (4) with K = 2d. Proof. Since det( J 2 + 2(n − 1) J 4 ) = 0, we may define



xs,0 = −n J 2 + 2(n − 1) J 4

− 1

J7

and we have obviously 0 ∈ / Qn . The first condition |Qn | = 2d guarantees that the solution xs to (14) is of the shape (4) with K = 2d, i.e.

xs (t ) = xs,0 +



e αt xs,α ,

(22)

α ∈Qn

for some convenient vectors xs,α ∈ ker(Pn (α )). Since Qn ∩ Q0 = ∅, the matrix P0 (α ) is invertible for each may define for j ∈ {1, . . . , n}





x j ,α = 2P0 (α )−1 J 4 − α 2 J 3 xs,α ,

α ∈ Qn so we

∀α ∈ Qn .

2 J 4 )−1 (2 J 4 xs,0

Let us set x j ,0 = −( J 2 − + J 7 ). Straightforward computations show that x j ,0 = n1 xs,0 and x j ,α = n1 xs,α ,  ∀α ∈ Qn . Therefore, we have proved that x j ,0 + α ∈Qn e αt x j ,α = n1 xs (t ) is a particular solution to (15). Hence, the general solution to (15) is given by the formula (4) with K = 2d, i.e.

x j (t ) =

1 n

xs (t ) +



e β t xj ,β

(23)

β∈Q0

for some convenient vectors xj ,β ∈ ker(P0 (β)).

2

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

113

Let us consider now the well-posedness of the Dirichlet problem for C.E.L. If all the kernels ker(Pn (α )) for α ∈ Qn and ker(P0 (β)) for β ∈ Q0 occurring in the previous proof are one-dimensional, then the computation of the vectors xs,α and xj ,β amounts to solving linear systems for the abscissas along these kernels. Proposition 3.2. Under the assumptions (21), for all (t 0 , t f ) ∈ R2 except for a meager set, there exists one and only one solution X ∈ (C 2 )n of C.E.L. together with Dirichlet conditions. Moreover, if J 5 = 0 and the following symmetric (nd) × (nd) matrices

⎛ J 2 J3 . . . 2 J3 ⎞ 1 .. ⎟ .. ⎜ . . ⎟ ⎜ 2 J3 J1 J8 = ⎜ ⎟ .. .. ⎝ .. . . 2 J3 ⎠ . 2 J3 . . . 2 J3 J1

and

⎛ J 2 J4 . . . 2 J4 ⎞ 2 .. ⎟ .. ⎜ . . ⎟ ⎜ 2 J4 J2 J9 = ⎜ ⎟ .. .. ⎝ .. . . 2 J4 ⎠ . 2 J4 . . . 2 J4 J2

are definite positive, the solution X is a global strict minimizer of Acont (X). Proof. Let us fix ν = 0 or ν = n. Since det(Pν (λ)) admits exactly 2d distinct roots, then dim ker(Pν (λ)) = 1, ∀λ ∈ Qν and the union of the various ker(Pν (λ)) spans Cd . This follows from the existence of the Smith form for regular QEP, see for example [14, Section 3] or [7,9]. Let us consider first the existence and unicity of xs , when ν = n. Since Proposition 3.1 holds, the formula (22) holds. We choose arbitrarily one eigenvector vα ∈ ker(Pn (α )) for all α ∈ Qn . We are looking for the 2d abscissas ζα ∈ C such that xs,α = ζα vα . The Dirichlet conditions





e αt0 ζα vα = xs (t 0 ) − xs,0 ,

α ∈Qn

e αt f ζα vα = xs (t f ) − xs,0

α ∈Qn

give rise to a linear system of 2d equations w.r.t. 2d unknowns ζα . We proceed in a similar way to compute the 2d abscissas of xj ,β , β ∈ Q0 , occurring in the function x j (t ) − n1 xs (t ), j ∈ {1, . . . , n}, being fixed. Hence, we get n + 1 uncoupled square systems of size 2d (the very last one being useless due to the definition of xs ). The determinant of each system is a polynomial w.r.t. the 4d quantities e αt0 and e αt f with coefficients in the subfield of R generated over the rationals by Q0 ∪ Qn together with the entries of the matrices J k , k ∈ {1, . . . , 5}. By using classical results of algebraic geometry, the set of couples (t 0 , t f ) such that one of these determinants vanishes is meager in R2 . Therefore, for all (t 0 , t f ) in a dense subset of {(x, y )R2 /x < y }, all these systems are Cramer and there is one and only one solution X(t ) of (14) and (15) obeying the convenient Dirichlet conditions. Let us prove the second part of the proposition. If X ∈ (C 2 )n is a critical point of Acont , i.e. X satisfies C.E.L., and if Y ∈ (C 2 )n vanishes at t 0 and t f , then we have

Acont (X + Y) − Acont (X) =

1

t f

t

2



˙ J 8 Y˙ + t Y J 9 Y dt . Y

t0

˙ ), Eqs. (14) and (15) are necessary As a consequence, because the integrand is a positive definite quadratic form w.r.t. (Y, Y and sufficient conditions for a global strict minimum of the action Acont to occur. 2 3.3. The case of D.E.L. Let us extend Proposition 3.1 to the D.E.L. case. As well as the study of autonomous dynamic differential systems leading to QEP, the study of autonomous difference equations leads to transcendental eigenvalue problem associated to the following complicated matrix



  P˜ ν (ε , λ) := − J 1 + 2(ν − 1) J 3

−2N k2N − N  N |k+| N

− J5

N  1 k=− N

ε

γk+ γ kλε e ε2

  (γk − γ−k )ekλε − J 2 + 2(ν − 1) J 4 .

(24)

This matrix is the coefficient of e λt a in (16) when we substitute e λt a for xs (t ), as shows the proof of the following result. Let us introduce the following subsets



 Q˜ ν := λ ∈ C/ det P˜ ν (ε , λ) = 0 ,

∀ν ∈ R.

114

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

Proposition 3.3. We assume that

|Q˜ n | = |Q˜ 0 | = 4Nd, Q˜ n ∩ Q˜ 0 = ∅,     det P˜ n (ε , 0) = 0 and det P˜ 0 (ε , 0) = 0.

(25)

Then there are solutions x˜ s and x˜ j to (16) and (17) respectively of the shape (4) inside the interval [t 0 + 2N ε , t f − 2N ε ]. Proof. The two last conditions (25) imply that 0 ∈ / Q˜ n . If we set ζ = e λ , we see that the quantity P˜ n (ε , λ)ζ 2N is a poly2Nd det(P˜ n (ε , λ)) = 0 gives rise to a polynomial equation w.r.t. ζ of degree nomial w.r.t. ζ of degree 4N. So the equation ζ 4Nd. Computation of the l.h.s. of (16) is performed by using the definition of 2 and 2 and we find

 −2N k2N − N  N |k+| N

+

N  k=− N

ε





1 2





γk+ γ χ (t )χ−k (t ) J 1 + 2(n − 1) J 3 x˜ s (t + kε) + J 2 + 2(n − 1) J 4 x˜ s (t )



1



χ−k (t ) (γk − γ−k ) J 5 x˜ s (t + kε) + n 2 1 J 6 + n J 7 = 0. ε

When t lies in the interval [t 0 + 2N ε , t f − 2N ε ], the various characteristic functions For some arbitrary but not variable t i ∈ [t 0 , t f ], we define the grid

(26)

χk (t ) occurring in (26) are equal to 1.

Gt i ,ε = {t i + mε , m ∈ Z} ∩ [t 0 + 2N ε , t f − 2N ε ]. Both restrictions of x˜ s (t ) classical theory of those order of the recurrence,

and x˜ j (t ) to Gt i ,ε are vector-valued sequences satisfying linear constant matricial recurrences. The systems shows that, provided the  characteristic equation admits a number of roots equal to the x˜ s (t ) has the shape x˜ s (t ) = x˜ s,0 + λ e λt x˜ s,λ , ∀t ∈ Gt i ,ε , for some vectors x˜ s,λ and x˜ s,0 . Here, the

˜ n |. order of recurrence is equal to 4Nd and it is also equal to the number of roots of the characteristic equation which is |Q So we may plug the previous formula into (26) and we find







e λt P˜ n (ε , λ)˜xs,λ + P˜ n (ε , 0)˜xs,0 = n 2 1 J 6 + n J 7 .

(27)

λ t −t

Because the values of the function e λt , on the grid Gt i ,ε , are the numbers e λt i ζ m with m = ε i ∈ N, all the functions e λt on this grid are linearly independent. Indeed, a linear relationship between these functions would give rise to a Vandermonde determinant w.r.t. to the associated distinct numbers ζ . Therefore, every non-constant function of t must vanish in (27), ˜ n . By assumption, P˜ n (ε , 0) is invertible which means that the “phases” λ occurring in x˜ s (t ) are exactly the roots α of Q and P˜ n (ε , α ) is singular. Thus, we may choose x˜ s,0 = nP˜ n (ε , 0)−1 ( J 7 + (2 1) J 6 ) and x˜ s,λ ∈ ker(P˜ n (ε , λ)). Finally, we have determined the general solution x˜ s to (16) on the grid Gt i ,ε , namely

x˜ s (t ) = x˜ s,0 +



e αt x˜ s,α .

(28)

α ∈ Q˜ n

Now, let us deal with x˜ j (t ). This function satisfies the following functional equation, which is similar to (26)

−



1

−2N k2N − N  N |k+| N

ε2

γk+ γ χ (t )χ−k (t )( J 1 − 2 J 3 )˜x j (t + kε)

− ( J 2 − 2 J 4 )˜x j (t ) − J 5 

N  k=− N

1

χ−k (t ) (γk − γ−k )˜x j (t + kε) ε

   = 2 J 3 2 2˜xs (t ) + J 4 xs (t ) − 2 1 J 6 + J 7 . 

(29)

Let us construct a particular solution to (29) for t ∈ Gt i ,ε . By using the previous expression for x˜ s (t ), the r.h.s. of (29) may be rewritten as









2 2 (21) J 3 + J 4 x˜ s,0 − 2 1 J 6 + J 7 + 2

 α ∈Q˜ n

e αt ( J 3 + θα J 4 )˜xs,α

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

115

where

1 

θα = e −αt 2 2e αt = Now, if we substitute

x˜ j (t ) = x˜ j ,0 +



ε2

γk+ j γ j ekαε .

k, j

e αt x˜ j ,α

α ∈Q˜ n

in (29), we note that the l.h.s. of (29) is equal to

P˜ 0 (ε , 0)˜x j ,0 +



e αt P˜ 0 (ε , α )˜x j ,α .

α ∈Q˜ n

Because det(P˜ 0 (ε , 0)) = 0, we may define

 









x˜ j ,0 = P˜ 0 (ε , 0)−1 2 2 (21) J 3 + J 4 x˜ s,0 − 2 1 J 6 + J 7 .

˜ n ∩ Q˜ 0 = ∅, the matrix P˜ 0 (ε , α ) is invertible for each Since Q

α ∈ Q˜ n and we may set

x˜ j ,α = 2P˜ 0 (ε , α )−1 ( J 3 + θα J 4 )˜xs,α .

˜ n . At last, since |Q˜ 0 | = 4Nd, we Similarly to the case of C.E.L., we readily prove that x˜ j ,0 = n1 x˜ s,0 and x˜ j ,α = n1 x˜ s,α , ∀α ∈ Q conclude that the general solution to (29) on the grid Gt i ,ε is given by x˜ j (t ) =

1 n



x˜ s (t ) +

e β t x˜ j ,β

(30)

˜0 β∈Q

for some convenient vectors x˜ j ,β in ker(P˜ 0 (ε , β)). If we drop the requirement that t lies in Gt i ,ε , the functions t → x˜ s (t ) and t → x˜ j (t ) may be extended by the preceding formulas to pseudo-periodic functions t → x˜ s (t ) and t → x˜ j (t ) respectively. Since the equations of motion are autonomous (independent w.r.t. t), these n + 1 functions are solutions to (16) and (17) respectively. Therefore, these functions are of the shape (4) with K = 4Nd and K = 8Nd respectively and the proof is complete. 2 Remark 3.1. The previous proof shows that for all t i ∈ [t 0 , t f ] and each solution X ∈ C np w of D.E.L., we may associate one and

˜ ∈ C np w which agrees with X on Gt ,ε . In contrast, there are solutions of D.E.L. which are only one pseudo-periodic solution X i not of the shape (4) since the delayed functional equations D.E.L. are not well-posed. Remark 3.2. Solving D.E.L. with Dirichlet conditions leads to n + 1 uncoupled linear systems, one of size   ˜ ˜ ˜ 0 dim ker(P0 (ε , β)). If those systems are Cramer, then the pseudoα ∈Q˜ n dim ker(Pν (ε , α )) and the n others of size β∈Q periodic solution to D.E.L. exists and is unique.

4. Convergence issues Let us present some convergence results for the QEP and its discrete generalization. Let us fix

ν , (γk )−N kN , ( J k )1k5 . ε

It is natural at first sight to ask if the matrix-valued function P˜ ν (ε , λ) tends to Pν (λ) locally uniformly w.r.t. λ ∈ C as tends to 0. Note that these matrices are independent on time. Next, we recall the Hausdorff metric





d H ( F 1 , F 2 ) := max max min |x − y |, max min |x − y | , x∈ F 1 y ∈ F 2

x∈ F 2 y ∈ F 1

defined for all nonempty finite subsets F 1 , F 2 ⊂ C. Thus, we naturally investigate the convergence, as Hausdorff sense, of

ε tends to 0, in the

  −1   −1 Q˜ ν = det P˜ ν (ε , .) {0} to Qν = det Pν (.) {0}. To proceed, we need the following lemma, similar to [13, Theorem 6.1]: Lemma 4.1. We assume that N  k=− N

γk = 0 and

N  k=− N

kγk = 1.

(31)

116

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

Then, we have the following results of convergence, locally uniformly on ]t 0 , t f [:

2e λt → λe λt ,

2 e λt → −λe λt and 2 2e λt → −λ2 e λt .

ε →0

ε →0

ε →0

(32)

Proof. By using (26), we see that N  1

2e λt = e λt

k=− N

2 2e λt = e λt

ε

γk ekλε χ−k (t ),

2 e λt = e λt

N  1 k=− N



1

ε2

−2N k2N − N  N |k+| N

ε

γ−k ekλε χ−k (t ),

γk+ γ ekλε χ (t )χ−k (t ).

Since the three properties (32) may be proved in a similar way, let us briefly give the proof of the third one. Let δ > 0 and ε < δ/2N. When t lies in [t 0 + 2N ε , t f − 2N ε ], all characteristic functions are equal to 1. Using the Taylor formula with Lagrange remainder, we get





e −λt 2 2e λt + λ2 e λt 



1

1

ε

|ε |

|21|2 + 2

    3   1 2 |k| |γk+ γ | k γk+ γ  + |ε |λ3  |21||2t | + λ2 1 + 2

where the symbols of summation have the same meaning as in the formula for 2 2e λt . Now, formula (31) are obviously equivalent to 21 = 0 and 2t = 1 inside [t 0 + 2N ε , t f − 2N ε ]. Because



k

2

γk+ γ =

N  

2

2

j + k − 2kj





N 

γ j γk = 2(21)

j ,k=− N

 k

2

γk − 2(2t )2 = −2,

k=− N

the first three terms in the previous upper bound are equal to 0 and the last one tends to zero uniformly on [t 0 + δ, t f − δ] as ε tends to 0. 2 We also need the following theorem of Cucker and Corbalan [6], which extends older results of Weber and Ostrowski to the case of perturbation of polynomials of distinct degrees. Theorem 4.1. Let P ( X ) = a0 X m + a1 X m−1 + · · · + am ∈ C[ X ]\{0}. Let ξ1 , . . . , ξr be its roots in C, with multiplicities μ1 , . . . , μr respectively, and let B1 , . . . , Br be disjoint disks centered at ξ1 , . . . , ξr with radii ε0 and contained in the open disk centered at 0 with radius 1/ε0 . Then, there is a δ ∈ R+ , such that, if |b j − a j | < δ for every 0  j  m, then the polynomial Q ( X ) = b0 X m + b1 X m−1 + · · · + bm has μ j roots (counted with multiplicity) in each B j and deg( Q ) − deg( P ) roots with absolute value greater than 1/ε0 . Hence, we must exclude the 4Nd − 2d divergent roots, as convergence mentioned above.

ε tends to 0, from the set Q˜ ν to prove the second result of

Theorem 4.2. Let us fix a Lagrangian L and ν ∈ C. We assume that γ− N γ N = 0 and that the properties (21), (25) and (31) hold. ˜ ν ∩ K tends to Qν in the Then, P˜ ν (ε , λ) → Pν (λ) locally uniformly in C. Moreover, for all compact neighborhood K ⊂ C of Qν , Q ε→0

Hausdorff sense as ε tends to 0.

ε = 0 and λ ∈ C,  2  P˜ ν (ε , λ) − Pν (λ) = J 1 + 2(ν − 1) J 3 λ + e −λt 2 2e λt   + J 5 −2λ + e −λt 2e λt − e −λt 2 e λt ,

Proof. With the help of (20) and (24), we have for all



(33)

provided t ∈ [t 0 + 2N ε , t f − 2N ε ]. We notice that all the characteristic functions are equal to 1 and both sides of the previous

formula are independent on t. Hence, the mapping λ → P˜ ν (ε , λ) − Pν (λ) tends to 0 uniformly on any compact subset of C, as ε tends to 0 thanks to Lemma 4.1. Let us deal now with Qν . Expanding in Taylor series the exponentials ekλε w.r.t. ε , we find a matrix-valued convergent Taylor series w.r.t. ε for P˜ ν (ε , λ). The coefficient of εm in P˜ ν (ε , λ) is a polynomial matrix w.r.t. λ, independent of t inside [t 0 + 2N ε , t f − 2N ε ]. Now, the determinant of such a convergent Taylor series is itself a convergent Taylor series. At this point we have established that ζ 2Nd det(P˜ ν (ε , λ)) is a polynomial of degree 4Nd w.r.t. ζ = e λε and admits a Taylor expansion w.r.t. ε starting at det(Pν (λ)).

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

117

˜ ν . Let δ as in TheoWe choose ε0 so that Qν ⊂ K˙ ⊂ K ⊂ B (0, 1/ε0 ) and small enough to separate the elements of Q rem 4.1. We choose next  so that, if 1  m  4Nd, the coefficient of ζ m in det(P˜ ν ) − det(Pν ) is less than δ . Now, we may formulate the conclusion of Theorem 4.1 as the following inclusion

  −1   Q˜ ν = det P˜ ν (ε , .) {0} ⊂ C\B(0, 1/ε0 ) ∪



 B(λ, ε0 ) .

λ∈Qν

˜ ν ∩ K , Qν ) < As a consequence, intersecting both sides with K we get d H (Q proof. 2

ε0 for all ε small enough. This ends the

Remark 4.1. The convergence of x˜ (t ) to x(t ) as ε tends to 0 implies more complicated issues. Indeed, not only the phases Q˜ n and Q˜ 0 have to tend to Qn and Q0 respectively but the amplitudes x˜ s,α , x˜ j,α , and x˜ j,β , where α ∈ Q˜ n and β ∈ Q˜ 0 , have

also to tend to the respective amplitudes xs,α , x j ,α , and xj ,β , where of the difficulties in the case n = 1.

α ∈ Q n and β ∈ Q 0 . We refer to [13] for an examination

5. Periodicity and choreographies We focus in this section on periodic and choreographic solutions. In celestial mechanics, a choreography is a particular solution of Newton’s equations in which the n-bodies chase each other on the same closed curve, see for instance [3]. Although we do not work with Newton’s potential, we may extend this interesting concept for general equations of motion. We define a choreography of n particles (x1 , . . . , xn ) in Rd as a T -periodic solution to (2) or (3) for which the trajectories , k ∈ {1, . . . , n}. In other words, there exists a C 2 T -periodic mapping differ one to the other by some delay of the shape kT n u : R/ T Z → Rd such that, for all j, x j (t ) = u(t + jT /n) and

¨ t+ u

jT n







  T T = Fj u t + , . . . , u(t + T ), u˙ t + , . . . , u˙ (t + T ) ,



−2 2u t +

jT n



n



= F˜ j u t +

T n



n

  T , . . . , u(t + T ), 2u t + , . . . , 2u(t + T ) . n

We have used here the fact that the translation in C ([t 0 , t f ]) or in C p w commute with the derivative or the operator 2 respectively. 1

Theorem 5.1. 1. Under the assumptions of Proposition 3.1, all the solutions to C.E.L. are periodic if and only if

Qn ∪ Q0 ⊂ i R,

∀λ , λ ∈ Qn ∪ Q0 , λ /λ ∈ Q,

(34)

2. Under the assumptions of Proposition 3.3, all the pseudo-periodic solutions to D.E.L. are periodic if and only if

Q˜ n ∪ Q˜ 0 ⊂ i R,

∀λ , λ ∈ Q˜ n ∪ Q˜ 0 , λ /λ ∈ Q.

(35)

3. If det(Pn (0)) = 0, det(P˜ n (ε , 0)) = 0 and

|Q0 | = 2d, |Q˜ 0 | = 4Nd,

∀λ , λ ∈ Q0 , λ /λ ∈ Q,

Q0 ⊂ i R , Q˜ 0 ⊂ i R , 

∀λ , λ ∈ Q˜ 0 , λ /λ ∈ Q, 







(36) (37)

then there are choreographic solutions x j (t ) and x˜ j (t ) to C.E.L. and D.E.L. Proof. 1. We first notice that if {u }=1,..., K , is a family of nonzero vectors in Cd , then the various functions t → e λ t u are linearly independent if and only if the λ are pairwise distinct.  KIt relies on the nonsingularity of the Vandermonde matrix V (λ1 , . . . , λ K ). As a consequence, the function u(t ) = =1 e λ t u is periodic if and only if for some T > 0 we have ∀, λ T ∈ 2i π Z which is equivalent to the requirements ∀, i λ ∈ R and ∀ j , k, λ j /λk ∈ Q . Therefore, the

T λ  period T of u(t ) is inf{ T > 0, 2i π ∈ Z, ∀}. Taking into account that the vectors xs,α and x j ,β , occurring in the proof of Proposition 3.1, may be chosen arbitrarily in the respective appropriate null spaces ker(Pn (α )) and ker(P0 (β)), the previous properties of periodicity apply to the set of solutions to C.E.L. and give condition (34). 2. Pseudo-periodic solutions to D.E.L. are of the shape (4) by using Proposition 3.3 and the previous arguments apply.

118

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

3. We first observe that for each choreographic solution of the shape (4), xs (t ) is necessarily constant. Indeed,

xs (t ) =

n 

u t+

j =1

jT n

 = nu0 +

K 

=1

u e λ t

n  

e

λ T n

j

(38)

.

j =1

n

λ T

n ) j = 0. Having this fact in mind, we may solve (15). Our By periodicity, we have λ T ∈ 2i π Z for all  so that j =1 (e assumptions imply that the solution x j (t ) to (15) may be written as (4) with K = 2d, since the underlying Quadratic

jT

Eigenvalue Problem satisfies |Q0 | = 2d and det(Pn (0)) = 0 (see Proposition 3.1). Plugging x j (t ) = u(t + n ) into (15) and using the linear independence of the summands (4), we see that (15) is satisfied if and only if u0 = Pn (0)−1 J 7 and for all λ ∈ Q0 , uλ ∈ ker(P0 (λ)). If we choose the vectors x j (t 0 ) and x˙ j (t 0 ) or x j (t 0 ) and x j (t f ) for all j according to the preceding explicit form for x j (t ), we have justified the existence of choreographic solutions to C.E.L. with constant xs . Let us deal now with D.E.L. As seen in the proof of Proposition 3.3, each solution x˜ j (t ) to (17) has the shape (4) with

˜ 0 and P˜ n (ε , 0). The remainder of the proof is entirely similar to C.E.L. First, (17) K = 4Nd due to our assumptions on Q ˜ 0 , uλ ∈ ker(P˜ 0 (ε , λ)). Second, convenient choice is satisfied if and only if u0 = P˜ n (ε , 0)−1 ( J 7 − 2 J 6 ) and for all λ ∈ Q of initial or boundary conditions guarantees the existence of choreographic solutions to D.E.L. with constant xs and this ends the proof. 2 Remark 5.1. We may convert the existence of choreographic solutions into a linear algebra problem. Let us add to the systems described at the end of Sections 3.2 and 3.3 the following equations

xs,α = 0,

T

x j ,α = e α ( j −1) n x1,α

and

T

xj ,β = e β( j −1) n x1,β ,

for all α ∈ Qn , β ∈ Q0 and j = 1, . . . , n. Due to (22) and (23), provided the Dirichlet problem is well-posed, we find a choreographic solution. Remark 5.2. In literature (see for instance [14, Section 3.10]), the problem of the existence of periodic solutions arises when one studies gyroscopic systems. The algebraic conditions P7 and P8 in [14, Table 1.1] amount to requiring that J 1 = t J 1 > 0, J 2 = t J 2 > 0, J 3 and J 4 symmetric and small enough compared to J 1 and J 2 respectively. 6. Numerical experiments on choreographies Experimental algorithms performed in this last section are implemented in Maple and Matlab. We deal with real symmetric matrices J 1 , . . . , J 4 , zero vectors J 6 , J 7 , small dimension systems (d = 2, 3) since it already displays the main features, and arbitrary number of particles. Furthermore, the existence of periodic or choreographic solutions requires that J 5 = 0. Let us give  some details on the choice of the matrices J i . Given J 1 , J 2 , J 3 , we set, if d = 2, J4 =

1 2

j1 j2 j3 j4

J 2 + 12 ( J 1 − 2 J 3 )

. Identifying the coefficients of the polynomial det(P0 (λ)) with those of (λ2 + β12 )(λ2 + β22 )

and requiring that J 4 = t J 4 , we get three equations on j 1 , j 2 , j 3 , the coefficient j 4 standing free. Thus, we may choose J 1 + 2(n − 1) J 3 and J 1 − 2 J 3 definite positive, J 2 − 2 J 4 and J 2 + 2(n − 1) J 4 definite negative. We present in Figs. 1 and 2 the graphs of two solutions to gyroscopic C.E.L., sharing the same matrices J 1 , J 2 , J 3 . On the left, a typical periodic curve obtained √ by considering β1 = 4i and β2 = 10i and on the right, a non-periodic curve. Incommensurability between 4i and 7i 2 explains the non-choreographic behavior of the curve, as mentioned in property (34). Let us deal now with D.E.L. For the sake of clarity, we shall denote by x j (t ) and y j , M (t ), ∀ j ∈ {1, . . . , n}, the unique solution to C.E.L. (15) and the unique pseudo-periodic extension to [t 0 , t f ] of the unique solution to D.E.L. (29) on the grid

Gt0 ,ε with ε =

t f −t 0 M

respectively. First, we give some hints to solve (29). When t ∈ [t 0 + 2N ε , t f − 2N ε ], we may compute y j (t 0 + 2N ε ) as a function of / [t 0 + 2N ε , t f − 2N ε ], some of the characteristic functions occurring y j (t 0 + kε ) with k varying from −2N to 2N − 1. If t ∈ in (29) vanish and solving (29) must be slightly modified, see more details in [13]. We consider the matrices J 1 = J2 =



5 −1 −1 5



, J3 =



81 18





72 27



,

and (β1 , β2 ) = (2i , 5i ). In a first experiment, we use an operator 2 such that N = 1 and



1 1 (γ−1 , γ0 , γ1 ) = − + ik, −2ik, + ik where k ∈ R. 2

2

(39)

In both settings, C.E.L. as well as D.E.L., all the particles have the same trajectory since we have x j +1 (t ) = x j (t + n1 T ) and y j +1, M (t ) = y j , M (t + n1 T˜ ). Fig. 3 depicts the curves of y j ,35 , y j ,42 , y j ,75 and x j . Since our algorithms are suitable for each dimension, we provide also an example of quadratic choreography with d = 3 in Fig. 4. We choose J 1 , . . . , J 4 such that det(P˜ 0 (ε , λ)) = (λ2 + 1)(λ2 + 4)(λ2 + 9) and an operator of the shape (39).

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

119

Fig. 1. C.E.L., β1 = 4i and β2 = 10i.

√

Fig. 2. C.E.L., β1 = 4i and β2 = 7i 2.

Figs. 3 and 4 illustrate roughly the phenomenon of convergence, as M increases, of the pseudo-periodic solution y j , M (t ) to D.E.L. to the solution x j (t ) to C.E.L., for all j, for all operators 2 satisfying (39). For the approximation to become meaningful, the number M = #(Gt0 ,ε ) − 1 must satisfy

M2 

1

|γ − 1 γ 1 |

  (t f − t 0 )2 ρ ( J 1 − 2 J 3 )−1 ( J 2 − 2 J 4 )

as one sees from Eq. (29) (ρ denotes here the spectral radius). Let us remark that this lower bound is independent on d and still holds for general operators 2. In the next experiment, we still work with the previous Lagrangian, with N = 1 but with an operator 2 such that γ0 = −(γ−1 + γ1 ). We choose M = 100 in order to avoid some erratic behavior observed, for example, in the first plot (M = 35) of Fig. 3. We provide the plot of − log(min(x j (t ) − y j (t )2 , 3M )) as a function of (γ−1 , γ1 ) ∈ R2 in Fig. 5. Two peaks occur at (γ−1 , γ1 ) = ±(1/2, 1/2) and reveal a good approximation of x j (t ) by y j (t ). The two previous pairs are better understood if we have a look to the error x j (t ) − y j (t ) with complex operators 2 defined as in (39). In Fig. 6, we give the plot of the 2-norm of the error with (γ−1 , γ0 , γ1 ) = (− 12 + ik, −2ik, 12 + ik) for several values of M. The choice of k = 0 in (39) seems to be in any case the better one. Let us give an additional comment on the convergence of solutions, completing Remark 4.1. Let us first recall that D.E.L. converges to C.E.L. if and only if 2 is of the shape (1) and checks 21 = 0 and 2t = 1, provided t ∈ [t 0 + 2N ε , t f − 2N ε ], see [12, Definition 6.1. and Theorem 6.3]. In this case, the condition (γ−1 , γ0 , γ1 ) = (− 12 + ik, −2ik, 12 + ik) is linked to the

˜ 0 ⊂ i R as mentioned in [13, Proposition 5.1] for the special case d = 1. Based on the preceding experiments, we inclusion Q

120

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

Fig. 3. Solution of D.E.L. for M = 35, 42, 75 and solution to C.E.L. in R2 .

Fig. 4. Solution of D.E.L. for M = 30 and solution to C.E.L. in R3 .

conjecture that under mild condition of non-resonance of the Lagrangian, the solution to D.E.L. converges to the solution to C.E.L., as ε tends to 0. 7. Conclusion The least action principle leads to ODE in classical setting, to recurrence equations in the context of variational integrators [8,11], to asymptotic difference equations in the calculus of variations over Hölder spaces [4], or to a combination of recurrence and differential equations in the time-scales calculus [2]. Our work formalizes the least action principle for actions in which the derivatives of the dynamic variables have been replaced with generalized discrete operators, and leads instead to functional difference equations [12,13]. The context of quadratic Lagrangians for n-body problems is a convenient framework to fulfill all the steps from the physical model to its numerical treatment. Indeed, solving the underlying Euler–Lagrange equations may be performed through linear algebra techniques such as those of quadratic or transcendental eigenvalue problems. We focused specifically on the pseudo-periodic solutions and emphasized on the well-posedness of the equations of motion. We linked resonance phenomena [1] for the equations of motion to features of the generalized eigenvalue problems and also to convergence issues. In scope of future work, we plan to implement the D.E.L. scheme in the case of gravitic interaction of particles, and especially in the choreographic case [3].

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

Fig. 5. − log(min(x j (t ) − y j (t )2 , 3M )) as function of

121

γ−1 and γ1 (M = 100).

Fig. 6. 2-Norm of the error for several values of M.

Acknowledgements We thank the reviewers for their thorough review and appreciate the comments and suggestions, which contributed to improving the quality of the publication. References [1] [2] [3] [4] [5]

V.I. Arnold, Mathematical Methods of Classical Mechanics, Springer-Verlag, Berlin, 1989. M. Bohner, Calculus of variations on time scales, Dynam. Systems Appl. 13 (3–4) (2004) 339–349. A. Chenciner, J. Fejoz, R. Montgomery, Rotating eight I: the three Γi families, Nonlinearity 18 (3) (2004) 1407–1424. J. Cresson, Non-differentiable variational principles, J. Math. Anal. Appl. 307 (1) (2005) 48–64. J. Cresson, G.F.F. Frederico, D.F.M. Torres, Constants of motion for non-differentiable quantum variational problems, Topol. Methods Nonlinear Anal. 33 (2) (2009) 217–232.

122

[6] [7] [8] [9] [10] [11] [12] [13] [14]

P. Ryckelynck, L. Smoch / Applied Numerical Mathematics 75 (2014) 108–122

F. Cucker, A.G. Corbalan, An alternate proof of the continuity of the roots of a polynomial, Amer. Math. Monthly 96 (1989) 342–345. I. Gohberg, P. Lancaster, L. Rodman, Matrix Polynomials, Academic Press, New York, 1982. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration, second edition, Springer-Verlag, Berlin, 2006. P. Lancaster, Lambda-Matrices and Vibrating Systems, Pergamon Press, Oxford, UK, 1966. P. Lancaster, Quadratic eigenvalue problems, Linear Algebra Appl. 150 (1991) 499–506. J.E. Marsden, M. West, Discrete mechanics and variational integrators, Acta Numer. 10 (2001) 357–514. P. Ryckelynck, L. Smoch, Discrete calculus of variations for quadratic Lagrangians, Commun. Math. Anal. 15 (1) (2013) 44–60. P. Ryckelynck, L. Smoch, Discrete calculus of variations for oscillatory quadratic Lagrangians. Convergence issues, arXiv:1106.5350. F. Tisseur, K. Meerbergen, The quadratic eigenvalue problem, SIAM Review 43 (2001) 235–286.