Higher Order Variational Integrators for Multibody System Dynamics with Constraints

Abstract

The continuous and discrete Euler-Lagrangian equations with holonomic constraints are presented based on continuous and discrete Hamiltonian Principle. Using Lagrangian polynomial to interpolate state variables and Gauss quadrature formula to approximate Hamiltonian action integral, the higher order variational Galerkin integrators for multibody system dynamics with constraints and the computation procedure are given. Numerical examples are provided to show the long-time behavior of the methods proposed in this paper via comparisons with traditional Runge-Kutta methods.

1. Introduction

Dynamics of multibody systems [1] are usually described by a set of differential-algebraic equations (DAEs) [2, 3], which are solved traditionally using numerical methods of ordinary differential equations (ODEs) and combining some constraint stabilization techniques. In the community of numerical analysis, numerical methods for ODEs are always designed in a small area with local stability analysis, which lead to a lot of difficulties for long-time simulation. Geometric numerical integration methods [3] originated from long-term simulation of molecular dynamics and planet dynamics in solar system with global numerical stability due to the property of structure-preserving provide alternatives to overcome these problems. These numerical methods are referred to as the ones preserving the invariants in their corresponding continuous systems, such as preserving symplectic, energy, momentum, symmetry, and orthogonality, in the original systems, thus revolutionizing the traditional concepts of numerical stability. It has been widely accepted that the more invariants conserved in a numerical method, the more stable it is.

Symplectic algorithms [4] and energy conservation methods [5] are two frequently used methods for numerical integration of constrained Hamiltonian systems, but they are designed for conservative systems and it is difficult to design an algorithm that is symplectic and energy preserving at the same time [6]. The variational integrators [7, 8] based on discrete variational principle [9] provides an excellent framework for design of geometric numerical integrators with many merits: it preserves momentum naturally with good behavior of approximate energy conservation, even preserves strictly symplectic-momentum-energy [10]; it can be easily extended to a large class of problems [11, 12], such as the construction of geometric structure-preserving numerical integrators for PDEs [13], nonsmooth collisions [14, 15], stochastic systems [16, 17], nonholonomic systems [18, 19], and constrained systems [20], dissipative systems [21], optimal control [22, 23], parameter optimization [24]; it can generate easily a large class of higher order methods based on polynomial interpolation and numerical quadrature formulas systematically [25 –29] and some Lie group variational integrators [30 –32] combining Lie group methods [33].

For the general conservative systems without constraints, [25, 26] propose the same framework of design higher order variational Galerkin integrators using numerical interpolation and quadrature techniques for state variables and action integral, respectively, during one integration step based on discrete variational principle. Motivated by these works, we will investigate higher order variational Galerkin integrators for multibody system dynamics with constraints for the long-term simulation purpose. Due to the space limit, we will focus on the cases with holonomic constraints, but the results can be easily extended to the cases with nonholonomic constraints and the cases with nonconservative forces.

The paper is organized as follows. In Section 2, the continuous and discrete Euler-Lagrangian equations with holonomic constraints are presented based on continuous and discrete Hamiltonian Principle, respectively; in Section 3, we derive the higher order variational Galerkin integrators for multibody system dynamics with constraints using Lagrangian polynomial to interpolate state variables and using Gauss quadrature formula to approximate Hamiltonian action integral; the computation procedure is given in Section 4; numerical examples are provided in Section 5 to show the long-time behavior of the methods proposed in this paper via comparisons with traditional methods. The last one is concluding remarks including summary and future works.

2. Discrete Euler-Lagrangian Equations with Constraints

The Hamilton Principle of a multibody system with holonomic constraints and conservative forces can be stated as

δ S = δ \int_{0}^{t_{f}} (L (q, \dot{q}, t) - λ^{T} Φ (q, t)) d t = 0,

(1)

where S is Hamilton action integral, L is Lagrangian, q, $\dot{q}$ are generalized displacement and velocity, respectively, and λ is Lagrange multiplier vector corresponding to the constraint

Φ (q, t) = 0 .

(2)

Via variational method, the following Euler-Lagrange equations of a constrained mechanical system can be got:

\begin{matrix} \frac{d}{d t} (\frac{\partial L (q, \dot{q}, t)}{\partial \dot{q}}) - \frac{\partial L (q, \dot{q}, t)}{\partial q} + Φ^{T} (q, t) λ = 0, \\ Φ (q, t) = 0 . \end{matrix}

(3)

It is a typical differential-algebraic equation which can be solved using different traditional numerical methods of DAEs [2, 3].

Following the variational integrators based on discrete variational principle, the time interval t∈[0,t_f] is divided into N subintervals with time step h = t_f/N equally along with N + 1 grids, 0 = t₀,t₁,t₂,…, t_N–1, t_N = t_f, t_i = ih, i = 0,1,2,…, N. The discrete generalized displacement and velocity and Lagrangian multiplier can be denoted as q_i = q(t_i), ${\dot{q}}_{i} = \dot{q} (t_{i})$ , λ_i = λ(t_i), i = 0,1,2,…, N. Equation (1) can be rewritten as the following discrete counterpart:

\begin{matrix} δ S_{d} (q_{0}, q_{1}, \dots, q_{N - 1}, q_{N}, λ_{0}, \dots, λ_{N - 2}, λ_{N - 1}) \\ = δ \sum_{i = 0}^{N - 1} L_{d} (q_{i}, q_{i + 1}, t_{i}) - δ \sum_{i = 0}^{N - 1} f_{d} (q_{i}, q_{i + 1}, λ_{i}, t_{i}) \\ = 0, \end{matrix}

(4)

where L_d(q_i,q_{i + 1},h) and f_d(q_i,q_{i + 1},λ_i,t_i) are defined as the discrete approximation of the integration of $L (q, \dot{q}, t)$ and λ^TΦ(q,t) on the interval t∈[t_i,t_{i + 1}],

\begin{matrix} L_{d} (q_{i}, q_{i + 1}, t_{i}) = \int_{t_{i}}^{t_{i + 1}} L (q, \dot{q}, t) d t, \\ f_{d} (q_{i}, q_{i + 1}, λ_{i}, t_{i}) = \int_{t_{i}}^{t_{i + 1}} λ^{T} Φ (q, t) d t . \end{matrix}

(5)

Using standard variational method, the discrete Euler-Lagrangian (DEL) equations are derived as

\begin{matrix} D_{1} L_{d} (q_{i}, q_{i + 1}, t_{i}) + D_{2} L_{d} (q_{i - 1}, q_{i}, t_{i - 1}) - h Φ_{q}^{T} (q_{i}, t_{i}) λ_{i} = 0, \\ Φ (q_{i + 1}, t_{i + 1}) = 0, \end{matrix}

(6)

where D_jL_d (j = 1, 2) is the partial derivative of L_d with respect to the jth variable.

3. Higher Order Variational Galerkin Integrators of Multibody Systems with Constraints

The key technique to improve the accuracy of variational integrators is the approximation of the discrete Hamiltonian action integral, which can be realized through general Galerkin methods [24, 25]. Using the same methods in this section, Hamiltonian action integral is approximated in a subinterval $[t_{i}, t_{i + 1}]$ with an introduced local time parameter $τ \in [0,1]$ , where

τ = \frac{t - t_{i}}{t_{i + 1} - t_{i}} .

(7)

Given s + 1 control points 0 = τ₀<τ₁<⋯<τ_s = 1, for k = 0, 1,…, s, Lagrange interpolation polynomial l_k(τ):[0, 1]→R is defined as

l_{k} (τ) = \prod_{0 \leq j \leq s, j \neq k} \frac{τ - τ_{j}}{τ_{k} - τ_{j}} .

(8)

The state variables q(t) on $[t_{i}, t_{i + 1}]$ can be interpolated as follows:

q_{d} (τ) = \sum_{k = 0}^{s} l_{k} (τ) q_{i, k},

(9)

where $q_{i, k} = q (t_{i} + k \cdot (h / s))$ . The derivative of q_d(τ) with respect to t provides an approximation of $\dot{q} (t)$ on $[t_{i}, t_{i + 1}]$ as

{\dot{q}}_{d} (τ) = \frac{1}{h} \sum_{k = 0}^{s} l_{k}^{'} (τ) q_{i, k},

(10)

where $l_{k}^{'} (τ) = (\partial / \partial τ) l_{k} (τ)$ . Then, Hamiltonian action integral can be approximated on $[t_{i}, t_{i + 1}]$ as follows:

\begin{array}{l} \int_{t_{i}}^{t_{i + 1}} (L (q (t), \dot{q} (t), t) - λ^{T} Φ (q, t)) d t \\ = \int_{0}^{1} h (L (q_{d} (τ), {\dot{q}}_{d} (τ), t_{i} + τ h) - λ_{i}^{T} Φ (q_{d} (τ), t_{i} + τ h)) d τ . \end{array}

(11)

Lagrange interpolation polynomials in (9) and comparable derivatives with respect to τ in (10) are listed in Table 1.

Table 1:

Lagrange interpolation polynomials and comparable derivatives with respect to τ.

s	l_k(τ)	l_k′(τ)

1	$\begin{matrix} l_{0} (τ) = 1 - τ \\ l_{1} (τ) = τ \end{matrix}$	$\begin{matrix} l_{0}^{'} (τ) = - 1 \\ l_{1}^{'} (τ) = 1 \end{matrix}$
2	$\begin{matrix} l_{0} (τ) = (2 τ - 1) (τ - 1) \\ l_{1} (τ) = - 4 τ (τ - 1) \\ l_{2} (τ) = τ (2 τ - 1) \end{matrix}$	$\begin{matrix} l_{0}^{'} (τ) = 4 τ - 3 \\ l_{1}^{'} (τ) = - 4 (2 τ - 1) \\ l_{2}^{'} (τ) = 4 τ - 1 \end{matrix}$

Using Gauss quadrature formula, which has maximal degree of accuracy 2n–1 for a fixed number n of quadrature points, we have

\begin{array}{l} \int_{t_{i}}^{t_{i + 1}} (L (q (t), \dot{q} (t), t) - λ^{T} Φ (q, t)) d t \\ = \sum_{r = 1}^{n} A_{r} h (L (q_{d} ({\tilde{τ}}_{r}), {\dot{q}}_{d} ({\tilde{τ}}_{r}), t_{i} + {\tilde{τ}}_{r} h) - λ_{i}^{T} Φ (q_{d} ({\tilde{τ}}_{r}), t_{i} + {\tilde{τ}}_{r} h)) \\ = L_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, t_{i}) - f_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, λ_{i}, t_{i}), \end{array}

(12)

where A_r is the weights and ${\tilde{τ}}_{r}$ is the quadrature point. Quadrature points and weights of the Gauss quadrature formula with respect to the interval $[- 1, 1]$ are listed in Table 2.

Table 2:

Quadrature points and weights of the Gauss quadrature formula with respect to the interval [−1, 1].

n	${\tilde{τ}}_{r}$	A _r

2	$\pm \sqrt{3} / 3$	1
3	$\pm \sqrt{15} / 5$	5/9
3	0	8/9

Apply discrete Hamilton's principle, the derivatives of the action with respect to q_i,k have to vanish for all i = 0, 1,…, N–1 and k = 0, 1,…, s. That leads to the discrete Euler-Lagrange equations

\begin{matrix} D_{1} L_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, t_{i}) - D_{1} f_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, λ_{i}, t_{i}) + D_{s + 1} L_{d} (q_{i - 1,0}, q_{i - 1,1}, \dots, q_{i - 1, s}, t_{i - 1}) - D_{s + 1} f_{d} (q_{i - 1,0}, q_{i - 1,1}, \dots, q_{i - 1, s}, λ_{i - 1}, t_{i - 1}) = 0 \\ D_{k} L_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, t_{i}) - D_{k} f_{d} (q_{i, 0}, q_{i, 1}, \dots, q_{i, s}, λ_{i}, t_{i}) = 0, k = 1, \dots, s \\ Φ (q_{i + 1}, t_{i + 1}) = 0 . \end{matrix}

(13)

4. Computation Procedure

The following computation procedure can be used to get the solution of the discrete Euler-Lagrange equations (13).

Step 1. Given the order of Lagrange interpolation polynomial s, the number of Gauss quadrature points n, fixed time step h, initial state variables q₀, and velocity variables ${\dot{q}}_{0}$ .

Step 2. Get q₁, ${\dot{q}}_{1}$ , and λ₁ by Euler method or Runge-Kutta method. Use q_0,0 = q₀, and q_0,s = q₁ to calculate q_0,k, k = 1,…, s–1.

Step 3. For i = 2 to N–1, solve (13) by Newton iteration method to get q_i,1,…,q_i,s and λi, thus q_{i + 1} = q_i,s, ${\dot{q}}_{i + 1} = {\dot{q}}_{d} (1) = (1 / h) \sum_{k = 0}^{s} l_{k}^{'} (1) q_{i, k}$ .

Step 4. Use the solution qi, ${\dot{q}}_{i + 1}$ , i = 0,1,…, N–1 to get the potential energy T, kinetic energy V, total energy of the system H = T + V, and constraints Φ.

5. Numerical Example

Consider a simple double pendulum system shown in Figure 1 the state variables are chosen as q = [x₁ y₁ x₂ y₂]^T, then the potential energy, kinetic energy of the system are

\begin{matrix} T = \frac{1}{2} m_{1} ({\dot{x}}_{1}^{2} + {\dot{y}}_{1}^{2}) + \frac{1}{2} m_{2} ({\dot{x}}_{2}^{2} + {\dot{y}}_{2}^{2}), \\ V = m_{1} g y_{1} + m_{2} g y_{2} . \end{matrix}

(14)

The mass matrix is a constant matrix as follows:

M = diag [\begin{bmatrix} m_{1} & m_{1} & m_{2} & m_{2} \end{bmatrix}] .

(15)

The constraint equations are

Φ (q) = (\begin{pmatrix} x_{1}^{2} + y_{1}^{2} - l_{1}^{2} \\ {(x_{2} - x_{1})}^{2} + {(y_{2} - y_{1})}^{2} - l_{2}^{2} \end{pmatrix}) = 0 .

(16)

Figure 1:

Double pendulum.

Given m₁ = m₂ = 1 kg, l₁ = l₂ = 1 m, with time step h = 0.01's, the results of the method presented in this paper are compared with the Runge-Kutta method in Table 3 for the terminal time 100's. Here, P1L-Q2G means that the order of Lagrange interpolation polynomial is 1 and the number of Gauss quadrature points is 2, and the same as P2L-Q2G and P2L-Q3G. ε(Η), ε(Φ) are the maximum errors of the total energy H and constraints Φ.

Table 3:

Comparison of the methods mentioned with time step h = 0.01.

Numerical method	Run time (s)	$\max_{} \| ε (H) \|$	$\max_{} \| ε (Φ) \|$

Runge-Kutta	15.0385	3.0045	0.2345
P1L-Q2G	16.0993	0.4049	1.9984 × 10⁻¹⁵
P2L-Q2G	24.8510	0.2270	8.8818 × 10⁻¹⁶
P2L-Q3G	31.2938	0.2273	1.3323 × 10⁻¹⁵

Using the traditional method such as Runge-Kutta method, to keep the errors of the total energy H and constraints Φ are all smaller, the time step should be chosen as about 0.001's in this example, and the total computation time is 99.7938's. But using the discrete variational methods presented in this paper, the time step can be chosen longer as 0.01's. Table 3 shows that, with this longer time step, the discrete variational methods can keep the total energy and constraints in higher accuracy than Runge-Kutta method, the cost instead is the longer computation time because of the iterations during the nonlinear equations solving process, but it is still shorter than the computation time of Runge-Kutta method with smaller time step to reach the same accuracy of the results.

It is also shown in Tables 3 and 4 that when the order of Lagrange interpolation polynomial is higher, the total energy and constraints can be kept better, but there are no big differences between 2 or 3 points in Gauss quadrature formula.

Table 4:

Comparison of the methods mentioned with time step h = 0.001.

Numerical method	Run time (s)	$\max_{} \| ε (H) \|$	$\max_{} \| ε (Φ) \|$

Runge-Kutta	99.7938	4.7378 × 10⁻⁴	3.6633 × 10⁻⁵
P1L-Q2G	121.4468	0.0405	9.6040 × 10⁻¹¹
P2L-Q2G	212.9882	0.0023	1.2912 × 10⁻⁹
P2L-Q3G	265.9505	0.0021	1.2912 × 10⁻⁹

Figure 2 shows the energies of the previously mentioned methods. With the time step h = 0.01's, the total energy is up and down around zero with little errors during the long-time simulation by the methods presented in this paper, while increases quickly during the simulation by Runge-Kutta method. The results of high order variational Galerkin integrators are better than low order integrators.

Figure 2:

Total energy of double pendulum by different method, h = 0.01.

In the case of nonconstant mass matrix, the methods of higher order variational Galerkin integrators are also applicable. For another example, Figure 3 shows a rotary rod slider system. OA is a rigid rod with uniform mass m_A and length l, which rotates round O in the plane OXY. B is a slider with the mass m_B, and the stiffness of the spring on it is k; the mass of the spring is ignored. Suppose only gravity in the plane OXY is considered for the system.

Figure 3:

A rotary rod slider system.

The state variables are chosen as q = [R_{x
ⁱ} R_{y
ⁱ} R_{x
^j} R_{y
^j}θ x_i]^T. Then, the mass matrix is a nonconstant matrix

M = [\begin{bmatrix} m_{B} & 0 & 0 & 0 & - m_{B} x_{i} \sin θ & m_{B} \cos θ \\ 0 & m_{B} & 0 & 0 & m_{B} x_{i} \cos θ & m_{B} \sin θ \\ 0 & 0 & m_{A} & 0 & 0 & 0 \\ 0 & 0 & 0 & m_{A} & 0 & 0 \\ - m_{B} x_{i} \sin θ & m_{B} x_{i} \cos θ & 0 & 0 & \frac{1}{3} m_{A} l^{2} + m_{B} x_{i}^{2} & 0 \\ m_{B} \cos θ & m_{B} \sin θ & 0 & 0 & 0 & m_{B} \end{bmatrix}] .

(17)

The constraint equations are

Φ = (\begin{pmatrix} R_{x^{j}} - \frac{l}{2} \cos θ \\ R_{y^{j}} - \frac{l}{2} \sin θ \\ (R_{x^{j}} - R_{x^{i}}) \cos θ + (R_{y^{j}} - R_{y^{i}}) \sin θ - x_{i} \\ (R_{x^{j}} - R_{x^{i}}) \sin θ - (R_{y^{j}} - R_{y^{i}}) \cos θ \end{pmatrix}) = 0 .

(18)

Given l₁ = 2m, m_A = 2 kg, m_B = 5 kg, k = 1000 N/m. With the initial step length h = 0.01's and the terminal time 100's, using the method P1L-Q2G, where the order of Lagrange interpolation polynomial is 1 and the number of Gauss quadrature points is 2, the results are also satisfying. The total computer time is 86.0034's; the maximum errors of the total energy H and constraints Φ are 0.2611 and 3.8858 × 10⁻¹⁶, which is illustrated in Figures 4 and 5.

Figure 4:

Energies of rotary rod slider system by method P1L-Q2G, h = 0.01.

Figure 5:

Constraints of rotary rod slider system by method P1L-Q2G, h = 0.01.

6. Conclusions

Based on Lagrange interpolation and Gauss quadrature, higher order variational Galerkin integrators have been developed for multibody systems with constraints. Numerical results validate the accuracy of the method for longer time step and long-time simulation. Nevertheless, as the computation is more complex, the computation cost is higher. Further study is needed to have more efficient integrators.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Footnotes

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grants nos. 11272166, 11002075, 10972110).

References

Shabana

A. A.

, Dynamics of Multibody Systems, Cambridge University Press, 3rd edition, 2005.

Laulusa

and Bauchau

O. A.

, “Review of classical approaches for constraint enforcement in multibody systems,” Journal of Computational and Nonlinear Dynamics, vol. 3, no. 1, Article ID 011004, 2008.

Bauchau

O. A.

and Laulusa

, “Review of contemporary approaches for constraint enforcement in multibody systems,” Journal of Computational and Nonlinear Dynamics, vol. 3, no. 1, Article ID 011005, 2008.

Jay

, “Symplectic partitioned Runge-Kutta methods for constrained Hamiltonian systems,” SIAM Journal on Numerical Analysis, vol. 33, no. 1, pp. 368–387, 1996.

Gonzalez

, “Mechanical systems subject to holonomic constraints: differential-algebraic formulations and conservative integration,” Physica D, vol. 132, no. 1–2, pp. 165–174, 1999.

and Marsden

J. E.

, “Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators,” Physics Letters A, vol. 133, no. 3, pp. 134–139, 1988.

Wendlandt

J. M.

and Marsden

J. E.

, “Mechanical integrators derived from a discrete variational principle,” Physica D, vol. 106, no. 3–4, pp. 223–246, 1997.

Marsden

J. E.

and West

, “Discrete Mechanics and Variational Integrators,” Acta Numerica, vol. 10, pp. 357–514, 2001.

Moser

and Veselov

A. P.

, “Discrete versions of some classical integrable systems and factorization of matrix polynomials,” Communications in Mathematical Physics, vol. 139, no. 2, pp. 217–243, 1991.

10.

Kane

Marsden

J. E.

, and Ortiz

, “Symplectic-energy-momentum preserving variational integrators,” Journal of Mathematical Physics, vol. 40, no. 7, pp. 3353–3371, 1999.

11.

Lew

Marsden

J. E.

Ortiz

, and West

, “Variational time integrators,” International Journal for Numerical Methods in Engineering, vol. 60, no. 1, pp. 153–212, 2004.

12.

Lew

Marsden

J. E.

Ortiz

, and West

, “An overview of variational integrators,” in Finite Element Methods: 1970s and Beyond, pp. 85–146, CIMNE, Barcelona, Spain, 2003.

13.

Lew

Marsden

J. E.

Ortiz

, and West

, “Asynchronous variational integrators,” Archive for Rational Mechanics and Analysis, vol. 167, no. 2, pp. 85–146, 2003.

14.

Pandolfi

Kane

Marsden

J. E.

, and Ortiz

, “Time-discretized variational formulation of non-smooth frictional contact,” International Journal for Numerical Methods in Engineering, vol. 53, no. 8, pp. 1801–1829, 2002.

15.

Fetecau

R. C.

Marsden

J. E.

Ortiz

, and West

, “Nonsmooth lagrangian mechanics and variational collision integrators,” SIAM Journal on Applied Dynamical Systems, vol. 2, no. 3, pp. 381–416, 2003.

16.

Bou-Rabee

and Owhadi

, “Stochastic variational integrators,” IMA Journal of Numerical Analysis, vol. 29, no. 2, pp. 421–443, 2009.

17.

Deng

, “Generalized multi-symplectic integrators for a class of Hamiltonian nonlinear wave PDEs,” Journal of Computational Physics, vol. 235, pp. 394–406, 2013.

18.

Cortés

and Martínez

, “Non-holonomic integrators,” Nonlin-earity, vol. 14, no. 5, pp. 1365–1392, 2001.

19.

McLachlan

and Perlmutter

, “Integrators for nonholonomic mechanical systems,” Journal of Nonlinear Science, vol. 16, no. 4, pp. 283–328, 2006.

20.

Leyendecker

Marsden

J. E.

, and Ortiz

, “Variational integrators for constrained dynamical systems,” ZAMM Zeitschrift für Angewandte Mathematik und Mechanik, vol. 88, no. 9, pp. 677–708, 2008.

21.

Kane

Marsden

J. E.

Ortiz

, and West

, “Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems,” International Journal for Numerical Methods in Engineering, vol. 49, no. 10, pp. 1295–1325, 2000.

22.

Leyendecker

Ober-Blöbaum

Marsden

J. E.

, and Ortiz

, “Discrete mechanics and optimal control for constrained systems,” Optimal Control Applications and Methods, vol. 31, no. 6, pp. 505–528, 2010.

23.

Ober-Blöbaum

Junge

, and Marsden

J. E.

, “Discrete mechanics and optimal control: an analysis,” Control, Optimisation and Calculus of Variations, vol. 17, no. 2, pp. 322–352, 2011.

24.

Witting

Ober-Blöbaum

, and Dellnitz

, “A variational approach to define robustness for parametric multi-objective optimization problems,” Journal of Global Optimization, vol. 57, no. 2, pp. 331–345, 2013.

25.

Leok

and Zhang

, “Discrete Hamiltonian variational integrators,” IMA Journal of Numerical Analysis, vol. 31, no. 4, pp. 1497–1532, 2011.

26.

Ober-Blöbaum

and Saake

, “Construction and analysis of higher order Galerkin variational integrators,” Numerical Analysis. In press. http://arxiv.org/abs/1304.1398.

27.

Leok

Shingel

, and Shingel

T. P.

, “rolongation-collocation variational integrators,” IMA Journal of Numerical Analysis, vol. 32, no. 3, pp. 1194–1216, 2012.

28.

Leyendecker

and Ober-Blöbaum

, “A variational approach to multirate integration for constrained systems,” in Multibody Dynamics, Samin

J.-C.

and Fisette

, Eds., vol. 28 of Computational Methods in Applied Sciences, pp. 97–121, Springer, 2013.

29.

Campos

C. M.

Junge

, and Ober-Blöbaum

, “Higher order variational time discretization of optimal control problems,” in Proceedings of the 20th International Symposium on Mathematical Theory of Networks and Systems, pp. 9–13, Melbourne, Australia, July 2012.

30.

Leok

and Shingel

, “General techniques for constructing variational integrators,” Frontiers of Mathematics in China, vol. 7, no. 2, pp. 273–303, 2012.

31.

Bou-Rabee

and Marsden

J. E.

, “Hamilton-Pontryagin integrators on lie groups part I: introduction and structure-preserving properties,” Foundations of Computational Mathematics, vol. 9, no. 2, pp. 197–219, 2009.

32.

Lee

Leok

, and McClamroch

N. H.

, “Lie group variational integrators for the full body problem,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 29–30, pp. 2907–2924, 2007.

33.

Iserles

Munthe-Kaas

H. Z.

Nørsett

S. P.

, and Zanna

, “Lie-Group methods,” Acta Numerica, vol. 9, pp. 215–365, 2000.