# Simpson's Variational Integrator for Systems with Quadratic Lagrangians

Juan Antonio Rojas-Quintero, François Dubois, José Guadalupe Cabrera-Díaz

► To cite this version:

Juan Antonio Rojas-Quintero, François Dubois, José Guadalupe Cabrera-Díaz. Simpson's Variational Integrator for Systems with Quadratic Lagrangians. Axioms, 2024, 13 (4), pp.255. 10.3390/axioms13040255. [https://hal.science/hal-04548211](https://hal.science/hal-04548211)

HAL Id: hal-04548211

https://hal.science/hal-04548211

Submitted on 16 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SIMPSON'S VARIATIONAL INTEGRATOR FOR SYSTEMS WITH QUADRATIC
LAGRANGIANS

JUAN ANTONIO ROJAS-QUINTERO 1,2,*, FRANÇOIS DUBOIS 3,4, AND JOSÉ GUADALUPE CABRERA-DÍAZ 5

ABSTRACT. This contribution proposes a variational symplectic integrator aimed at linear sys-
tems issued from the least action principle. An internal quadratic finite elements interpolation of
the state is performed at each time step. Then, the action is approximated by Simpson's quadrature
formula. The implemented scheme is implicit, symplectic, and conditionally stable. It is applied to
the time integration of systems with quadratic Lagrangians. The example of the linearized double
pendulum is treated. Our method is compared with Newmark's variational integrator. The exact
solution of the linearized double pendulum example is used for benchmarking. Simulation results
illustrate the precision and convergence of the proposed integrator.

**Keywords.** Ordinary differential equations; oscillator; numerical analysis; symplectic scheme.
MSC. 34A30, 65L05, 65P10.

1. INTRODUCTION

*Simpson's quadrature* is the name that is generally given to a numerical approximation of
definite integrals that is exact for polynomials up to the third degree:

$$
(1) \qquad \int_0^1 \psi(\theta) d\theta \approx \frac{1}{6} \left( \psi(0) + 4\psi\left(\frac{1}{2}\right) + \psi(1) \right).
$$

It is well-known that this rule was found by Bonaventura Cavalieri (1598–1647), known to James Gregory (1638–1675) [1], and even used by Johannes Kepler (1571–1630) to approximate the volume of barrels [2]. However, Thomas Simpson (1710–1761) is usually credited for this rule. As such, formula (1) is also widely known as *Simpson's 1/3 rule*. It corresponds to a special case of Newton-Cotes's formula [1] and coincides with the classical fourth-order Runge-Kutta method [1, 3].

Generally, numerical methods involving Simpson's quadrature estimate a definite integral
by using quadratic polynomials to approximate the integrand on a sequence of intervals. This
general idea is at the foundation of numerous methods that can be applied to solve engineering
problems such as the low-thrust orbit transfer problem [4], or the gait optimization of a bipedal
walking robot [5]. Recently, much attention has been brought to fractional calculus, for which
solvers based on Simpson's quadrature (adapted to fractional form) have been developed [6].
Some applications involve solving initial value problems of fractional differential systems [7]
or the solution of fractional equations affected by noisy signals [8]. Another recent application
of Simpson's quadrature involves the solution of partial integro-differential equations [9].

Our contribution is aimed at solving differential equations characterizing the motions of
mechanical systems. It is well known that the motions of a mechanical system are the extremals
of the variational principle of least action [10]. This principle is one of the most general laws of
theoretical physics and is foundational for characterizing a system evolution in the form of
differential equations. It is valid across disciplines such as classical and quantum mechanics,
cosmology, electromagnetism, optics, and relativity [10–14]. As such, this variational principle
is closely involved in the development of the finite-element method [15], which is used for
space and time integration of differential equations [16].

<footnote n="1">1 CONAHCYT/Tecnológico Nacional de México/I. T. Ensenada, Ensenada, 22780, B.C., Mexico</footnote>

<footnote n="2">2 IMT Atlantique, LS2N, UMR CNRS 6004, 44307 Nantes, France</footnote>

<footnote n="3">³ Université Paris-Saclay, Laboratoire de Mathématiques d'Orsay, Orsay, 91400, France</footnote>

<footnote n="4">⁴ Conservatoire National des Arts et Métiers, LMSSC, Paris, 75141, France</footnote>

<footnote n="5">⁵ Tecnológico Nacional de México/I. T. Ensenada, Ensenada, 22780, B.C., Mexico</footnote>

<footnote>E-mail addresses: * jarojas@conahcyt.mx.</footnote>

<footnote>Date: April 5, 2024 (author version). Published in *Axioms*. 2024; 13(4):255. doi:10.3390/axioms13040255.</footnote>
SIMPSON'S VARIATIONAL INTEGRATOR FOR SYSTEMS WITH QUADRATIC
LAGRANGIANS

JUAN ANTONIO ROJAS-QUINTERO 1,2,*, FRANÇOIS DUBOIS 3,4, AND JOSÉ GUADALUPE CABRERA-DÍAZ 5

ABSTRACT. This contribution proposes a variational symplectic integrator aimed at linear sys-
tems issued from the least action principle. An internal quadratic finite elements interpolation of
the state is performed at each time step. Then, the action is approximated by Simpson's quadrature
formula. The implemented scheme is implicit, symplectic, and conditionally stable. It is applied to
the time integration of systems with quadratic Lagrangians. The example of the linearized double
pendulum is treated. Our method is compared with Newmark's variational integrator. The exact
solution of the linearized double pendulum example is used for benchmarking. Simulation results
illustrate the precision and convergence of the proposed integrator.

**Keywords.** Ordinary differential equations; oscillator; numerical analysis; symplectic scheme.
MSC. 34A30, 65L05, 65P10.

1. INTRODUCTION

*Simpson's quadrature* is the name that is generally given to a numerical approximation of
definite integrals that is exact for polynomials up to the third degree:

$$
(1) \qquad \int_0^1 \psi(\theta) d\theta \approx \frac{1}{6} \left( \psi(0) + 4\psi\left(\frac{1}{2}\right) + \psi(1) \right).
$$

It is well-known that this rule was found by Bonaventura Cavalieri (1598–1647), known to James Gregory (1638–1675) [1], and even used by Johannes Kepler (1571–1630) to approximate the volume of barrels [2]. However, Thomas Simpson (1710–1761) is usually credited for this rule. As such, formula (1) is also widely known as *Simpson's 1/3 rule*. It corresponds to a special case of Newton-Cotes's formula [1] and coincides with the classical fourth-order Runge-Kutta method [1, 3].

Generally, numerical methods involving Simpson's quadrature estimate a definite integral
by using quadratic polynomials to approximate the integrand on a sequence of intervals. This
general idea is at the foundation of numerous methods that can be applied to solve engineering
problems such as the low-thrust orbit transfer problem [4], or the gait optimization of a bipedal
walking robot [5]. Recently, much attention has been brought to fractional calculus, for which
solvers based on Simpson's quadrature (adapted to fractional form) have been developed [6].
Some applications involve solving initial value problems of fractional differential systems [7]
or the solution of fractional equations affected by noisy signals [8]. Another recent application
of Simpson's quadrature involves the solution of partial integro-differential equations [9].

Our contribution is aimed at solving differential equations characterizing the motions of
mechanical systems. It is well known that the motions of a mechanical system are the extremals
of the variational principle of least action [10]. This principle is one of the most general laws of
theoretical physics and is foundational for characterizing a system evolution in the form of
differential equations. It is valid across disciplines such as classical and quantum mechanics,
cosmology, electromagnetism, optics, and relativity [10–14]. As such, this variational principle
is closely involved in the development of the finite-element method [15], which is used for
space and time integration of differential equations [16].

<footnote n="1">1 CONAHCYT/Tecnológico Nacional de México/I. T. Ensenada, Ensenada, 22780, B.C., Mexico</footnote>

<footnote n="2">2 IMT Atlantique, LS2N, UMR CNRS 6004, 44307 Nantes, France</footnote>

<footnote n="3">³ Université Paris-Saclay, Laboratoire de Mathématiques d'Orsay, Orsay, 91400, France</footnote>

<footnote n="4">⁴ Conservatoire National des Arts et Métiers, LMSSC, Paris, 75141, France</footnote>

<footnote n="5">⁵ Tecnológico Nacional de México/I. T. Ensenada, Ensenada, 22780, B.C., Mexico</footnote>

<footnote>E-mail addresses: * jarojas@conahcyt.mx.</footnote>

<footnote>Date: April 5, 2024 (author version). Published in *Axioms*. 2024; 13(4):255. doi:10.3390/axioms13040255.</footnote>
Numerical schemes for dynamical systems issued from the principle of least action are typically referred to as variational [17–20]. The general idea resides in performing a discretization at the least action principle level. As a result, the evolution equations deriving from this discretized principle characterize the system evolution but are also a numerical scheme. It is well known that such numerical methods are endowed with interesting characteristics; one characteristic is the property of being symplectic [18–21]. One remarkable example of such methods is Newmark’s integrator [17, 19], which is very popular for solving problems in dynamics of structures [22, 23], and has recently been geometrized to solve motion equations of sliding rods [24] and soft robots [25].

A symplectic scheme based on Simpson’s rule has been proposed by the authors in [26], for the linear and scalar case of the harmonic oscillator. The scheme uses a quadratic finite elements interpolation. The method was adapted to the monodimensional non-linear pendulum system in [27]. In this work, Simpson’s symplectic scheme is further studied as an alternative to Newmark’s method. It is generalized to the case of multiple degrees of freedom systems characterized by quadratic Lagrangians. Obtained results confirm the convergence rate previously observed in [26]. The new stability condition on the step size is revealed to be similar to the one previously obtained in [26]. A simplecticity analysis that applies to the multi-degrees-of-freedom case, along with the expression of a related conserved quadratic form is provided in this contribution.

We begin our study by detailing Newmark’s classical scheme, deriving it from variational principles in section 2. Then, Simpson’s alternative scheme is detailed and derived from variational principles in section 3. Section 4 analyzes the symplectic property of Simpson’s scheme. A proof that applies to both Newmark’s and Simpson’s schemes is provided. To compare both methods in a case study, a two-degrees of freedom system is presented. Therefore, the exact solution to the linearized double pendulum is provided in section 5. This exact solution serves for benchmarking purposes in our comparisons. Section 6 presents and comments on the obtained numerical results. Simpson’s scheme convergence is revealed to be of the fourth order. The manuscript ends with a brief discussion and concluding remarks in section 7.

## 2. NEWMARK’S SCHEME

### 2.1. Discrete action

Let us derive the classical, symplectic variational integrator based on Newmark’s scheme [17, 19, 22, 23]. The continuous action is defined as

$$ (2) \qquad S_c = \int_0^T \mathcal{L} \left( \frac{d\mathbf{q}(t)}{dt}, \mathbf{q}(t) \right) dt $$

where $\mathcal{L}$ is the system Lagrangian. We focus on dynamical systems for which the Lagrangian can be expressed quadratically as

$$ (3) \qquad \mathcal{L} = \frac{1}{2} \dot{\mathbf{q}}^T M \dot{\mathbf{q}} - \frac{1}{2} \mathbf{q}^T K \mathbf{q}, $$

where $M$ and $K$ are symmetric, positive-definite $n$-dimensional matrices with constant coefficients; $\mathbf{q} \in \mathbb{R}^n$.

We can discretize (2) by splitting the simulation interval $[0, T]$ into $N$ elements using a time step $h = T/N$. An approximation $q_j$ of $q(t_j)$ is calculated at each instance $t_j = jh$. The following action $S_d$ represents the discrete version of (2):

$$ (4) \qquad S_d = \sum_{j=1}^{N-1} L_d(\mathbf{q}_j, \mathbf{q}_{j+1}), $$

where $L_d(\mathbf{q}_\ell, \mathbf{q}_r)$ is the discrete form of the Lagrangian (3). Subscripts $\ell$ and $r$ stand for “left” and “right” values, respectively. Let us consider a centered finite difference approximation

$$ \frac{d\mathbf{q}}{dt} \approx \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h}, $$

and a midpoint quadrature

$$ \int_{0}^{h} \varphi(\mathbf{q}(t)) dt \approx h \varphi \left( \frac{\mathbf{q}_{\ell} + \mathbf{q}_{r}}{2} \right). $$
Numerical schemes for dynamical systems issued from the principle of least action are typically referred to as variational [17–20]. The general idea resides in performing a discretization at the least action principle level. As a result, the evolution equations deriving from this discretized principle characterize the system evolution but are also a numerical scheme. It is well known that such numerical methods are endowed with interesting characteristics; one characteristic is the property of being symplectic [18–21]. One remarkable example of such methods is Newmark’s integrator [17, 19], which is very popular for solving problems in dynamics of structures [22, 23], and has recently been geometrized to solve motion equations of sliding rods [24] and soft robots [25].

A symplectic scheme based on Simpson’s rule has been proposed by the authors in [26], for the linear and scalar case of the harmonic oscillator. The scheme uses a quadratic finite elements interpolation. The method was adapted to the monodimensional non-linear pendulum system in [27]. In this work, Simpson’s symplectic scheme is further studied as an alternative to Newmark’s method. It is generalized to the case of multiple degrees of freedom systems characterized by quadratic Lagrangians. Obtained results confirm the convergence rate previously observed in [26]. The new stability condition on the step size is revealed to be similar to the one previously obtained in [26]. A simplecticity analysis that applies to the multi-degrees-of-freedom case, along with the expression of a related conserved quadratic form is provided in this contribution.

We begin our study by detailing Newmark’s classical scheme, deriving it from variational principles in section 2. Then, Simpson’s alternative scheme is detailed and derived from variational principles in section 3. Section 4 analyzes the symplectic property of Simpson’s scheme. A proof that applies to both Newmark’s and Simpson’s schemes is provided. To compare both methods in a case study, a two-degrees of freedom system is presented. Therefore, the exact solution to the linearized double pendulum is provided in section 5. This exact solution serves for benchmarking purposes in our comparisons. Section 6 presents and comments on the obtained numerical results. Simpson’s scheme convergence is revealed to be of the fourth order. The manuscript ends with a brief discussion and concluding remarks in section 7.

## 2. NEWMARK’S SCHEME

### 2.1. Discrete action

Let us derive the classical, symplectic variational integrator based on Newmark’s scheme [17, 19, 22, 23]. The continuous action is defined as

$$ (2) \qquad S_c = \int_0^T \mathcal{L} \left( \frac{d\mathbf{q}(t)}{dt}, \mathbf{q}(t) \right) dt $$

where $\mathcal{L}$ is the system Lagrangian. We focus on dynamical systems for which the Lagrangian can be expressed quadratically as

$$ (3) \qquad \mathcal{L} = \frac{1}{2} \dot{\mathbf{q}}^T M \dot{\mathbf{q}} - \frac{1}{2} \mathbf{q}^T K \mathbf{q}, $$

where $M$ and $K$ are symmetric, positive-definite $n$-dimensional matrices with constant coefficients; $\mathbf{q} \in \mathbb{R}^n$.

We can discretize (2) by splitting the simulation interval $[0, T]$ into $N$ elements using a time step $h = T/N$. An approximation $q_j$ of $q(t_j)$ is calculated at each instance $t_j = jh$. The following action $S_d$ represents the discrete version of (2):

$$ (4) \qquad S_d = \sum_{j=1}^{N-1} L_d(\mathbf{q}_j, \mathbf{q}_{j+1}), $$

where $L_d(\mathbf{q}_\ell, \mathbf{q}_r)$ is the discrete form of the Lagrangian (3). Subscripts $\ell$ and $r$ stand for “left” and “right” values, respectively. Let us consider a centered finite difference approximation

$$ \frac{d\mathbf{q}}{dt} \approx \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h}, $$

and a midpoint quadrature

$$ \int_{0}^{h} \varphi(\mathbf{q}(t)) dt \approx h \varphi \left( \frac{\mathbf{q}_{\ell} + \mathbf{q}_{r}}{2} \right). $$
The discrete Lagrangian becomes

$$L_d(\mathbf{q}_\ell, \mathbf{q}_r) = \frac{h}{2} \left[ \left( \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h} \right)^T M \left( \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h} \right) \right] - \frac{h}{2} \left[ \left( \frac{\mathbf{q}_\ell + \mathbf{q}_r}{2} \right)^T K \left( \frac{\mathbf{q}_\ell + \mathbf{q}_r}{2} \right) \right].$$

## 2.2. Discrete Euler-Lagrange equations

The discrete action (4) being a sum, only two terms contain the variables $\mathbf{q}_j$:

$$S_d = \dots + L_d(\mathbf{q}_{j-1}, \mathbf{q}_j) + L_d(\mathbf{q}_j, \mathbf{q}_{j+1}) + \dots$$

So when the discrete action is stationary ($\delta S_d = 0$ for arbitrary variations $\delta \mathbf{q}_j$ of the states $\mathbf{q}_j$), only two terms remain. Necessarily,

$$(5) \qquad \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_{j-1}, \mathbf{q}_j) + \frac{\partial L_d}{\partial \mathbf{q}_\ell} (\mathbf{q}_j, \mathbf{q}_{j+1}) = 0.$$

The generalized momenta $\mathbf{p}_j \in \mathbb{R}^n$ are defined, on the right, as

$$(6) \qquad \mathbf{p}_j = \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_\ell, \mathbf{q}_r).$$

Therefore, the first term of equation (5) identifies as $\mathbf{p}_j$, so applying (6) on (5) leads to

$$(7) \qquad \mathbf{p}_j = -\frac{\partial L_d}{\partial \mathbf{q}_\ell} (\mathbf{q}_j, \mathbf{q}_{j+1}) = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) + \frac{h}{2} K \left( \frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2} \right).$$

Then, $\mathbf{p}_{j+1}$ is constructed following (6):

$$(8) \qquad \mathbf{p}_{j+1} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{2} K \left( \frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2} \right).$$

Using (7) and (8), it can be established that

$$(9) \qquad \frac{\mathbf{p}_{j+1} - \mathbf{p}_j}{h} = -K \left( \frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2} \right); \quad \frac{\mathbf{p}_j + \mathbf{p}_{j+1}}{2} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right).$$

Equations (9) are consistent with $\frac{d\mathbf{p}}{dt} = -K\mathbf{q}$ and $\mathbf{p} = M\frac{d\mathbf{q}}{dt}$, respectively.

## 2.3. Newmark's scheme

The system (9) can then be arranged in linear form as

$$(10) \qquad A_n \eta_{j+1} = B_n \eta_j,$$

where $\eta \in \mathbb{R}^{2n}$, $\eta = (\mathbf{p}, \mathbf{q})^T$ and

$$(11) \qquad A_n = \begin{pmatrix} \mathbb{I}_n & -X_n \\ \mathbb{I}_n & Y_n \end{pmatrix}; \quad B_n = \begin{pmatrix} -\mathbb{I}_n & -X_n \\ \mathbb{I}_n & -Y_n \end{pmatrix}; \quad X_n = \frac{2}{h}M; \quad Y_n = \frac{h}{2}K;$$

$\mathbb{I}_n$ being the $n$-dimensional identity matrix.

Newmark's symplectic scheme is obtained by matrix inversion of (10). We can establish that

$$(12) \qquad \eta_{j+1} = \Phi_n \eta_j, \quad \Phi_n = A_n^{-1} B,$$

where matrices $A$ and $B$ are defined in (11) above. It has been observed that this particular variant of Newmark's method is unconditionally stable and second-order convergent [19].

## 3. SIMPSON'S SCHEME

Newmark's scheme, presented in section 2, uses a midpoint quadrature for the numerical integration of a regular function. This quadrature is exact only for polynomials up to the first degree. A better precision is obtained with Simpson's quadrature (1), which is exact for polynomials up to the third degree. Notice how formula (1) introduces a midpoint. This midpoint will be regarded as an additional degree of freedom in our proposed integrator.

Let us now derive a symplectic scheme based on this integration rule. As with Newmark's scheme, the continuous action is defined by (2) and the Lagrangian has the structure of (3).
The discrete Lagrangian becomes

$$L_d(\mathbf{q}_\ell, \mathbf{q}_r) = \frac{h}{2} \left[ \left( \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h} \right)^T M \left( \frac{\mathbf{q}_r - \mathbf{q}_\ell}{h} \right) \right] - \frac{h}{2} \left[ \left( \frac{\mathbf{q}_\ell + \mathbf{q}_r}{2} \right)^T K \left( \frac{\mathbf{q}_\ell + \mathbf{q}_r}{2} \right) \right].$$

## 2.2. Discrete Euler-Lagrange equations

The discrete action (4) being a sum, only two terms contain the variables $\mathbf{q}_j$:

$$S_d = \dots + L_d(\mathbf{q}_{j-1}, \mathbf{q}_j) + L_d(\mathbf{q}_j, \mathbf{q}_{j+1}) + \dots$$

So when the discrete action is stationary ($\delta S_d = 0$ for arbitrary variations $\delta \mathbf{q}_j$ of the states $\mathbf{q}_j$), only two terms remain. Necessarily,

$$(5) \qquad \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_{j-1}, \mathbf{q}_j) + \frac{\partial L_d}{\partial \mathbf{q}_\ell} (\mathbf{q}_j, \mathbf{q}_{j+1}) = 0.$$

The generalized momenta $\mathbf{p}_j \in \mathbb{R}^n$ are defined, on the right, as

$$(6) \qquad \mathbf{p}_j = \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_\ell, \mathbf{q}_r).$$

Therefore, the first term of equation (5) identifies as $\mathbf{p}_j$, so applying (6) on (5) leads to

$$(7) \qquad \mathbf{p}_j = -\frac{\partial L_d}{\partial \mathbf{q}_\ell}(\mathbf{q}_j, \mathbf{q}_{j+1}) = M\left(\frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h}\right) + \frac{h}{2}K\left(\frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2}\right).$$

Then, $\mathbf{p}_{j+1}$ is constructed following (6):

$$(8) \qquad \mathbf{p}_{j+1} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{2} K \left( \frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2} \right).$$

Using (7) and (8), it can be established that

$$(9) \qquad \frac{\mathbf{p}_{j+1} - \mathbf{p}_j}{h} = -K \left( \frac{\mathbf{q}_j + \mathbf{q}_{j+1}}{2} \right); \quad \frac{\mathbf{p}_j + \mathbf{p}_{j+1}}{2} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right).$$

Equations (9) are consistent with $\frac{d\mathbf{p}}{dt} = -K\mathbf{q}$ and $\mathbf{p} = M\frac{d\mathbf{q}}{dt}$, respectively.

## 2.3. Newmark's scheme

The system (9) can then be arranged in linear form as

$$(10) \qquad A_n \eta_{j+1} = B_n \eta_j,$$

where $\eta \in \mathbb{R}^{2n}$, $\eta = (\mathbf{p}, \mathbf{q})^T$ and

$$(11) \qquad A_n = \begin{pmatrix} \mathbb{I}_n & -X_n \\ \mathbb{I}_n & Y_n \end{pmatrix}; \quad B_n = \begin{pmatrix} -\mathbb{I}_n & -X_n \\ \mathbb{I}_n & -Y_n \end{pmatrix}; \quad X_n = \frac{2}{h}M; \quad Y_n = \frac{h}{2}K;$$

$\mathbb{I}_n$ being the $n$-dimensional identity matrix.

Newmark's symplectic scheme is obtained by matrix inversion of (10). We can establish that

$$(12) \qquad \eta_{j+1} = \Phi_n \eta_j, \quad \Phi_n = A_n^{-1} B,$$

where matrices $A$ and $B$ are defined in (11) above. It has been observed that this particular variant of Newmark's method is unconditionally stable and second-order convergent [19].

## 3. SIMPSON'S SCHEME

Newmark's scheme, presented in section 2, uses a midpoint quadrature for the numerical integration of a regular function. This quadrature is exact only for polynomials up to the first degree. A better precision is obtained with Simpson's quadrature (1), which is exact for polynomials up to the third degree. Notice how formula (1) introduces a midpoint. This midpoint will be regarded as an additional degree of freedom in our proposed integrator.

Let us now derive a symplectic scheme based on this integration rule. As with Newmark's scheme, the continuous action is defined by (2) and the Lagrangian has the structure of (3).
### 3.1. Quadratic finite elements interpolation

An internal interpolation can be performed at each time step, for $t \in [0, h]$, using quadratic finite elements [16, 28]. We use the following compact basis functions for $0 \le \theta \le 1$:

$$ (13) \qquad \varphi_0(\theta) = (1-\theta)(1-2\theta), \quad \varphi_{1/2}(\theta) = 4\theta(1-\theta), \quad \varphi_1(\theta) = \theta(2\theta-1). $$

At $t = h\theta$, the states $\mathbf{q}(t) \in P_2$ are approximated with the above basis functions as

$$ (14) \qquad \mathbf{q}(t) = \mathbf{q}_{\ell}\varphi_0(\theta) + \mathbf{q}_m\varphi_{1/2}(\theta) + \mathbf{q}_r\varphi_1(\theta). $$

Note that $\mathbf{q}(0) = \mathbf{q}_{\ell}, \mathbf{q}(\frac{h}{2}) = \mathbf{q}_m$ and $\mathbf{q}(h) = \mathbf{q}_r$; here, subscript $m$ stands for "middle". This means that the finite elements (13) are well adapted to the internal degree of freedom at $h/2$. Then, by time differentiation,

$$ \begin{aligned} \frac{d\mathbf{q}}{dt} &= \frac{1}{h} \left( \mathbf{q}_{\ell} \frac{d\varphi_0}{d\theta} + \mathbf{q}_m \frac{d\varphi_{1/2}}{d\theta} + \mathbf{q}_r \frac{d\varphi_1}{d\theta} \right) \\ &= \frac{1}{h} (\mathbf{q}_{\ell}(4\theta - 3) + 4\mathbf{q}_m(1 - 2\theta) + \mathbf{q}_r(4\theta - 1)) \\ &= \mathbf{g}_{\ell}(1 - \theta) + \mathbf{g}_r\theta \end{aligned} $$

where derivatives $\mathbf{g}_{\ell}, \mathbf{g}_r \in \mathbb{R}^n$ are given by Gear's scheme [29]. Gear's scheme is used as the differentiation approximation for $\mathbf{q}(t) \in P_2$ as

$$ \begin{aligned} \mathbf{g}_{\ell} &= \frac{d\mathbf{q}}{dt}(0) = \frac{1}{h}(-3\mathbf{q}_{\ell} + 4\mathbf{q}_m - \mathbf{q}_r), \\ \mathbf{g}_m &= \frac{d\mathbf{q}}{dt}(\frac{h}{2}) = \frac{\mathbf{g}_{\ell} + \mathbf{g}_r}{2} = \frac{\mathbf{q}_r - \mathbf{q}_{\ell}}{h}, \\ \mathbf{g}_r &= \frac{d\mathbf{q}}{dt}(h) = \frac{1}{h}(\mathbf{q}_{\ell} - 4\mathbf{q}_m + 3\mathbf{q}_r), \end{aligned} \quad (15) $$

where $\mathbf{g}_m \in \mathbb{R}^n$. The above confirms that a first-order centered finite difference is recovered by $\mathbf{g}_m$, which is the derivative at the middle of the discretization interval.

The interpolation is used within an interval of length $h$ by splitting the range $[0, T]$ into $N$ pieces, giving a fixed step size of $h = T/N$. At each discrete time instance $t_j = jh$, we have

$$ \begin{aligned} \mathbf{q}_j &\approx \mathbf{q}(t_j), & \forall 0 \le j \le N; \\ \mathbf{q}_{j+1/2} &\approx \mathbf{q}(t_j + \frac{h}{2}), & \forall 0 \le j \le N-1. \end{aligned} $$

Taking (14), $\mathbf{q}(t)$ is a quadratic polynomial vector function within the interval $[t_j, t_{j+1}]$ with

$$ t = t_j + \theta h, \quad \mathbf{q}_{\ell} = \mathbf{q}_j, \quad \mathbf{q}_m = \mathbf{q}_{j+1/2}, \quad \mathbf{q}_r = \mathbf{q}_{j+1}. $$

### 3.2. Discrete Lagrangian

Let us recall that the continuous action is defined by (2) and that the Lagrangian is defined by (3). In the present case, the discrete action sum $\Sigma_d$ for a motion $t \mapsto \mathbf{q}(t)$ between the initial time and a given final time $T > 0$, is discretized with $N$ regular intervals as

$$ (16) \qquad \Sigma_d = \sum_{j=1}^{N-1} L_h(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}), $$

where $L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r)$ is the discrete form

$$ L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r) \approx \int_0^h \mathcal{L} \, dt, $$

of the Lagrangian (3). Using Simpson's rule (1), the polynomial approximations (14) of the states, and derivatives (15), the discrete Lagrangian of a linear system is expressed as

$$ L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r) = \frac{h}{2} \left[ \frac{1}{6} \mathbf{g}_{\ell}^T M \mathbf{g}_{\ell} + \frac{2}{3} \mathbf{g}_m^T M \mathbf{g}_m + \frac{1}{6} \mathbf{g}_r^T M \mathbf{g}_r \right] \\ - \frac{h}{2} \left[ \frac{1}{6} \mathbf{q}_{\ell}^T K \mathbf{q}_{\ell} + \frac{2}{3} \mathbf{q}_m^T K \mathbf{q}_m + \frac{1}{6} \mathbf{q}_r^T K \mathbf{q}_r \right]. $$
### 3.1. Quadratic finite elements interpolation

An internal interpolation can be performed at each time step, for $t \in [0, h]$, using quadratic finite elements [16, 28]. We use the following compact basis functions for $0 \le \theta \le 1$:

$$ (13) \qquad \varphi_0(\theta) = (1-\theta)(1-2\theta), \quad \varphi_{1/2}(\theta) = 4\theta(1-\theta), \quad \varphi_1(\theta) = \theta(2\theta-1). $$

At $t = h\theta$, the states $\mathbf{q}(t) \in P_2$ are approximated with the above basis functions as

$$ (14) \qquad \mathbf{q}(t) = \mathbf{q}_{\ell}\varphi_0(\theta) + \mathbf{q}_m\varphi_{1/2}(\theta) + \mathbf{q}_r\varphi_1(\theta). $$

Note that $\mathbf{q}(0) = \mathbf{q}_{\ell}, \mathbf{q}(\frac{h}{2}) = \mathbf{q}_m$ and $\mathbf{q}(h) = \mathbf{q}_r$; here, subscript $m$ stands for "middle". This means that the finite elements (13) are well adapted to the internal degree of freedom at $h/2$. Then, by time differentiation,

$$ \begin{aligned} \frac{d\mathbf{q}}{dt} &= \frac{1}{h} \left( \mathbf{q}_{\ell} \frac{d\varphi_0}{d\theta} + \mathbf{q}_m \frac{d\varphi_{1/2}}{d\theta} + \mathbf{q}_r \frac{d\varphi_1}{d\theta} \right) \\ &= \frac{1}{h} (\mathbf{q}_{\ell}(4\theta - 3) + 4\mathbf{q}_m(1 - 2\theta) + \mathbf{q}_r(4\theta - 1)) \\ &= \mathbf{g}_{\ell}(1 - \theta) + \mathbf{g}_r\theta \end{aligned} $$

where derivatives $\mathbf{g}_{\ell}, \mathbf{g}_r \in \mathbb{R}^n$ are given by Gear's scheme [29]. Gear's scheme is used as the differentiation approximation for $\mathbf{q}(t) \in P_2$ as

$$ \begin{aligned} \mathbf{g}_{\ell} &= \frac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}(0) = \frac{1}{h}(-3\mathbf{q}_{\ell} + 4\mathbf{q}_m - \mathbf{q}_r), \\ \mathbf{g}_m &= \frac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}\left(\frac{h}{2}\right) = \frac{\mathbf{g}_{\ell} + \mathbf{g}_r}{2} = \frac{\mathbf{q}_r - \mathbf{q}_{\ell}}{h}, \\ \mathbf{g}_r &= \frac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}(h) = \frac{1}{h}(\mathbf{q}_{\ell} - 4\mathbf{q}_m + 3\mathbf{q}_r), \end{aligned} \quad (15) $$

where $\mathbf{g}_m \in \mathbb{R}^n$. The above confirms that a first-order centered finite difference is recovered by $\mathbf{g}_m$, which is the derivative at the middle of the discretization interval.

The interpolation is used within an interval of length $h$ by splitting the range $[0, T]$ into $N$ pieces, giving a fixed step size of $h = T/N$. At each discrete time instance $t_j = jh$, we have

$$ \begin{aligned} \mathbf{q}_j &\approx \mathbf{q}(t_j), & \forall 0 \le j \le N; \\ \mathbf{q}_{j+1/2} &\approx \mathbf{q}\left(t_j + \frac{h}{2}\right), & \forall 0 \le j \le N-1. \end{aligned} $$

Taking (14), $\mathbf{q}(t)$ is a quadratic polynomial vector function within the interval $[t_j, t_{j+1}]$ with

$$ t = t_j + \theta h, \quad \mathbf{q}_{\ell} = \mathbf{q}_j, \quad \mathbf{q}_m = \mathbf{q}_{j+1/2}, \quad \mathbf{q}_r = \mathbf{q}_{j+1}. $$

### 3.2. Discrete Lagrangian

Let us recall that the continuous action is defined by (2) and that the Lagrangian is defined by (3). In the present case, the discrete action sum $\Sigma_d$ for a motion $t \mapsto \mathbf{q}(t)$ between the initial time and a given final time $T > 0$, is discretized with $N$ regular intervals as

$$ (16) \qquad \Sigma_d = \sum_{j=1}^{N-1} L_h(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}), $$

where $L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r)$ is the discrete form

$$ L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r) \approx \int_0^h \mathcal{L} \, dt, $$

of the Lagrangian (3). Using Simpson's rule (1), the polynomial approximations (14) of the states, and derivatives (15), the discrete Lagrangian of a linear system is expressed as

$$ L_h(\mathbf{q}_{\ell}, \mathbf{q}_m, \mathbf{q}_r) = \frac{h}{2} \left[ \frac{1}{6} \mathbf{g}_{\ell}^T M \mathbf{g}_{\ell} + \frac{2}{3} \mathbf{g}_m^T M \mathbf{g}_m + \frac{1}{6} \mathbf{g}_r^T M \mathbf{g}_r \right] \\ - \frac{h}{2} \left[ \frac{1}{6} \mathbf{q}_{\ell}^T K \mathbf{q}_{\ell} + \frac{2}{3} \mathbf{q}_m^T K \mathbf{q}_m + \frac{1}{6} \mathbf{q}_r^T K \mathbf{q}_r \right]. $$
### 3.3. Discrete Euler-Lagrange equations

Recall that Simpson's rule introduces an internal degree of freedom in the middle of the interpolation interval. The discrete action (16) is a sum where only two terms contain the variables $\mathbf{q}_j$ and $\mathbf{q}_{j+1/2}$:

$$ \Sigma_d = \dots + L_h(\mathbf{q}_{j-1}, \mathbf{q}_{j-1/2}, \mathbf{q}_j) + L_h(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) + \dots $$

Maupertuis's stationary action principle [10] implies that $\delta\Sigma_d = 0$ for an arbitrary variation of the internal degree of freedom $\delta\mathbf{q}_{j+1/2} \in [t_j, t_{j+1}]$. Considering Gear's scheme (15), we have

$$ \frac{\partial g_{\ell}^{i}}{\partial q_{m}^{k}} = \frac{4}{h'}, \quad \frac{\partial g_{m}^{i}}{\partial q_{m}^{k}} = 0, \quad \frac{\partial g_{r}^{i}}{\partial q_{m}^{k}} = -\frac{4}{h'}, \quad \forall i = k, 0 \text{ otherwise,} $$

where $g^i$ is the i-th component of $\mathbf{g}$ and $q^k$ is the k-th component of $\mathbf{q}$.

When $\delta\Sigma_d = 0$, $\frac{\partial L_h}{\partial \mathbf{q}_m}(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) = 0$ by necessity. This conforms to the discrete Euler-Lagrange equations at the middle of the interpolation interval:

$$ (17) \qquad \frac{4}{3h} M \mathbf{g}_j - \frac{4}{3h} M \mathbf{g}_{j+1} - \frac{4}{3} K \mathbf{q}_{j+1/2} = 0. $$

However, $\mathbf{g}_j - \mathbf{g}_{j+1} = \frac{1}{h}(-4\mathbf{q}_j + 8\mathbf{q}_{j+1/2} - 4\mathbf{q}_{j+1})$, so equation (17) becomes

$$ (18) \qquad M \left( 4 \frac{\mathbf{q}_j - 2\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1}}{h^2} \right) + K \mathbf{q}_{j+1/2} = 0. $$

This last equation is consistent with $M \frac{d^2\mathbf{q}}{dt^2} + K\mathbf{q} = 0$.

Additionally, for an arbitrary variation $\delta\mathbf{q}_j$, the Euler-Lagrange equations are given by the necessary condition that

$$ (19) \qquad \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_{j-1}, \mathbf{q}_{j-1/2}, \mathbf{q}_j) + \frac{\partial L_d}{\partial \mathbf{q}_\ell} (\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) = 0. $$

The generalized momenta $\mathbf{p}_j$ are defined, on the right, as

$$ (20) \qquad \mathbf{p}_j = \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_\ell, \mathbf{q}_m, \mathbf{q}_r). $$

Therefore, the first term of equation (19) identifies as $\mathbf{p}_j$, and it can established that

$$ \begin{align} \mathbf{p}_j &= -\frac{\partial L_d}{\partial \mathbf{q}_\ell} (\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) \nonumber \\ (21) \qquad &= -\frac{h}{2} \left[ -\frac{3}{3h} M \mathbf{g}_j - \frac{4}{3h} M \mathbf{g}_{j+1/2} + \frac{1}{3h} M \mathbf{g}_{j+1} \right] + \frac{h}{6} K \mathbf{q}_j \nonumber \\ &= -\frac{1}{6h} M \left( 14 \mathbf{q}_j - 16 \mathbf{q}_{j+1/2} + 2 \mathbf{q}_{j+1} \right) + \frac{h}{6} K \mathbf{q}_j, \end{align} $$

because $-3\mathbf{q}_j - 4\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1} = \frac{1}{h}(14\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 2\mathbf{q}_{j+1})$. Equation (18) is then multiplied by $h/3$ and the result is added to equation (21). This eliminates $\mathbf{q}_{j+1/2}$ from the first term of the right hand side:

$$ (22) \qquad \mathbf{p}_j = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{6} K \left( -2\mathbf{q}_{j+1/2} - \mathbf{q}_j \right). $$

Then, $\mathbf{p}_{j+1}$ is calculated according to equation (20):

$$ \begin{align} \mathbf{p}_{j+1} &= \frac{\partial L_d}{\partial \mathbf{q}_r} (\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) \tag{23} \\ &= \frac{h}{2} \left[ -\frac{1}{3h} M \mathbf{g}_j + \frac{4}{3h} M \mathbf{g}_{j+1/2} + \frac{3}{3h} M \mathbf{g}_{j+1} \right] - \frac{h}{6} K \mathbf{q}_{j+1} \nonumber \\ &= \frac{1}{6h} M \left( 2\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 14\mathbf{q}_{j+1} \right) - \frac{h}{6} K \mathbf{q}_{j+1} \nonumber \end{align} $$
### 3.3. Discrete Euler-Lagrange equations

Recall that Simpson's rule introduces an internal degree of freedom in the middle of the interpolation interval. The discrete action (16) is a sum where only two terms contain the variables $\mathbf{q}_j$ and $\mathbf{q}_{j+1/2}$:

$$ \Sigma_d = \dots + L_h(\mathbf{q}_{j-1}, \mathbf{q}_{j-1/2}, \mathbf{q}_j) + L_h(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) + \dots $$

Maupertuis's stationary action principle [10] implies that $\delta\Sigma_d = 0$ for an arbitrary variation of the internal degree of freedom $\delta\mathbf{q}_{j+1/2} \in [t_j, t_{j+1}]$. Considering Gear's scheme (15), we have

$$ \frac{\partial g_{\ell}^{i}}{\partial q_{m}^{k}} = \frac{4}{h'}, \quad \frac{\partial g_{m}^{i}}{\partial q_{m}^{k}} = 0, \quad \frac{\partial g_{r}^{i}}{\partial q_{m}^{k}} = -\frac{4}{h'}, \quad \forall i = k, 0 \text{ otherwise,} $$

where $g^i$ is the i-th component of $\mathbf{g}$ and $q^k$ is the k-th component of $\mathbf{q}$.

When $\delta\Sigma_d = 0$, $\frac{\partial L_h}{\partial \mathbf{q}_m}(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) = 0$ by necessity. This conforms to the discrete Euler-Lagrange equations at the middle of the interpolation interval:

$$ (17) \qquad \frac{4}{3h} M \mathbf{g}_j - \frac{4}{3h} M \mathbf{g}_{j+1} - \frac{4}{3} K \mathbf{q}_{j+1/2} = 0. $$

However, $\mathbf{g}_j - \mathbf{g}_{j+1} = \frac{1}{h}(-4\mathbf{q}_j + 8\mathbf{q}_{j+1/2} - 4\mathbf{q}_{j+1})$, so equation (17) becomes

$$ (18) \qquad M \left( 4 \frac{\mathbf{q}_j - 2\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1}}{h^2} \right) + K \mathbf{q}_{j+1/2} = 0. $$

This last equation is consistent with $M \frac{d^2\mathbf{q}}{dt^2} + K\mathbf{q} = 0$.

Additionally, for an arbitrary variation $\delta\mathbf{q}_j$, the Euler-Lagrange equations are given by the necessary condition that

$$ (19) \qquad \frac{\partial L_d}{\partial \mathbf{q}_r}(\mathbf{q}_{j-1}, \mathbf{q}_{j-1/2}, \mathbf{q}_j) + \frac{\partial L_d}{\partial \mathbf{q}_\ell}(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) = 0. $$

The generalized momenta $\mathbf{p}_j$ are defined, on the right, as

$$ (20) \qquad \mathbf{p}_j = \frac{\partial L_d}{\partial \mathbf{q}_r}(\mathbf{q}_\ell, \mathbf{q}_m, \mathbf{q}_r). $$

Therefore, the first term of equation (19) identifies as $\mathbf{p}_j$, and it can established that

$$ \begin{align} \mathbf{p}_j &= -\frac{\partial L_d}{\partial \mathbf{q}_\ell}(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) \nonumber \\ (21) \qquad &= -\frac{h}{2} \left[ -\frac{3}{3h} M \mathbf{g}_j - \frac{4}{3h} M \mathbf{g}_{j+1/2} + \frac{1}{3h} M \mathbf{g}_{j+1} \right] + \frac{h}{6} K \mathbf{q}_j \nonumber \\ &= -\frac{1}{6h} M \left( 14 \mathbf{q}_j - 16 \mathbf{q}_{j+1/2} + 2 \mathbf{q}_{j+1} \right) + \frac{h}{6} K \mathbf{q}_j, \end{align} $$

because $-3\mathbf{q}_j - 4\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1} = \frac{1}{h}(14\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 2\mathbf{q}_{j+1})$. Equation (18) is then multiplied by $h/3$ and the result is added to equation (21). This eliminates $\mathbf{q}_{j+1/2}$ from the first term of the right hand side:

$$ (22) \qquad \mathbf{p}_j = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{6} K \left( -2\mathbf{q}_{j+1/2} - \mathbf{q}_j \right). $$

Then, $\mathbf{p}_{j+1}$ is calculated according to equation (20):

$$ \begin{align} \mathbf{p}_{j+1} &= \frac{\partial L_d}{\partial \mathbf{q}_r}(\mathbf{q}_j, \mathbf{q}_{j+1/2}, \mathbf{q}_{j+1}) \tag{23} \\ &= \frac{h}{2} \left[ -\frac{1}{3h} M \mathbf{g}_j + \frac{4}{3h} M \mathbf{g}_{j+1/2} + \frac{3}{3h} M \mathbf{g}_{j+1} \right] - \frac{h}{6} K \mathbf{q}_{j+1} \nonumber \\ &= \frac{1}{6h} M \left( 2\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 14\mathbf{q}_{j+1} \right) - \frac{h}{6} K \mathbf{q}_{j+1} \nonumber \end{align} $$
because $-\mathbf{g}_j + 4\mathbf{g}_{j+1/2} + 3\mathbf{g}_{j+1} = \frac{1}{h}(2\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 14\mathbf{q}_{j+1})$. Equation (18) is then multiplied by $-h^{3/3}$ and the result is added to equation (23). This eliminates $\mathbf{q}_{j+1/2}$ from the first term of the right hand side:

$$ (24) \qquad \mathbf{p}_{j+1} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{6} K \left( 2\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1} \right). $$

Using (22) and (24) we can establish that

$$ (25) \qquad \left\{ \begin{aligned} \frac{\mathbf{p}_{j+1} - \mathbf{p}_j}{h} &= -K \left( \frac{1}{6}\mathbf{q}_j + \frac{2}{3}\mathbf{q}_{j+1/2} + \frac{1}{6}\mathbf{q}_{j+1} \right) \\ \frac{\mathbf{p}_j + \mathbf{p}_{j+1}}{2} &= \left( M - \frac{h^2}{12}K \right) \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right). \end{aligned} \right. $$

Equations (25) are consistent with $\frac{d\mathbf{p}}{dt} = -K\mathbf{q}$, and $\mathbf{p} = M\frac{d\mathbf{q}}{dt}$, respectively. Note that the term $\frac{h^2}{12}K$ in the second equation above vanishes as $h \to 0$.

### 3.4. First variant of Simpson's scheme

The system composed of equations (18) and (25) can be rearranged as

$$ (26) \qquad \left\{ \begin{aligned} & L \mathbf{q}_{j+1/2} - \frac{1}{2} \mathbf{q}_{j+1} = \frac{1}{2} \mathbf{q}_j \\ & \mathbf{p}_{j+1} - \left( \frac{2}{h} M - \frac{h}{6} K \right) \mathbf{q}_{j+1} = - \mathbf{p}_j - \left( \frac{2}{h} M - \frac{h}{6} K \right) \mathbf{q}_j \\ & \frac{2h}{3} K \mathbf{q}_{j+1/2} + \mathbf{p}_{j+1} + \frac{h}{6} K \mathbf{q}_{j+1} = \mathbf{p}_j - \frac{h}{6} K \mathbf{q}_j \end{aligned} \right. $$

where

$$ L = \left( \mathbb{I}_n - \frac{h^2}{8} M^{-1} K \right). $$

The system (26) can then be arranged in linear form:

$$ (27) \qquad A_{\sigma} \begin{pmatrix} \mathbf{q}_{j+1/2} \\ \eta_{j+1} \end{pmatrix} = B_{\sigma} \eta_{j+1}, $$

where $\eta \in \mathbb{R}^{2n}, \eta = (\mathbf{p}, \mathbf{q})^T$;

$$ (28) \qquad A = \begin{pmatrix} L & 0 & -\frac{1}{2}\mathbb{I}_n \\ 0 & \mathbb{I}_n & -\left(\frac{2}{h}M - \frac{h}{6}K\right) \\ \frac{2h}{3}K & \mathbb{I}_n & \frac{h}{6}K \end{pmatrix}; \quad B = \begin{pmatrix} 0 & \frac{1}{2}\mathbb{I}_n \\ -\mathbb{I}_n & -\left(\frac{2}{h}M - \frac{h}{6}K\right) \\ \mathbb{I}_n & -\frac{h}{6}K \end{pmatrix}. $$

The first variant of Simpson's scheme is obtained by matrix inversion of (27). We can establish that

$$ (29) \qquad \begin{pmatrix} \mathbf{q}_{j+1/2} \\ \eta_{j+1} \end{pmatrix} = A_{\sigma}^{-1} B_{\sigma} \eta_{j+1}, $$

where matrices $A_{\sigma}$ and $B_{\sigma}$ are defined in (28) above.

### 3.5. Second variant of Simpson's scheme

Simpson's scheme internal degree of freedom can be eliminated using the first equation of (26):

$$ \mathbf{q}_{j+1/2} = \frac{1}{2} L^{-1} (\mathbf{q}_{j+1} + \mathbf{q}_j). $$

This equation approximates the middle point when $h \to 0$ because then $L \to \mathbb{I}_n$. Substituting this value into the third equation of (26) leads to

$$ \mathbf{p}_{j+1} + \frac{h}{3} \left( KL^{-1} + \frac{1}{2}K \right) \mathbf{q}_{j+1} = \mathbf{p}_j - \frac{h}{3} \left( KL^{-1} + \frac{1}{2}K \right) \mathbf{q}_j, $$

and the second equation of (26) remains unchanged. Therefore, the internal degree of freedom is successfully eliminated so that now

$$ (30) \qquad A_s \eta_{j+1} = B_s \eta_j $$
because $-\mathbf{g}_j + 4\mathbf{g}_{j+1/2} + 3\mathbf{g}_{j+1} = \frac{1}{h}(2\mathbf{q}_j - 16\mathbf{q}_{j+1/2} + 14\mathbf{q}_{j+1})$. Equation (18) is then multiplied by $-h^{3/3}$ and the result is added to equation (23). This eliminates $\mathbf{q}_{j+1/2}$ from the first term of the right hand side:

$$ (24) \qquad \mathbf{p}_{j+1} = M \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right) - \frac{h}{6} K \left( 2\mathbf{q}_{j+1/2} + \mathbf{q}_{j+1} \right). $$

Using (22) and (24) we can establish that

$$ (25) \qquad \left\{ \begin{aligned} \frac{\mathbf{p}_{j+1} - \mathbf{p}_j}{h} &= -K \left( \frac{1}{6}\mathbf{q}_j + \frac{2}{3}\mathbf{q}_{j+1/2} + \frac{1}{6}\mathbf{q}_{j+1} \right) \\ \frac{\mathbf{p}_j + \mathbf{p}_{j+1}}{2} &= \left( M - \frac{h^2}{12}K \right) \left( \frac{\mathbf{q}_{j+1} - \mathbf{q}_j}{h} \right). \end{aligned} \right. $$

Equations (25) are consistent with $\frac{d\mathbf{p}}{dt} = -K\mathbf{q}$, and $\mathbf{p} = M\frac{d\mathbf{q}}{dt}$, respectively. Note that the term $\frac{h^2}{12}K$ in the second equation above vanishes as $h \to 0$.

### 3.4. First variant of Simpson's scheme

The system composed of equations (18) and (25) can be rearranged as

$$ (26) \qquad \left\{ \begin{aligned} & L \mathbf{q}_{j+1/2} - \frac{1}{2} \mathbf{q}_{j+1} = \frac{1}{2} \mathbf{q}_j \\ & \mathbf{p}_{j+1} - \left( \frac{2}{h} M - \frac{h}{6} K \right) \mathbf{q}_{j+1} = - \mathbf{p}_j - \left( \frac{2}{h} M - \frac{h}{6} K \right) \mathbf{q}_j \\ & \frac{2h}{3} K \mathbf{q}_{j+1/2} + \mathbf{p}_{j+1} + \frac{h}{6} K \mathbf{q}_{j+1} = \mathbf{p}_j - \frac{h}{6} K \mathbf{q}_j \end{aligned} \right. $$

where

$$ L = \left( \mathbb{I}_n - \frac{h^2}{8} M^{-1} K \right). $$

The system (26) can then be arranged in linear form:

$$ (27) \qquad A_{\sigma} \begin{pmatrix} \mathbf{q}_{j+1/2} \\ \eta_{j+1} \end{pmatrix} = B_{\sigma} \eta_{j+1}, $$

where $\eta \in \mathbb{R}^{2n}$, $\eta = (\mathbf{p}, \mathbf{q})^T$;

$$ (28) \qquad A = \begin{pmatrix} L & 0 & -\frac{1}{2}\mathbb{I}_n \\ 0 & \mathbb{I}_n & -\left(\frac{2}{h}M - \frac{h}{6}K\right) \\ \frac{2h}{3}K & \mathbb{I}_n & \frac{h}{6}K \end{pmatrix}; \quad B = \begin{pmatrix} 0 & \frac{1}{2}\mathbb{I}_n \\ -\mathbb{I}_n & -\left(\frac{2}{h}M - \frac{h}{6}K\right) \\ \mathbb{I}_n & -\frac{h}{6}K \end{pmatrix}. $$

The first variant of Simpson's scheme is obtained by matrix inversion of (27). We can establish that

$$ (29) \qquad \begin{pmatrix} \mathbf{q}_{j+1/2} \\ \eta_{j+1} \end{pmatrix} = A_{\sigma}^{-1} B_{\sigma} \eta_{j+1}, $$

where matrices $A_{\sigma}$ and $B_{\sigma}$ are defined in (28) above.

### 3.5. Second variant of Simpson's scheme

Simpson's scheme internal degree of freedom can be eliminated using the first equation of (26):

$$ \mathbf{q}_{j+1/2} = \frac{1}{2} L^{-1} (\mathbf{q}_{j+1} + \mathbf{q}_j). $$

This equation approximates the middle point when $h \to 0$ because then $L \to \mathbb{I}_n$. Substituting this value into the third equation of (26) leads to

$$ \mathbf{p}_{j+1} + \frac{h}{3} \left( KL^{-1} + \frac{1}{2}K \right) \mathbf{q}_{j+1} = \mathbf{p}_j - \frac{h}{3} \left( KL^{-1} + \frac{1}{2}K \right) \mathbf{q}_j, $$

and the second equation of (26) remains unchanged. Therefore, the internal degree of freedom is successfully eliminated so that now

$$ (30) \qquad A_s \eta_{j+1} = B_s \eta_j $$
where $\eta \in \mathbb{R}^{2n}; \eta = (\mathbf{p}, \mathbf{q})^T$;

$$ (31) \qquad A_s = \begin{pmatrix} \mathbb{I}_n & -X_s \\ \mathbb{I}_n & Y_s \end{pmatrix}; \quad B_s = \begin{pmatrix} -\mathbb{I}_n & -X_s \\ \mathbb{I}_n & -Y_s \end{pmatrix}; \quad X_s = \frac{2}{h}M - \frac{h}{6}K; \quad Y_s = \frac{h}{3}\left(KL^{-1} + \frac{1}{2}K\right). $$

The second variant of Simpson's symplectic scheme is obtained by matrix inversion of (30). We can establish that

$$ (32) \qquad \eta_{j+1} = \Phi_s \eta_j, \quad \Phi_s = A_s^{-1} B_s, $$

where matrices $A_s$ and $B_s$ are defined in (31) above. Note that schemes (29) and (32) are equivalent. However, this second variant eliminates the internal degree of freedom in the middle of the interval.

The symplecticity of Simpson's scheme (32) has not yet been proven. However, one can appreciate the similarity with Newmark's scheme by comparing (11) and (31). The symplectic property of both schemes is analyzed in the next section.

## 4. SYMPLECTICITY OF NEWMARK'S AND SIMPSON'S SCHEMES

The symplectic property of both Newmark's scheme (12) and Simpson's scheme (32) is now analyzed.

### 4.1. Symplectic property

A symplecticity proof is obtained by verifying that

$$ (33) \qquad \Phi^T \mathcal{J} \Phi = \mathcal{J}; \quad \mathcal{J} = \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & 0 \end{pmatrix}. $$

$\Phi$ corresponds to the scheme transformation matrix and characterizes a discrete time evolution of the system. $\mathcal{J}$ is sometimes referred to as the canonical matrix for Hamiltonian systems [30] and has the property that $\mathcal{J}^{-1} = \mathcal{J}^T = -\mathcal{J}$. When equation (33) holds, it means that $\Phi$ is an area-preserving transformation and that the scheme (12) is symplectic (see e. g. [18–20,31] for more details on this demonstration).

**Proposition 1.** *An implicit scheme of the type*

$$ \eta_{j+1} = A^{-1}B\eta_j; \quad \eta = (\mathbf{p}, \mathbf{q})^T, $$

is symplectic if

$$ (34) \qquad A = \begin{pmatrix} \mathbb{I}_n & -X \\ \mathbb{I}_n & Y \end{pmatrix}, \quad B = \begin{pmatrix} -\mathbb{I}_n & -X \\ \mathbb{I}_n & -Y \end{pmatrix}, $$

are square, partitioned, invertible matrices; and blocks X and Y are symmetric and positive-definite.

*Proof of Proposition 1.* Let us first explicit the transformation $A^{-1}$. Since A is square and partitioned, its inversion is done using auxiliary variables $\alpha$ and $\beta$. Let us establish that

$$ (35) \qquad A \begin{pmatrix} \mathbf{p} \\ \mathbf{q} \end{pmatrix} = \begin{pmatrix} \mathbf{p} - X\mathbf{q} \\ \mathbf{p} + Y\mathbf{q} \end{pmatrix} = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} $$

Subtracting both equations above gives

$$ (36) \qquad \mathbf{q} = Z^{-1}(\beta - \alpha); \quad Z = X + Y. $$

Since Z is the sum of two symmetric, positive-definite matrices, it is invertible. Equation (36) is then substituted into the first equation of (35):

$$ (37) \qquad \mathbf{p} = (\mathbb{I}_n - XZ^{-1})\alpha + XZ^{-1}\beta. $$

Matrix A (34) is inverted in (36) and (37) as:

$$ A^{-1} = \begin{pmatrix} \mathbb{I}_n - XZ^{-1} & XZ^{-1} \\ -Z^{-1} & Z^{-1} \end{pmatrix}. $$
where $\eta \in \mathbb{R}^{2n}; \eta = (\mathbf{p}, \mathbf{q})^T$;

$$ (31) \qquad A_s = \begin{pmatrix} \mathbb{I}_n & -X_s \\ \mathbb{I}_n & Y_s \end{pmatrix}; \quad B_s = \begin{pmatrix} -\mathbb{I}_n & -X_s \\ \mathbb{I}_n & -Y_s \end{pmatrix}; \quad X_s = \frac{2}{h}M - \frac{h}{6}K; \quad Y_s = \frac{h}{3}\left(KL^{-1} + \frac{1}{2}K\right). $$

The second variant of Simpson's symplectic scheme is obtained by matrix inversion of (30). We can establish that

$$ (32) \qquad \eta_{j+1} = \Phi_s \eta_j, \quad \Phi_s = A_s^{-1} B_s, $$

where matrices $A_s$ and $B_s$ are defined in (31) above. Note that schemes (29) and (32) are equivalent. However, this second variant eliminates the internal degree of freedom in the middle of the interval.

The symplecticity of Simpson's scheme (32) has not yet been proven. However, one can appreciate the similarity with Newmark's scheme by comparing (11) and (31). The symplectic property of both schemes is analyzed in the next section.

## 4. SYMPLECTICITY OF NEWMARK'S AND SIMPSON'S SCHEMES

The symplectic property of both Newmark's scheme (12) and Simpson's scheme (32) is now analyzed.

### 4.1. Symplectic property

A symplecticity proof is obtained by verifying that

$$ (33) \qquad \Phi^T \mathcal{J} \Phi = \mathcal{J}; \quad \mathcal{J} = \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & 0 \end{pmatrix}. $$

$\Phi$ corresponds to the scheme transformation matrix and characterizes a discrete time evolution of the system. $\mathcal{J}$ is sometimes referred to as the canonical matrix for Hamiltonian systems [30] and has the property that $\mathcal{J}^{-1} = \mathcal{J}^T = -\mathcal{J}$. When equation (33) holds, it means that $\Phi$ is an area-preserving transformation and that the scheme (12) is symplectic (see e. g. [18–20,31] for more details on this demonstration).

**Proposition 1.** *An implicit scheme of the type*

$$ \eta_{j+1} = A^{-1}B\eta_j; \quad \eta = (\mathbf{p}, \mathbf{q})^T, $$

is symplectic if

$$ (34) \qquad A = \begin{pmatrix} \mathbb{I}_n & -X \\ \mathbb{I}_n & Y \end{pmatrix}, \quad B = \begin{pmatrix} -\mathbb{I}_n & -X \\ \mathbb{I}_n & -Y \end{pmatrix}, $$

are square, partitioned, invertible matrices; and blocks X and Y are symmetric and positive-definite.

*Proof of Proposition 1.* Let us first explicit the transformation $A^{-1}$. Since A is square and partitioned, its inversion is done using auxiliary variables $\alpha$ and $\beta$. Let us establish that

$$ (35) \qquad A \begin{pmatrix} \mathbf{p} \\ \mathbf{q} \end{pmatrix} = \begin{pmatrix} \mathbf{p} - X\mathbf{q} \\ \mathbf{p} + Y\mathbf{q} \end{pmatrix} = \begin{pmatrix} \alpha \\ \beta \end{pmatrix} $$

Subtracting both equations above gives

$$ (36) \qquad \mathbf{q} = Z^{-1}(\beta - \alpha); \quad Z = X + Y. $$

Since Z is the sum of two symmetric, positive-definite matrices, it is invertible. Equation (36) is then substituted into the first equation of (35):

$$ (37) \qquad \mathbf{p} = (\mathbb{I}_n - XZ^{-1})\alpha + XZ^{-1}\beta. $$

Matrix A (34) is inverted in (36) and (37) as:

$$ A^{-1} = \begin{pmatrix} \mathbb{I}_n - XZ^{-1} & XZ^{-1} \\ -Z^{-1} & Z^{-1} \end{pmatrix}. $$
Then, it suffices to verify property (33) with $\Phi = A^{-1}B$. Thus,

$$ 
\begin{aligned}
\Phi^T \mathcal{J} \Phi &= B^T A^{-T} \mathcal{J} A^{-1} B = B^T \begin{pmatrix} -Z^{-T} & -(\mathbb{I}_n - Z^{-T} X^T) \\ Z^{-T} & -Z^{-T} X^T \end{pmatrix} A^{-1} B \\
&= B^T \begin{pmatrix} 0 & Z^{-1} \\ -Z^{-1} & 0 \end{pmatrix} B = \begin{pmatrix} Z^{-1} & Z^{-1} \\ -YZ^{-1} & XZ^{-1} \end{pmatrix} B \\
&= \begin{pmatrix} -Z^{-1} + Z^{-1} & -Z^{-1}(X+Y) \\ (Y+X)Z^{-1} & YZ^{-1}X - XZ^{-1}Y \end{pmatrix} \\
&= \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & (X^{-1}ZY^{-1})^{-1} - (Y^{-1}ZX^{-1})^{-1} \end{pmatrix} = \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & 0 \end{pmatrix} = \mathcal{J},
\end{aligned}
$$

because $X = X^T$, $Y = Y^T$, so $Z^{-1} = Z^{-T}$.

<figure><img src="image_1.png" /></figure>

### 4.2. Symplectic property of Newmark's scheme

**Proposition 2.** Newmark's scheme (12) is symplectic.

*Proof of Proposition 2.* In Newmark's scheme formulation (12), matrices $A_n$ and $B_n$ (11) are of the form (34); because $X_n$ and $Y_n$ (11) are symmetric and positive-definite. By Proposition 1, Newmark's scheme is symplectic. □

This confirms the classical result (e. g. [19]) on Newmark's scheme symplecticity.

### 4.3. Symplectic property of Simpson's scheme

To prove that Simpson's scheme is symplectic, we first need to prove that $A_s$ and $B_s$ have the structure of (34). For this, blocks $X_s$ and $Y_s$ are required to be symmetric and positive-definite.

**Proposition 3.** Matrix $Y_s$ (31) is symmetric.

*Proof of Proposition 3.* Since $M$ and $K$ are symmetric, $Y_s$ is symmetric if and only if its first term

$$W = KL^{-1}$$

is symmetric as well. $W$ is symmetric if $W^{-1}$ is symmetric. As

$$W^{-1} = \left( \mathbb{I}_n - \frac{h^2}{8} M^{-1} K \right) K^{-1} = K^{-1} - \frac{h^2}{8} M^{-1},$$

is symmetric, $Y_s$ is also symmetric. □

<figure><img src="image_2.png" /></figure>

For $Y_s$ to be positive-definite, the part $KL^{-1}$ must be positive-definite. Since $KL^{-1}$ is symmetric by Proposition 3, a condition on the step size $h$ is required.

Let us introduce the smallest and largest eigenvalues of matrices $M$ and $K$:

$$ 
\begin{aligned}
(38) \qquad & 0 < \mu \|\varphi\|^2 \le \varphi^T M \varphi \le m \|\varphi\|^2 \\
& 0 < \kappa \|\varphi\|^2 \le \varphi^T K \varphi \le k \|\varphi\|^2
\end{aligned}
$$

where $(\mu, \kappa)$ are the smallest and $(m, k)$ are the largest eigenvalues of matrices $M$ and $K$, respectively; $\varphi \ne 0$ is an eigenvector. Taking $M^{1/2}\varphi = \psi$ and then $K^{1/2}\varphi = \psi$, (38) becomes

$$ 
\begin{aligned}
\frac{1}{m} \|\psi\|^2 &\le \psi^T M \psi \le \frac{1}{\mu} \|\psi\|^2 \\
\frac{1}{k} \|\psi\|^2 &\le \psi^T K \psi \le \frac{1}{\kappa} \|\psi\|^2,
\end{aligned}
$$

and so

$$ (39) \qquad \psi^T (LK^{-1}) \psi \ge \left( \frac{1}{k} - \frac{h^2}{8} \frac{1}{\mu} \right) \|\psi\|^2. $$
Then, it suffices to verify property (33) with $\Phi = A^{-1}B$. Thus,

$$ 
\begin{aligned}
\Phi^T \mathcal{J} \Phi &= B^T A^{-T} \mathcal{J} A^{-1} B = B^T \begin{pmatrix} -Z^{-T} & -(\mathbb{I}_n - Z^{-T} X^T) \\ Z^{-T} & -Z^{-T} X^T \end{pmatrix} A^{-1} B \\
&= B^T \begin{pmatrix} 0 & Z^{-1} \\ -Z^{-1} & 0 \end{pmatrix} B = \begin{pmatrix} Z^{-1} & Z^{-1} \\ -YZ^{-1} & XZ^{-1} \end{pmatrix} B \\
&= \begin{pmatrix} -Z^{-1} + Z^{-1} & -Z^{-1}(X+Y) \\ (Y+X)Z^{-1} & YZ^{-1}X - XZ^{-1}Y \end{pmatrix} \\
&= \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & (X^{-1}ZY^{-1})^{-1} - (Y^{-1}ZX^{-1})^{-1} \end{pmatrix} = \begin{pmatrix} 0 & -\mathbb{I}_n \\ \mathbb{I}_n & 0 \end{pmatrix} = \mathcal{J},
\end{aligned}
$$

because $X = X^T$, $Y = Y^T$, so $Z^{-1} = Z^{-T}$.

<figure><img src="image_2.png" /></figure>

## 4.2. Symplectic property of Newmark's scheme

**Proposition 2.** Newmark's scheme (12) is symplectic.

*Proof of Proposition 2.* In Newmark's scheme formulation (12), matrices $A_n$ and $B_n$ (11) are of the form (34); because $X_n$ and $Y_n$ (11) are symmetric and positive-definite. By Proposition 1, Newmark's scheme is symplectic.

<figure><img src="image_3.png" /></figure>

This confirms the classical result (e. g. [19]) on Newmark's scheme symplecticity.

## 4.3. Symplectic property of Simpson's scheme

To prove that Simpson's scheme is symplectic, we first need to prove that $A_s$ and $B_s$ have the structure of (34). For this, blocks $X_s$ and $Y_s$ are required to be symmetric and positive-definite.

**Proposition 3.** Matrix $Y_s$ (31) is symmetric.

*Proof of Proposition 3.* Since $M$ and $K$ are symmetric, $Y_s$ is symmetric if and only if its first term

$$W = KL^{-1}$$

is symmetric as well. $W$ is symmetric if $W^{-1}$ is symmetric. As

$$W^{-1} = (\mathbb{I}_n - \frac{h^2}{8} M^{-1} K) K^{-1} = K^{-1} - \frac{h^2}{8} M^{-1},$$

is symmetric, $Y_s$ is also symmetric.

<figure><img src="image_4.png" /></figure>

For $Y_s$ to be positive-definite, the part $KL^{-1}$ must be positive-definite. Since $KL^{-1}$ is symmetric by Proposition 3, a condition on the step size $h$ is required.

Let us introduce the smallest and largest eigenvalues of matrices $M$ and $K$:

$$ 
\begin{aligned}
(38) \qquad & 0 < \mu \|\varphi\|^2 \le \varphi^T M \varphi \le m \|\varphi\|^2 \\
& 0 < \kappa \|\varphi\|^2 \le \varphi^T K \varphi \le k \|\varphi\|^2
\end{aligned}
$$

where $(\mu, \kappa)$ are the smallest and $(m, k)$ are the largest eigenvalues of matrices $M$ and $K$, respectively; $\varphi \ne 0$ is an eigenvector. Taking $M^{1/2}\varphi = \psi$ and then $K^{1/2}\varphi = \psi$, (38) becomes

$$ 
\begin{aligned}
\frac{1}{m} \|\psi\|^2 &\le \psi^T M \psi \le \frac{1}{\mu} \|\psi\|^2 \\
\frac{1}{k} \|\psi\|^2 &\le \psi^T K \psi \le \frac{1}{\kappa} \|\psi\|^2,
\end{aligned}
$$

and so

$$ (39) \qquad \psi^T (LK^{-1}) \psi \ge \left(\frac{1}{k} - \frac{h^2}{8} \frac{1}{\mu}\right) \|\psi\|^2. $$
The above expression is positive for

$$ (40) \qquad 0 < \frac{k}{\mu} h^2 < 8. $$

This inequality is a sufficient stability condition for Simpson's scheme. Let us remark that $k/\mu$ corresponds to the maximum eigenvalue of the dynamical matrix inverse $M^{-1}K$ and is associated with the maximum characteristic eigenfrequency of the system (see [22]) by

$$ \frac{k}{\mu} = \omega_{\max}^2. $$

The stability condition (40) can also be stated as

$$ 0 < \omega_{\max} h < 2\sqrt{2}. $$

This condition is similar to the stability condition characterizing the mono-dimensional case for Simpson's scheme [26].

**Proposition 4.** Matrix $Y_s$ (31) is positive-definite if $0 < \frac{k}{\mu}h^2 < 8$.

*Proof of Proposition 4.* When the condition (40) is met, equation (39) becomes

$$ \psi^T (LK^{-1}) \psi > 0, \quad \forall \psi \neq \mathbf{0}, $$

and $LK^{-1}$ is positive-definite. Therefore, $KL^{-1}$ is also positive-definite, and recalling Proposition 3, it is symmetric. Consequently,

$$ \psi^T Y_s \psi > 0, \quad \forall \psi \neq \mathbf{0}, $$

and $Y_s$ is positive-definite.

<figure><img src="image_3.png" /></figure>

Now, only the positive-definiteness of block $X_s$ (31) remains to be proven.

**Proposition 5.** Matrix $X_s$ (31) is positive-definite if $0 < \frac{k}{\mu}h^2 < 8$.

*Proof of Proposition 5.* From inequalities (38),

$$ \varphi^T M \varphi - \frac{h^2}{12} \varphi^T K \varphi \geq \left( \mu - \frac{h^2}{12} k \right) \|\varphi\|^2, $$

and substituting the condition (40) on the first term of the right-hand side of the above inequality,

$$ \varphi^T \left( M - \frac{h^2}{12} K \right) \varphi \geq \frac{\mu}{3} \|\varphi\|^2 > 0, \quad \forall \varphi \neq \mathbf{0}. $$

Therefore, matrix $X_s$ is positive-definite.

<figure><img src="image_4.png" /></figure>

**Proposition 6.** Simpson's scheme (32) is symplectic.

*Proof of Proposition 6.* For the second variant of Simpson's scheme (32), matrices $A_s$ and $B_s$ (31) are of the form (34); because $X_s$ is symmetric and positive-definite by Proposition 5, and $Y_s$ is symmetric and positive-definite by Propositions 3 and 4.

<figure><img src="image_5.png" /></figure>

These results prove that the proposed Simpson's scheme is symplectic.

### 4.4. Conservation of a discrete quadratic form

Symplectic integrators usually do not preserve the energy quantity. This has been summarized in [32] and outlined in [19]. The goal is to verify that Simpson's scheme preserves some quadratic form. It is required that some quadratic function $\phi(\mathbf{p}, \mathbf{q})$ verifies

$$ \phi(\mathbf{p}_{j+1}, \mathbf{q}_{j+1}) = \phi(\mathbf{p}_j, \mathbf{q}_j), $$

where $(\mathbf{p}_{j+1}, \mathbf{q}_{j+1})$ and $(\mathbf{p}_j, \mathbf{q}_j)$ satisfy the dynamics of Simpson's scheme (31)-(32).
The above expression is positive for

$$ (40) \qquad 0 < \frac{k}{\mu} h^2 < 8. $$

This inequality is a sufficient stability condition for Simpson's scheme. Let us remark that $k/\mu$ corresponds to the maximum eigenvalue of the dynamical matrix inverse $M^{-1}K$ and is associated with the maximum characteristic eigenfrequency of the system (see [22]) by

$$ \frac{k}{\mu} = \omega_{\max}^2. $$

The stability condition (40) can also be stated as

$$ 0 < \omega_{\max} h < 2\sqrt{2}. $$

This condition is similar to the stability condition characterizing the mono-dimensional case for Simpson's scheme [26].

**Proposition 4.** Matrix $Y_s$ (31) is positive-definite if $0 < \frac{k}{\mu}h^2 < 8$.

*Proof of Proposition 4.* When the condition (40) is met, equation (39) becomes

$$ \psi^T (LK^{-1}) \psi > 0, \quad \forall \psi \neq \mathbf{0}, $$

and $LK^{-1}$ is positive-definite. Therefore, $KL^{-1}$ is also positive-definite, and recalling Proposition 3, it is symmetric. Consequently,

$$ \psi^T Y_s \psi > 0, \quad \forall \psi \neq \mathbf{0}, $$

and $Y_s$ is positive-definite.

<figure><img src="image_5.png" /></figure>

Now, only the positive-definiteness of block $X_s$ (31) remains to be proven.

**Proposition 5.** Matrix $X_s$ (31) is positive-definite if $0 < \frac{k}{\mu}h^2 < 8$.

*Proof of Proposition 5.* From inequalities (38),

$$ \varphi^T M \varphi - \frac{h^2}{12} \varphi^T K \varphi \geq \left( \mu - \frac{h^2}{12} k \right) \|\varphi\|^2, $$

and substituting the condition (40) on the first term of the right-hand side of the above inequality,

$$ \varphi^T \left( M - \frac{h^2}{12} K \right) \varphi \geq \frac{\mu}{3} \|\varphi\|^2 > 0, \quad \forall \varphi \neq \mathbf{0}. $$

Therefore, matrix $X_s$ is positive-definite.

<figure><img src="image_6.png" /></figure>

**Proposition 6.** Simpson's scheme (32) is symplectic.

*Proof of Proposition 6.* For the second variant of Simpson's scheme (32), matrices $A_s$ and $B_s$ (31) are of the form (34); because $X_s$ is symmetric and positive-definite by Proposition 5, and $Y_s$ is symmetric and positive-definite by Propositions 3 and 4.

<figure><img src="image_7.png" /></figure>

These results prove that the proposed Simpson's scheme is symplectic.

### 4.4. Conservation of a discrete quadratic form

Symplectic integrators usually do not preserve the energy quantity. This has been summarized in [32] and outlined in [19]. The goal is to verify that Simpson's scheme preserves some quadratic form. It is required that some quadratic function $\phi(\mathbf{p}, \mathbf{q})$ verifies

$$ \phi(\mathbf{p}_{j+1}, \mathbf{q}_{j+1}) = \phi(\mathbf{p}_j, \mathbf{q}_j), $$

where $(\mathbf{p}_{j+1}, \mathbf{q}_{j+1})$ and $(\mathbf{p}_j, \mathbf{q}_j)$ satisfy the dynamics of Simpson's scheme (31)-(32).
**Proposition 7.** Given an implicit scheme of the type

$$
(41) \qquad \eta_{j+1} = A^{-1}B \eta_j; \qquad \eta = (\mathbf{p}, \mathbf{q})^T,
$$

where

$$
A = \begin{pmatrix} \mathbb{I}_n & -X \\ \mathbb{I}_n & Y \end{pmatrix}, \quad B = \begin{pmatrix} -\mathbb{I}_n & -X \\ \mathbb{I}_n & -Y \end{pmatrix},
$$

are square, partitioned, invertible matrices, and blocks X and Y are symmetric and positive-definite;
there exists a quadratic form

$$
(42) \qquad \phi(\mathbf{p}, \mathbf{q}) = \frac{1}{2} \mathbf{p}^T \zeta \mathbf{p} + \frac{1}{2} \mathbf{q}^T \zeta \mathbf{q},
$$

which is conserved if

$$
\zeta = (X + Y)^{-1}; \quad \zeta = (X^{-1} + Y^{-1})^{-1}.
$$

Proof of Proposition 7. Let us expand the scheme (41):

$$
\mathbf{p}_{j+1} - X\mathbf{q}_{j+1} = -\mathbf{p}_j - X\mathbf{q}_j,
$$

$$
\mathbf{p}_{j+1} + X\mathbf{q}_{j+1} = \mathbf{p}_j - Y\mathbf{q}_j.
$$

The above can also be written as

$$
\mathbf{p}_{j+1} + \mathbf{p}_j = X (\mathbf{q}_{j+1} - \mathbf{q}_j)
$$

$$
\mathbf{p}_{j+1} - \mathbf{p}_j = -Y (\mathbf{q}_{j+1} + \mathbf{q}_j).
$$

Therefore, by multiplying $\zeta$ by the first equation above on the left, and by the second equation above on the right,

$$
(43) \qquad
\begin{aligned}
(\mathbf{p}_{j+1} + \mathbf{p}_j)^T \zeta (\mathbf{p}_{j+1} - \mathbf{p}_j) &= (X (\mathbf{q}_{j+1} - \mathbf{q}_j))^T \zeta (-Y) (\mathbf{q}_{j+1} + \mathbf{q}_j) \\
&= -(\mathbf{q}_{j+1} - \mathbf{q}_j)^T X \zeta Y (\mathbf{q}_{j+1} + \mathbf{q}_j).
\end{aligned}
$$

Since

$$
(X\zeta Y)^{-1} = Y^{-1}(X + Y)X^{-1} = X^{-1} + Y^{-1} = \zeta^{-1}
$$

is symmetric and positive-definite, it is deduced that $\zeta = X\zeta Y$ is symmetric and positive-definite. Following from (43),

$$
\mathbf{p}_{j+1}{}^T \zeta \mathbf{p}_{j+1} + \mathbf{q}_{j+1}{}^T \zeta \mathbf{q}_{j+1} = \mathbf{p}_j{}^T \zeta \mathbf{p}_j + \mathbf{q}_j{}^T \zeta \mathbf{q}_j,
$$

and the property is proven since $\phi(\mathbf{p}_{j+1}, \mathbf{q}_{j+1}) = \phi(\mathbf{p}_j, \mathbf{q}_j)$. $\square$

<figure><img src="image_6.png" /></figure>

By Proposition 7, and condition (40), Simpson's scheme is conditionally stable.

5. LINEAR DOUBLE PENDULUM MODEL AND EXACT SOLUTION

This section presents a case study for subsequent numerical experiments.

5.1. Lagrangian

Let us model the system depicted by Figure 1. It is a two-degree-of-freedom dynamical
system composed of two masses ($m_1, m_2$) linked together by two massless thin rigid rods of
respective fixed lengths ($l_1, l_2$). Each joint articulates the system with one rotational degree of
freedom. Masses coordinates are given by

$$
(x_1, y_1) = (l_1 \sin q_1, -l_1 \cos q_1)
$$

$$
(x_2, y_2) = (l_1 \sin q_1 + l_2 \sin q_2, -l_1 \cos q_1 - l_2 \cos q_2),
$$

and their velocities are obtained by time differentiation considering that $q_i = q_i(t)$. The system
kinetic energy is then given by

$$
T = \frac{1}{2} m_1 (x_1^2 + y_1^2) + \frac{1}{2} m_2 (x_2^2 + y_2^2),
$$

where an overdot indicates time differentiation. Potential energy is calculated as

$$
V = -m_1 g l_1 \cos q_1 - m_2 g (l_1 \cos q_1 + l_2 \cos q_2),
$$
**Proposition 7.** Given an implicit scheme of the type

$$
(41) \qquad \eta_{j+1} = A^{-1}B \eta_j; \qquad \eta = (\mathbf{p}, \mathbf{q})^T,
$$

where

$$
A = \begin{pmatrix} \mathbb{I}_n & -X \\ \mathbb{I}_n & Y \end{pmatrix}, \quad B = \begin{pmatrix} -\mathbb{I}_n & -X \\ \mathbb{I}_n & -Y \end{pmatrix},
$$

are square, partitioned, invertible matrices, and blocks X and Y are symmetric and positive-definite;
there exists a quadratic form

$$
(42) \qquad \phi(\mathbf{p}, \mathbf{q}) = \frac{1}{2} \mathbf{p}^T \zeta \mathbf{p} + \frac{1}{2} \mathbf{q}^T \zeta \mathbf{q},
$$

which is conserved if

$$
\zeta = (X + Y)^{-1}; \quad \zeta = (X^{-1} + Y^{-1})^{-1}.
$$

Proof of Proposition 7. Let us expand the scheme (41):

$$
\mathbf{p}_{j+1} - X\mathbf{q}_{j+1} = -\mathbf{p}_j - X\mathbf{q}_j,
$$

$$
\mathbf{p}_{j+1} + X\mathbf{q}_{j+1} = \mathbf{p}_j - Y\mathbf{q}_j.
$$

The above can also be written as

$$
\mathbf{p}_{j+1} + \mathbf{p}_j = X (\mathbf{q}_{j+1} - \mathbf{q}_j)
$$

$$
\mathbf{p}_{j+1} - \mathbf{p}_j = -Y (\mathbf{q}_{j+1} + \mathbf{q}_j).
$$

Therefore, by multiplying $\zeta$ by the first equation above on the left, and by the second equation above on the right,

$$
(43) \qquad
\begin{aligned}
(\mathbf{p}_{j+1} + \mathbf{p}_j)^T \zeta (\mathbf{p}_{j+1} - \mathbf{p}_j) &= (X (\mathbf{q}_{j+1} - \mathbf{q}_j))^T \zeta (-Y) (\mathbf{q}_{j+1} + \mathbf{q}_j) \\
&= -(\mathbf{q}_{j+1} - \mathbf{q}_j)^T X \zeta Y (\mathbf{q}_{j+1} + \mathbf{q}_j).
\end{aligned}
$$

Since

$$
(X\zeta Y)^{-1} = Y^{-1}(X + Y)X^{-1} = X^{-1} + Y^{-1} = \zeta^{-1}
$$

is symmetric and positive-definite, it is deduced that $\zeta = X\zeta Y$ is symmetric and positive-definite. Following from (43),

$$
\mathbf{p}_{j+1}{}^T \zeta \mathbf{p}_{j+1} + \mathbf{q}_{j+1}{}^T \zeta \mathbf{q}_{j+1} = \mathbf{p}_j{}^T \zeta \mathbf{p}_j + \mathbf{q}_j{}^T \zeta \mathbf{q}_j,
$$

and the property is proven since $\phi(\mathbf{p}_{j+1}, \mathbf{q}_{j+1}) = \phi(\mathbf{p}_j, \mathbf{q}_j)$.

<figure><img src="image_8.png" /></figure>

By Proposition 7, and condition (40), Simpson's scheme is conditionally stable.

5. LINEAR DOUBLE PENDULUM MODEL AND EXACT SOLUTION

This section presents a case study for subsequent numerical experiments.

**5.1. Lagrangian**

Let us model the system depicted by Figure 1. It is a two-degree-of-freedom dynamical
system composed of two masses ($m_1, m_2$) linked together by two massless thin rigid rods of
respective fixed lengths ($l_1, l_2$). Each joint articulates the system with one rotational degree of
freedom. Masses coordinates are given by

$$
(x_1, y_1) = (l_1 \sin q_1, -l_1 \cos q_1)
$$

$$
(x_2, y_2) = (l_1 \sin q_1 + l_2 \sin q_2, -l_1 \cos q_1 - l_2 \cos q_2),
$$

and their velocities are obtained by time differentiation considering that $q_i = q_i(t)$. The system
kinetic energy is then given by

$$
T = \frac{1}{2} m_1 (x_1^2 + y_1^2) + \frac{1}{2} m_2 (x_2^2 + y_2^2),
$$

where an overdot indicates time differentiation. Potential energy is calculated as

$$
V = -m_1 g l_1 \cos q_1 - m_2 g (l_1 \cos q_1 + l_2 \cos q_2),
$$
and finally the system Lagrangian $\mathcal{L} = T - V$ can be explicit as

$$ (44) \qquad \mathcal{L} = \frac{1}{2}(m_1 + m_2)l_1^2 \dot{q}_1^2 + \frac{1}{2}m_2 l_2^2 \dot{q}_2^2 + m_2 l_1 l_2 \dot{q}_1 \dot{q}_2 \cos(q_1 - q_2) \\ \qquad \qquad + (m_1 + m_2)g l_1 \cos q_1 + m_2 g l_2 \cos q_2. $$

<figure><img src="image_7.png" /><figcaption>FIGURE 1. Double pendulum system subject to the gravity action. The system is composed of two masses ($m_1, m_2$) linked together by two massless thin rigid rods of respective fixed lengths ($l_1, l_2$). Each joint articulates the system with one rotational degree of freedom. Masses are located by the generalized coordinates $\mathbf{q} = (q_1, q_2)$.</figcaption></figure>

Small oscillations take place when $q_i(t)$ are small and around the stable equilibrium. This equilibrium corresponds to the system's resting position when it is aligned with the vertical axis pointing downwards. Such motions can be described by linear equations. In this situation, the Lagrangian (44) takes a simpler form provided that the following approximations take place:

$$ (45) \qquad \begin{aligned} \cos q_1 &\approx 1 - \frac{q_1^2}{2}; \\ \cos q_2 &\approx 1 - \frac{q_2^2}{2}; \\ \cos(q_1 - q_2) &\approx 1 - \frac{(q_1 - q_2)^2}{2}. \end{aligned} $$

Using (45), the linear form $L_L$ of Lagrangian (44) becomes

$$ (46) \qquad L_L = \frac{1}{2}(m_1 + m_2)l_1(l_1\dot{q}_1^2 + 2g - gq_1^2) + \frac{1}{2}m_2l_2(l_2\dot{q}_2^2 + 2l_1\dot{q}_1\dot{q}_2 + 2g - gq_2^2), $$

where the second term of $\cos(q_1 - q_2)$ approximation in (45) vanishes when multiplying the product $\dot{q}_1\dot{q}_2$. Generalized momenta are defined as

$$ p_i = \frac{\partial L_L}{\partial \dot{q}_i}. $$

According to the Lagrangian (46), we have

$$ \begin{aligned} p_1 &= (m_1 + m_2)l_1^2 \dot{q}_1 + m_2 l_1 l_2 \dot{q}_2, \\ p_2 &= m_2 l_2 (l_1 \dot{q}_1 + l_2 \dot{q}_2). \end{aligned} $$

Motion equations are then obtained by applying Euler-Lagrange equations $\frac{d}{dt} \frac{\partial L_L}{\partial \dot{q}_i} - \frac{\partial L_L}{\partial q_i} = 0$:
and finally the system Lagrangian $\mathcal{L} = T - V$ can be explicit as

$$ (44) \qquad \mathcal{L} = \frac{1}{2}(m_1 + m_2)l_1^2 \dot{q}_1^2 + \frac{1}{2}m_2 l_2^2 \dot{q}_2^2 + m_2 l_1 l_2 \dot{q}_1 \dot{q}_2 \cos(q_1 - q_2) \\ \qquad \qquad + (m_1 + m_2)g l_1 \cos q_1 + m_2 g l_2 \cos q_2. $$

<figure><img src="image_9.png" /><figcaption>FIGURE 1. Double pendulum system subject to the gravity action. The system is composed of two masses ($m_1, m_2$) linked together by two massless thin rigid rods of respective fixed lengths ($l_1, l_2$). Each joint articulates the system with one rotational degree of freedom. Masses are located by the generalized coordinates $\mathbf{q} = (q_1, q_2)$.</figcaption></figure>

Small oscillations take place when $q_i(t)$ are small and around the stable equilibrium. This equilibrium corresponds to the system's resting position when it is aligned with the vertical axis pointing downwards. Such motions can be described by linear equations. In this situation, the Lagrangian (44) takes a simpler form provided that the following approximations take place:

$$ (45) \qquad \begin{aligned} \cos q_1 &\approx 1 - \frac{q_1^2}{2}; \\ \cos q_2 &\approx 1 - \frac{q_2^2}{2}; \\ \cos(q_1 - q_2) &\approx 1 - \frac{(q_1 - q_2)^2}{2}. \end{aligned} $$

Using (45), the linear form $L_L$ of Lagrangian (44) becomes

$$ (46) \qquad L_L = \frac{1}{2}(m_1 + m_2)l_1(l_1 \dot{q}_1^2 + 2g - gq_1^2) + \frac{1}{2}m_2 l_2(l_2 \dot{q}_2^2 + 2l_1 \dot{q}_1 \dot{q}_2 + 2g - gq_2^2), $$

where the second term of $\cos(q_1 - q_2)$ approximation in (45) vanishes when multiplying the product $\dot{q}_1 \dot{q}_2$. Generalized momenta are defined as

$$ p_i = \frac{\partial L_L}{\partial \dot{q}_i}. $$

According to the Lagrangian (46), we have

$$ \begin{aligned} p_1 &= (m_1 + m_2)l_1^2 \dot{q}_1 + m_2 l_1 l_2 \dot{q}_2, \\ p_2 &= m_2 l_2 (l_1 \dot{q}_1 + l_2 \dot{q}_2). \end{aligned} $$

Motion equations are then obtained by applying Euler-Lagrange equations $\frac{d}{dt} \frac{\partial L_L}{\partial \dot{q}_i} - \frac{\partial L_L}{\partial q_i} = 0$:
$$ (47) \qquad \begin{aligned} & (m_1 + m_2)l_1^2 \ddot{q}_1 + m_2 l_1 l_2 \ddot{q}_2 + (m_1 + m_2)g l_1 q_1 = 0, \\ & m_2 l_1 l_2 \ddot{q}_1 + m_2 l_2^2 \ddot{q}_2 + m_2 g l_2 q_2 = 0. \end{aligned} $$

## 5.2. Exact solution

Equations (47) can also be established as a linear system of the form

$$ (48) \qquad M\ddot{q} + Kq = 0. $$

where

$$ M = \begin{pmatrix} (m_1 + m_2)l_1^2 & m_2 l_1 l_2 \\ m_2 l_1 l_2 & m_2 l_2^2 \end{pmatrix}; \quad q = \begin{pmatrix} q_1 \\ q_2 \end{pmatrix}; \quad K = \begin{pmatrix} (m_1 + m_2)g l_1 & 0 \\ 0 & m_2 g l_2 \end{pmatrix}. $$

The general solution of (48) is of the form

$$ q(t) = \mathrm{Re} \left( \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \end{bmatrix} \right) \cdot e^{i\omega t}, $$

where $\mathbf{x}_1$ and $\mathbf{x}_2$ are eigenvectors and $\omega$ denotes the oscillation frequency. Two characteristic frequencies ($\omega_1, \omega_2$) are determined by the solution of the auxiliary equation $\det(K - \omega^2 M) = 0$:

$$ (m_1 + m_2)g^2 - (m_1 + m_2)g(l_1 + l_2)\omega^2 + m_1 l_1 l_2 \omega^4 = 0. $$

Let us focus on the case where

$$ l_1 = l_2 = l. $$

In this particular case, oscillation frequencies are given by

$$ (49) \qquad \omega_{1,2} = \omega_0 \sqrt{\left(1 + \mu_r \pm \sqrt{\mu_r(1 + \mu_r)}\right)}, $$

with a mass ratio $\mu_r = m_2/m_1$ and frequency $\omega_0 = \sqrt{g/l}$.

Eigenvectors $\mathbf{x}_1$ and $\mathbf{x}_2$ are then obtained by solving $(K - \omega_i^2 M)\mathbf{x}_i = 0$ for $i = 1$ and $i = 2$,

$$ m_1 l \begin{bmatrix} (1 + \mu_r)(g - \omega_i^2 l) & -\omega_i^2 \mu_r l \\ -\omega_i^2 \mu_r l & \mu_r(g - \omega_i^2 \mu_r l) \end{bmatrix} \mathbf{x}_i = 0. $$

Solving the above system gives

$$ (50) \qquad \mathbf{x}_1 = \begin{bmatrix} 1 \\ -\sqrt{\frac{1+\mu_r}{\mu_r}} \end{bmatrix}, \quad \mathbf{x}_2 = \begin{bmatrix} 1 \\ \sqrt{\frac{1+\mu_r}{\mu_r}} \end{bmatrix}. $$

Finally, using (50), the general solution of (48) (or (47)) can be established as

$$ (51) \qquad \mathbf{q}(t) = c_1 \mathbf{x}_1 \cos(\omega_1 t + \varphi_1) + c_2 \mathbf{x}_2 \cos(\omega_2 t + \varphi_2), $$

where constants ($c_1, c_2, \varphi_1, \varphi_2$) are given by the chosen initial conditions on positions and velocities.

## 6. SIMULATION RESULTS

We will now assess the precision and convergence of our proposed integrator, previously described in section 3. It will be compared with Newmark's symplectic scheme, described in section 2. Some results obtained with Runge-Kutta's explicit fourth-order integrator, described in [3] and labeled as "RK4" throughout the rest of the document, are also given. Note that a thorough comparison with this classical integrator is beyond the scope of our contribution. Results are provided for reference since RK4 is among the most popular methods available. For benchmarking purposes, we apply these methods to the solution of the linear double pendulum (depicted by Figure 1), which has an exact solution described in the previous section 5.

The results presented in this section are for a simulated motion of this linearized double pendulum. Computations were carried out using Wolfram's Mathematica software [33]. Figure plots were then created using exported data with the "pgfplots" package from LATEX. Table 1 specifies the constants and initial conditions used for all of our simulations. Using these values
$$ (47) \qquad \begin{aligned} & (m_1 + m_2)l_1^2 \ddot{q}_1 + m_2 l_1 l_2 \ddot{q}_2 + (m_1 + m_2)g l_1 q_1 = 0, \\ & m_2 l_1 l_2 \ddot{q}_1 + m_2 l_2^2 \ddot{q}_2 + m_2 g l_2 q_2 = 0. \end{aligned} $$

## 5.2. Exact solution

Equations (47) can also be established as a linear system of the form

$$ (48) \qquad M\ddot{q} + Kq = 0. $$

where

$$ M = \begin{pmatrix} (m_1 + m_2)l_1^2 & m_2 l_1 l_2 \\ m_2 l_1 l_2 & m_2 l_2^2 \end{pmatrix}; \quad q = \begin{pmatrix} q_1 \\ q_2 \end{pmatrix}; \quad K = \begin{pmatrix} (m_1 + m_2)g l_1 & 0 \\ 0 & m_2 g l_2 \end{pmatrix}. $$

The general solution of (48) is of the form

$$ q(t) = \mathrm{Re} \left( \begin{bmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \end{bmatrix} \right) \cdot e^{i\omega t}, $$

where $\mathbf{x}_1$ and $\mathbf{x}_2$ are eigenvectors and $\omega$ denotes the oscillation frequency. Two characteristic frequencies ($\omega_1, \omega_2$) are determined by the solution of the auxiliary equation $\det(K - \omega^2 M) = 0$:

$$ (m_1 + m_2)g^2 - (m_1 + m_2)g(l_1 + l_2)\omega^2 + m_1 l_1 l_2 \omega^4 = 0. $$

Let us focus on the case where

$$ l_1 = l_2 = l. $$

In this particular case, oscillation frequencies are given by

$$ (49) \qquad \omega_{1,2} = \omega_0 \sqrt{\left(1 + \mu_r \pm \sqrt{\mu_r(1 + \mu_r)}\right)}, $$

with a mass ratio $\mu_r = m_2/m_1$ and frequency $\omega_0 = \sqrt{g/l}$.

Eigenvectors $\mathbf{x}_1$ and $\mathbf{x}_2$ are then obtained by solving $(K - \omega_i^2 M)\mathbf{x}_i = 0$ for $i = 1$ and $i = 2$,

$$ m_1 l \begin{bmatrix} (1 + \mu_r)(g - \omega_i^2 l) & -\omega_i^2 \mu_r l \\ -\omega_i^2 \mu_r l & \mu_r(g - \omega_i^2 \mu_r l) \end{bmatrix} \mathbf{x}_i = 0. $$

Solving the above system gives

$$ (50) \qquad \mathbf{x}_1 = \begin{bmatrix} 1 \\ -\sqrt{\frac{1+\mu_r}{\mu_r}} \end{bmatrix}, \quad \mathbf{x}_2 = \begin{bmatrix} 1 \\ \sqrt{\frac{1+\mu_r}{\mu_r}} \end{bmatrix}. $$

Finally, using (50), the general solution of (48) (or (47)) can be established as

$$ (51) \qquad \mathbf{q}(t) = c_1 \mathbf{x}_1 \cos(\omega_1 t + \varphi_1) + c_2 \mathbf{x}_2 \cos(\omega_2 t + \varphi_2), $$

where constants ($c_1, c_2, \varphi_1, \varphi_2$) are given by the chosen initial conditions on positions and velocities.

## 6. SIMULATION RESULTS

We will now assess the precision and convergence of our proposed integrator, previously described in section 3. It will be compared with Newmark's symplectic scheme, described in section 2. Some results obtained with Runge-Kutta's explicit fourth-order integrator, described in [3] and labeled as "RK4" throughout the rest of the document, are also given. Note that a thorough comparison with this classical integrator is beyond the scope of our contribution. Results are provided for reference since RK4 is among the most popular methods available. For benchmarking purposes, we apply these methods to the solution of the linear double pendulum (depicted by Figure 1), which has an exact solution described in the previous section 5.

The results presented in this section are for a simulated motion of this linearized double pendulum. Computations were carried out using Wolfram's Mathematica software [33]. Figure plots were then created using exported data with the "pgfplots" package from LATEX. Table 1 specifies the constants and initial conditions used for all of our simulations. Using these values
and following (51) with null initial phases ($\varphi_1, \varphi_2$), the exact solution that serves as the main reference in our comparisons is

$$ \mathbf{q}_{\text{ex}}(t) = \frac{\pi}{12} \left( \frac{-\cos(\omega_1 t) + \cos(\omega_2 t)}{\sqrt{2}} \right), $$

where $\omega_1$ and $\omega_2$ are given by (49).

<table><caption>TABLE 1. Constants and initial conditions used for numerical simulations.</caption><thead><tr><th>Constants</th><th>Initial conditions</th></tr></thead><tbody><tr><td>$\mu_r = 1$</td><td>$\mathbf{q}(0) = (0, \pi/6)^T \text{ rad}$</td></tr><tr><td>$\omega_0 = 2\pi \text{ s}^{-1}$</td><td>$\mathbf{p}(0) = (0, 0)^T \text{ kg m}^2 \text{ s}^{-1}$</td></tr></tbody></table>

Frequency $\omega_0$ is used to show results in terms of an oscillation period $\bar{t}$ such that

$$ \bar{t} = \frac{1}{\omega_0}. $$

Therefore, both the total simulation duration $T$ and step size $h$ are given in terms of $\bar{t}$. It is to be noted that the presented results from Simpson's scheme were obtained using the second variant (see section 3.5). Hence the absence of the middle value at each interpolation interval. However, both variants provided lead to the same result at each node.

## 6.1. Configuration parameters and generalized momenta

We begin by comparing the configuration parameter solutions $\mathbf{q}$ obtained with the proposed Simpson's rule-based variational integrator, against those given by Newmark's method. The proposed integrator uses quadratic finite elements for interpolation and Simpson's rule (see section 3). It is expected to be more precise than Newmark's method which uses a centered finite difference and the midpoint integration rule (see section 2). Figure 2 shows the configuration parameters provided by each method, compared against the exact solution, during one period $\bar{t}$. Simpson's integrator is already more precise than Newmark's scheme. Runge-Kutta's solution is also close to the exact one but not as much as Simpson's solution.

<figure><img src="image_8.png" /><figcaption>FIGURE 2. Configuration parameters $\mathbf{q}$ evolution for the linear double pendulum. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$. Simpson's integrator tracks the exact solution with more precision than Newmark's method and Runge-Kutta's integrator.</figcaption></figure>

Figure 2 shows that Simpson's integrator is more precise than both Newmark's and Runge-Kutta's integrators on a short simulation ($T = 1 \bar{t}$). However, Simpson's solutions correctly
and following (51) with null initial phases ($\varphi_1, \varphi_2$), the exact solution that serves as the main reference in our comparisons is

$$ \mathbf{q}_{\text{ex}}(t) = \frac{\pi}{12} \left( \frac{-\cos(\omega_1 t) + \cos(\omega_2 t)}{\sqrt{2}} \right), $$

where $\omega_1$ and $\omega_2$ are given by (49).

<table><caption>TABLE 1. Constants and initial conditions used for numerical simulations.</caption><thead><tr><th>Constants</th><th>Initial conditions</th></tr></thead><tbody><tr><td>$\mu_r = 1$</td><td>$\mathbf{q}(0) = (0, \pi/6)^T \text{ rad}$</td></tr><tr><td>$\omega_0 = 2\pi \text{ s}^{-1}$</td><td>$\mathbf{p}(0) = (0, 0)^T \text{ kg m}^2 \text{ s}^{-1}$</td></tr></tbody></table>

Frequency $\omega_0$ is used to show results in terms of an oscillation period $\bar{t}$ such that

$$ \bar{t} = \frac{1}{\omega_0}. $$

Therefore, both the total simulation duration $T$ and step size $h$ are given in terms of $\bar{t}$. It is to be noted that the presented results from Simpson's scheme were obtained using the second variant (see section 3.5). Hence the absence of the middle value at each interpolation interval. However, both variants provided lead to the same result at each node.

## 6.1. Configuration parameters and generalized momenta

We begin by comparing the configuration parameter solutions $\mathbf{q}$ obtained with the proposed Simpson's rule-based variational integrator, against those given by Newmark's method. The proposed integrator uses quadratic finite elements for interpolation and Simpson's rule (see section 3). It is expected to be more precise than Newmark's method which uses a centered finite difference and the midpoint integration rule (see section 2). Figure 2 shows the configuration parameters provided by each method, compared against the exact solution, during one period $\bar{t}$. Simpson's integrator is already more precise than Newmark's scheme. Runge-Kutta's solution is also close to the exact one but not as much as Simpson's solution.

<figure><img src="image_10.png" /><figcaption>FIGURE 2. Configuration parameters $\mathbf{q}$ evolution for the linear double pendulum. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$. Simpson's integrator tracks the exact solution with more precision than Newmark's method and Runge-Kutta's integrator.</figcaption></figure>

Figure 2 shows that Simpson's integrator is more precise than both Newmark's and Runge-Kutta's integrators on a short simulation ($T = 1 \bar{t}$). However, Simpson's solutions correctly
follow the exact ones for longer simulations on both the configuration parameters and generalized momenta, as shown by Figure 3.

<figure><img src="image_9.png" /><figcaption>FIGURE 3. Configuration parameters **q** and generalized momenta **p** evolution for the linear double pendulum during ten periods. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$. Simpson's solutions correctly follow the exact solution for longer simulations.</figcaption></figure>
follow the exact ones for longer simulations on both the configuration parameters and generalized momenta, as shown by Figure 3.

<figure><img src="image_11.png" /><figcaption>FIGURE 3. Configuration parameters **q** and generalized momenta **p** evolution for the linear double pendulum during ten periods. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$. Simpson's solutions correctly follow the exact solution for longer simulations.</figcaption></figure>
Newmark's integrator precision can be increased by refining the step size. With $h = 0.01 \bar{t}$, solutions improve but still deviate from the exact solution after a couple of periods. Simpson's solutions correctly follow the exact solution for longer simulations.

## 6.2. Phase portraits

With a step size of $h = 0.1 \bar{t}$, Newmark's solutions deviations are particularly visible when tracing the motion phase portrait. Figure 4 shows the exact phase portraits topped by both Newmark's and Simpson's solutions. Notice how Simpson's phase portrait clearly follows the exact one throughout the motion. Total simulation time was limited to $T = 3 \bar{t}$ for visualization purposes.

<figure><img src="image_10.png" /><figcaption>FIGURE 4. Phase portraits evolution for the linear double pendulum. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$ and three periods are shown ($T = 3 \bar{t}$). Simpson's phase portrait clearly follows the exact one.</figcaption></figure>

## 6.3. Energy conservation

The following function gives the system energy

$$H(\mathbf{p}, \mathbf{q}) = \frac{1}{2} \mathbf{p}^T M^{-1} \mathbf{p} + \frac{1}{2} \mathbf{q}^T K \mathbf{q}.$$

It has been previously observed (e. g. [19]) that Newmark's integrator exactly preserves the system energy. This is not the case for our proposed integrator based on Simpson's rule, and is characteristic of most symplectic methods [19, 32]. In the case of Simpson's scheme, the second equation of (25), introduces the small and vanishing quantity $-\frac{h^2}{12} K \left( \frac{\mathbf{q}_{i+1} - \mathbf{q}_i}{h} \right)$ into the discrete momentum equation. Consequently, one could assume that the exact system energy may not be conserved but a good energy behavior can be expected as outlined in [32].

Figure 5 shows that Simpson's solutions lead to a non-conserved energy $H(\mathbf{p}, \mathbf{q})$. Nevertheless, the maximum relative error with respect to the initial value is extremely small even for $h = 0.1 \bar{t}$, as evaluated values are in the order of $10^{-3}$. Notice that the energy error from Simpson's solutions does not grow with time. Instead, it oscillates in a bounded fashion. Note that the error drops by four orders of magnitude when dividing the step size by ten (see Figure 5).
Newmark's integrator precision can be increased by refining the step size. With $h = 0.01 \bar{t}$, solutions improve but still deviate from the exact solution after a couple of periods. Simpson's solutions correctly follow the exact solution for longer simulations.

## 6.2. Phase portraits

With a step size of $h = 0.1 \bar{t}$, Newmark's solutions deviations are particularly visible when tracing the motion phase portrait. Figure 4 shows the exact phase portraits topped by both Newmark's and Simpson's solutions. Notice how Simpson's phase portrait clearly follows the exact one throughout the motion. Total simulation time was limited to $T = 3 \bar{t}$ for visualization purposes.

<figure><img src="image_12.png" /><figcaption>FIGURE 4. Phase portraits evolution for the linear double pendulum. Initial conditions are specified in Table 1. The step size is fixed as $h = 0.1 \bar{t}$ and three periods are shown ($T = 3 \bar{t}$). Simpson's phase portrait clearly follows the exact one.</figcaption></figure>

## 6.3. Energy conservation

The following function gives the system energy

$$H(\mathbf{p}, \mathbf{q}) = \frac{1}{2} \mathbf{p}^T M^{-1} \mathbf{p} + \frac{1}{2} \mathbf{q}^T K \mathbf{q}.$$

It has been previously observed (e. g. [19]) that Newmark's integrator exactly preserves the system energy. This is not the case for our proposed integrator based on Simpson's rule, and is characteristic of most symplectic methods [19, 32]. In the case of Simpson's scheme, the second equation of (25), introduces the small and vanishing quantity $-\frac{h^2}{12} K \left(\frac{\mathbf{q}_{i+1}-\mathbf{q}_i}{h}\right)$ into the discrete momentum equation. Consequently, one could assume that the exact system energy may not be conserved but a good energy behavior can be expected as outlined in [32].

Figure 5 shows that Simpson's solutions lead to a non-conserved energy $H(\mathbf{p}, \mathbf{q})$. Nevertheless, the maximum relative error with respect to the initial value is extremely small even for $h = 0.1 \bar{t}$, as evaluated values are in the order of $10^{-3}$. Notice that the energy error from Simpson's solutions does not grow with time. Instead, it oscillates in a bounded fashion. Note that the error drops by four orders of magnitude when dividing the step size by ten (see Figure 5).
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 5. As expected, the classical RK4 integrator does not preserve the system energy. Relative error grows with simulation length. Newmark's integrator exactly preserves the system energy. Simpson's integrator does not, but the relative error is extremely small. Notice that such an error does not grow with time but remains bounded. Relative error drops by four orders of magnitude when dividing the step size by ten, showcasing the quality of the proposed integrator and its good energy behavior.</figcaption></figure>

Proposition 7 shows that Simpson's scheme preserves a quadratic form given by the function $\phi(\mathbf{p}, \mathbf{q})$ (42). Matrices $\zeta_s$ and $\zeta_s$ (where subscript s stands for Simpson) are according to Proposition 7 as

$$ \zeta_s = (X_s + Y_s)^{-1} = \left[ \frac{2}{h} M + \frac{h}{3} KL^{-1} \right]^{-1}; $$

$$ \zeta_s = (X_s^{-1} + Y_s^{-1})^{-1} = \left[ \left( \frac{2}{h} M - \frac{h}{6} K \right)^{-1} + \frac{3}{h} \left( KL^{-1} + \frac{1}{2} K \right)^{-1} \right]^{-1}. $$

Figure 6 plots the absolute error on function $\phi(\mathbf{p}, \mathbf{q})$ (42), by Simpson's scheme. The absolute error with respect to the initial value is minimal, in the order of $10^{-15}$, and may come from accumulated rounding errors. Note that this absolute error magnitude changes very little when refining the step size.

<figure><img src="image_13.png" /><figcaption>FIGURE 6. Simpson's scheme preserves the quadratic form $\phi(\mathbf{p}, \mathbf{q})$ (42). Absolute errors are minimal and may come from accumulated rounding errors.</figcaption></figure>

## 6.4. Convergence

The error $e(t) = \mathbf{q}(t) - \mathbf{q}_{\text{ex}}(t)$ and its convergence rate is measured following the prescriptions found in [16]. The schemes precision was evaluated using an $\ell^\infty$ error norm $\|e\|_\infty = \|\mathbf{q} - \mathbf{q}_{\text{ex}}\|_\infty = \sup_n |\mathbf{q}_n - \mathbf{q}_{\text{ex}_n}|$. Several simulations were held for decreasing values
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 5. As expected, the classical RK4 integrator does not preserve the system energy. Relative error grows with simulation length. Newmark's integrator exactly preserves the system energy. Simpson's integrator does not, but the relative error is extremely small. Notice that such an error does not grow with time but remains bounded. Relative error drops by four orders of magnitude when dividing the step size by ten, showcasing the quality of the proposed integrator and its good energy behavior.</figcaption></figure>

Proposition 7 shows that Simpson's scheme preserves a quadratic form given by the function $\phi(\mathbf{p}, \mathbf{q})$ (42). Matrices $\zeta_s$ and $\zeta_s$ (where subscript s stands for Simpson) are according to Proposition 7 as

$$ \zeta_s = (X_s + Y_s)^{-1} = \left[ \frac{2}{h} M + \frac{h}{3} KL^{-1} \right]^{-1}; $$

$$ \zeta_s = (X_s^{-1} + Y_s^{-1})^{-1} = \left[ \left( \frac{2}{h} M - \frac{h}{6} K \right)^{-1} + \frac{3}{h} \left( KL^{-1} + \frac{1}{2} K \right)^{-1} \right]^{-1}. $$

Figure 6 plots the absolute error on function $\phi(\mathbf{p}, \mathbf{q})$ (42), by Simpson's scheme. The absolute error with respect to the initial value is minimal, in the order of $10^{-15}$, and may come from accumulated rounding errors. Note that this absolute error magnitude changes very little when refining the step size.

<figure><img src="image_15.png" /><figcaption>FIGURE 6. Simpson's scheme preserves the quadratic form $\phi(\mathbf{p}, \mathbf{q})$ (42). Absolute errors are minimal and may come from accumulated rounding errors.</figcaption></figure>

## 6.4. Convergence

The error $e(t) = \mathbf{q}(t) - \mathbf{q}_{\text{ex}}(t)$ and its convergence rate is measured following the prescriptions found in [16]. The schemes precision was evaluated using an $\ell^\infty$ error norm $\|e\|_\infty = \|\mathbf{q} - \mathbf{q}_{\text{ex}}\|_\infty = \sup_n |\mathbf{q}_n - \mathbf{q}_{\text{ex}_n}|$. Several simulations were held for decreasing values
of $h$ between $h = 0.1 \bar{t}$ and $h = 0.001 \bar{t}$. The $\|e\|_{\infty}$ norm was calculated for each case. These errors were plotted in Figure 7, on the logarithmic scale.

Convergence rates are expressed as the power of the step size. These rates correspond to the slope of the error logarithm, as a function of the logarithm of $h$ (see Figure 7). These trials confirm previous analyses on Newmark's method [17, 19]: it is second-order convergent. Unsurprisingly, Runge-Kutta's integrator is fourth-order convergent. The results also confirm the analysis held in [26] on the convergence rate of Simpson's scheme: it is fourth-order convergent. This rate is two degrees higher than the order of the chosen quadratic interpolation. This is known as superconvergence and is closely related to the mesh uniformity [16].

<figure><img src="image_14.png" /><figcaption>FIGURE 7. Integrators convergence. The $\|e\|_{\infty}$ norm was calculated for several simulations. Each simulation used a fixed step size, which was decreased from $h = 0.1 \bar{t}$ to $h = 0.001 \bar{t}$. The convergence order corresponds to the slope of the error norm logarithm regression line. Runge-Kutta's classical integrator convergence is in $h^4$ as expected (fourth order). Newmark's integrator convergence is in the order of $h^2$ (second order). Simpson's integrator convergence is in the order of $h^4$ (fourth order).</figcaption></figure>

An important question is if convergence rates hold with growing simulation lengths $T$. Table 2 shows the convergence order evolution of Newmark's and Simpson's schemes with a growing simulation length, using a step size of $h = 0.1 \bar{t}$. It can be observed that Newmark's scheme convergence order decays to zero for a 1000-period simulation. RK4's convergence rate also degrades as the simulation duration increases, although not as much as Newmark's
of $h$ between $h = 0.1 \bar{t}$ and $h = 0.001 \bar{t}$. The $\|e\|_{\infty}$ norm was calculated for each case. These errors were plotted in Figure 7, on the logarithmic scale.

Convergence rates are expressed as the power of the step size. These rates correspond to the slope of the error logarithm, as a function of the logarithm of $h$ (see Figure 7). These trials confirm previous analyses on Newmark's method [17, 19]: it is second-order convergent. Unsurprisingly, Runge-Kutta's integrator is fourth-order convergent. The results also confirm the analysis held in [26] on the convergence rate of Simpson's scheme: it is fourth-order convergent. This rate is two degrees higher than the order of the chosen quadratic interpolation. This is known as superconvergence and is closely related to the mesh uniformity [16].

<figure><img src="image_16.png" /><figcaption>FIGURE 7. Integrators convergence. The $\|e\|_{\infty}$ norm was calculated for several simulations. Each simulation used a fixed step size, which was decreased from $h = 0.1 \bar{t}$ to $h = 0.001 \bar{t}$. The convergence order corresponds to the slope of the error norm logarithm regression line. Runge-Kutta's classical integrator convergence is in $h^4$ as expected (fourth order). Newmark's integrator convergence is in the order of $h^2$ (second order). Simpson's integrator convergence is in the order of $h^4$ (fourth order).</figcaption></figure>

An important question is if convergence rates hold with growing simulation lengths $T$. Table 2 shows the convergence order evolution of Newmark's and Simpson's schemes with a growing simulation length, using a step size of $h = 0.1 \bar{t}$. It can be observed that Newmark's scheme convergence order decays to zero for a 1000-period simulation. RK4's convergence rate also degrades as the simulation duration increases, although not as much as Newmark's
method. Simpson's scheme preserves its convergence order for higher simulation periods. Error norms $||e||_\infty$ are shown explicitly. Simpson's scheme loses precision according to one order of magnitude, each time the simulation length is multiplied by ten. Table 2 exposes a normal numerical behavior of the analyzed schemes since errors accumulate over long simulations.

<table><caption>TABLE 2. Convergence order with respect to simulation length for motion simulations held on a linearized double pendulum (see Figure 1). Initial conditions are specified in Table 1. Error norms $||e||_\infty$ for momenta p and states q increase with simulation length T. Newmark's scheme converges to zero as T increases. RK4's convergence rate also decays with an increasing simulation length. Simpson's scheme preserves its convergence rate for higher simulation durations.</caption><thead><tr><th rowspan="2">Simulation length T</th><th rowspan="2">Number of meshes</th><th colspan="3">||e||<sub>∞</sub> error norm values</th><th rowspan="2">Convergence order</th></tr><tr><th>10</th><th>20</th><th>40</th></tr></thead><tbody><tr><td rowspan="6">1f</td><td>Newmark p</td><td>7.51 × 10<sup>-2</sup></td><td>2.30 × 10<sup>-2</sup></td><td>6.06 × 10<sup>-3</sup></td><td>1.81</td></tr><tr><td>RK4 p</td><td>1.39 × 10<sup>-2</sup></td><td>8.00 × 10<sup>-4</sup></td><td>5.40 × 10<sup>-5</sup></td><td>4.01</td></tr><tr><td>Simpson p</td><td>6.40 × 10<sup>-4</sup></td><td>4.16 × 10<sup>-5</sup></td><td>2.57 × 10<sup>-6</sup></td><td><b>3.98</b></td></tr><tr><td>Newmark q</td><td>3.42 × 10<sup>-1</sup></td><td>9.61 × 10<sup>-2</sup></td><td>2.51 × 10<sup>-2</sup></td><td>1.88</td></tr><tr><td>RK4 q</td><td>4.83 × 10<sup>-2</sup></td><td>3.40 × 10<sup>-3</sup></td><td>2.00 × 10<sup>-4</sup></td><td>3.91</td></tr><tr><td>Simpson q</td><td>2.01 × 10<sup>-3</sup></td><td>1.41 × 10<sup>-4</sup></td><td>8.76 × 10<sup>-6</sup></td><td><b>3.92</b></td></tr><tr><td rowspan="6">10f</td><td>Newmark p</td><td>2.73 × 10<sup>-1</sup></td><td>2.06 × 10<sup>-1</sup></td><td>7.82 × 10<sup>-1</sup></td><td>0.90</td></tr><tr><td>RK4 p</td><td>8.22 × 10<sup>-2</sup></td><td>9.90 × 10<sup>-3</sup></td><td>1.37 × 10<sup>-2</sup></td><td>1.29</td></tr><tr><td>Simpson p</td><td>7.20 × 10<sup>-3</sup></td><td>4.33 × 10<sup>-4</sup></td><td>2.68 × 10<sup>-5</sup></td><td><b>4.03</b></td></tr><tr><td>Newmark q</td><td>6.94 × 10<sup>-1</sup></td><td>6.57 × 10<sup>-1</sup></td><td>2.44 × 10<sup>-1</sup></td><td>0.75</td></tr><tr><td>RK4 q</td><td>2.84 × 10<sup>-1</sup></td><td>3.29 × 10<sup>-2</sup></td><td>1.57 × 10<sup>-2</sup></td><td>2.09</td></tr><tr><td>Simpson q</td><td>2.35 × 10<sup>-2</sup></td><td>1.41 × 10<sup>-3</sup></td><td>9.06 × 10<sup>-5</sup></td><td><b>4.01</b></td></tr><tr><td rowspan="6">100f</td><td>Newmark p</td><td>5.21 × 10<sup>-1</sup></td><td>4.92 × 10<sup>-1</sup></td><td>2.23 × 10<sup>-1</sup></td><td>0.61</td></tr><tr><td>RK4 p</td><td>1.08 × 10<sup>-1</sup></td><td>7.86 × 10<sup>-2</sup></td><td>6.40 × 10<sup>-3</sup></td><td>2.03</td></tr><tr><td>Simpson p</td><td>7.05 × 10<sup>-2</sup></td><td>4.39 × 10<sup>-3</sup></td><td>2.72 × 10<sup>-4</sup></td><td><b>4.01</b></td></tr><tr><td>Newmark q</td><td>1.02</td><td>9.64 × 10<sup>-1</sup></td><td>6.65 × 10<sup>-1</sup></td><td>0.31</td></tr><tr><td>RK4 q</td><td>3.28 × 10<sup>-1</sup></td><td>2.650 × 10<sup>-1</sup></td><td>2.16 × 10<sup>-2</sup></td><td>1.96</td></tr><tr><td>Simpson q</td><td>2.37 × 10<sup>-1</sup></td><td>1.47 × 10<sup>-2</sup></td><td>9.14 × 10<sup>-4</sup></td><td><b>4.01</b></td></tr><tr><td rowspan="6">1000f</td><td>Newmark p</td><td>5.45 × 10<sup>-1</sup></td><td>5.51 × 10<sup>-1</sup></td><td>5.48 × 10<sup>-1</sup></td><td>0.00</td></tr><tr><td>RK4 p</td><td>3.26 × 10<sup>-1</sup></td><td>1.19 × 10<sup>-1</sup></td><td>5.95 × 10<sup>-2</sup></td><td>1.23</td></tr><tr><td>Simpson p</td><td>1.90 × 10<sup>-1</sup></td><td>4.38 × 10<sup>-2</sup></td><td>2.74 × 10<sup>-3</sup></td><td><b>3.06</b></td></tr><tr><td>Newmark q</td><td>1.02</td><td>1.03</td><td>1.03</td><td>0.01</td></tr><tr><td>RK4 q</td><td>5.81 × 10<sup>-1</sup></td><td>3.97 × 10<sup>-1</sup></td><td>2.00 × 10<sup>-1</sup></td><td>0.77</td></tr><tr><td>Simpson q</td><td>6.38 × 10<sup>-1</sup></td><td>1.47 × 10<sup>-1</sup></td><td>9.22 × 10<sup>-3</sup></td><td><b>3.06</b></td></tr></tbody></table>

## 7. CONCLUDING REMARKS AND PERSPECTIVES

In this contribution, Newmark's method has been recalled. It is a widely used integrator in certain fields of the engineering sciences and it is symplectic. This method has been used for benchmarking purposes in our work, where an alternative variational integrator based on Simpson's integration rule has been proposed.
method. Simpson's scheme preserves its convergence order for higher simulation periods. Error norms $||e||_\infty$ are shown explicitly. Simpson's scheme loses precision according to one order of magnitude, each time the simulation length is multiplied by ten. Table 2 exposes a normal numerical behavior of the analyzed schemes since errors accumulate over long simulations.

<table><caption>TABLE 2. Convergence order with respect to simulation length for motion simulations held on a linearized double pendulum (see Figure 1). Initial conditions are specified in Table 1. Error norms $||e||_\infty$ for momenta p and states q increase with simulation length T. Newmark's scheme converges to zero as T increases. RK4's convergence rate also decays with an increasing simulation length. Simpson's scheme preserves its convergence rate for higher simulation durations.</caption><thead><tr><th rowspan="2">Simulation length T</th><th rowspan="2">Number of meshes</th><th rowspan="2"></th><th colspan="3">||e||<sub>∞</sub> error norm values</th><th rowspan="2">Convergence order</th></tr><tr><th>10</th><th>20</th><th>40</th></tr></thead><tbody><tr><td rowspan="6">1<span style="text-decoration: overline;">t</span></td><td>Newmark</td><td>p</td><td>7.51 × 10<sup>-2</sup></td><td>2.30 × 10<sup>-2</sup></td><td>6.06 × 10<sup>-3</sup></td><td>1.81</td></tr><tr><td>RK4</td><td>p</td><td>1.39 × 10<sup>-2</sup></td><td>8.00 × 10<sup>-4</sup></td><td>5.40 × 10<sup>-5</sup></td><td>4.01</td></tr><tr><td>Simpson</td><td>p</td><td>6.40 × 10<sup>-4</sup></td><td>4.16 × 10<sup>-5</sup></td><td>2.57 × 10<sup>-6</sup></td><td><strong>3.98</strong></td></tr><tr><td>Newmark</td><td>q</td><td>3.42 × 10<sup>-1</sup></td><td>9.61 × 10<sup>-2</sup></td><td>2.51 × 10<sup>-2</sup></td><td>1.88</td></tr><tr><td>RK4</td><td>q</td><td>4.83 × 10<sup>-2</sup></td><td>3.40 × 10<sup>-3</sup></td><td>2.00 × 10<sup>-4</sup></td><td>3.91</td></tr><tr><td>Simpson</td><td>q</td><td>2.01 × 10<sup>-3</sup></td><td>1.41 × 10<sup>-4</sup></td><td>8.76 × 10<sup>-6</sup></td><td><strong>3.92</strong></td></tr><tr><td rowspan="6">10<span style="text-decoration: overline;">t</span></td><td>Newmark</td><td>p</td><td>2.73 × 10<sup>-1</sup></td><td>2.06 × 10<sup>-1</sup></td><td>7.82 × 10<sup>-1</sup></td><td>0.90</td></tr><tr><td>RK4</td><td>p</td><td>8.22 × 10<sup>-2</sup></td><td>9.90 × 10<sup>-3</sup></td><td>1.37 × 10<sup>-2</sup></td><td>1.29</td></tr><tr><td>Simpson</td><td>p</td><td>7.20 × 10<sup>-3</sup></td><td>4.33 × 10<sup>-4</sup></td><td>2.68 × 10<sup>-5</sup></td><td><strong>4.03</strong></td></tr><tr><td>Newmark</td><td>q</td><td>6.94 × 10<sup>-1</sup></td><td>6.57 × 10<sup>-1</sup></td><td>2.44 × 10<sup>-1</sup></td><td>0.75</td></tr><tr><td>RK4</td><td>q</td><td>2.84 × 10<sup>-1</sup></td><td>3.29 × 10<sup>-2</sup></td><td>1.57 × 10<sup>-2</sup></td><td>2.09</td></tr><tr><td>Simpson</td><td>q</td><td>2.35 × 10<sup>-2</sup></td><td>1.41 × 10<sup>-3</sup></td><td>9.06 × 10<sup>-5</sup></td><td><strong>4.01</strong></td></tr><tr><td rowspan="6">100<span style="text-decoration: overline;">t</span></td><td>Newmark</td><td>p</td><td>5.21 × 10<sup>-1</sup></td><td>4.92 × 10<sup>-1</sup></td><td>2.23 × 10<sup>-1</sup></td><td>0.61</td></tr><tr><td>RK4</td><td>p</td><td>1.08 × 10<sup>-1</sup></td><td>7.86 × 10<sup>-2</sup></td><td>6.40 × 10<sup>-3</sup></td><td>2.03</td></tr><tr><td>Simpson</td><td>p</td><td>7.05 × 10<sup>-2</sup></td><td>4.39 × 10<sup>-3</sup></td><td>2.72 × 10<sup>-4</sup></td><td><strong>4.01</strong></td></tr><tr><td>Newmark</td><td>q</td><td>1.02</td><td>9.64 × 10<sup>-1</sup></td><td>6.65 × 10<sup>-1</sup></td><td>0.31</td></tr><tr><td>RK4</td><td>q</td><td>3.28 × 10<sup>-1</sup></td><td>2.650 × 10<sup>-1</sup></td><td>2.16 × 10<sup>-2</sup></td><td>1.96</td></tr><tr><td>Simpson</td><td>q</td><td>2.37 × 10<sup>-1</sup></td><td>1.47 × 10<sup>-2</sup></td><td>9.14 × 10<sup>-4</sup></td><td><strong>4.01</strong></td></tr><tr><td rowspan="6">1000<span style="text-decoration: overline;">t</span></td><td>Newmark</td><td>p</td><td>5.45 × 10<sup>-1</sup></td><td>5.51 × 10<sup>-1</sup></td><td>5.48 × 10<sup>-1</sup></td><td>0.00</td></tr><tr><td>RK4</td><td>p</td><td>3.26 × 10<sup>-1</sup></td><td>1.19 × 10<sup>-1</sup></td><td>5.95 × 10<sup>-2</sup></td><td>1.23</td></tr><tr><td>Simpson</td><td>p</td><td>1.90 × 10<sup>-1</sup></td><td>4.38 × 10<sup>-2</sup></td><td>2.74 × 10<sup>-3</sup></td><td><strong>3.06</strong></td></tr><tr><td>Newmark</td><td>q</td><td>1.02</td><td>1.03</td><td>1.03</td><td>0.01</td></tr><tr><td>RK4</td><td>q</td><td>5.81 × 10<sup>-1</sup></td><td>3.97 × 10<sup>-1</sup></td><td>2.00 × 10<sup>-1</sup></td><td>0.77</td></tr><tr><td>Simpson</td><td>q</td><td>6.38 × 10<sup>-1</sup></td><td>1.47 × 10<sup>-1</sup></td><td>9.22 × 10<sup>-3</sup></td><td><strong>3.06</strong></td></tr></tbody></table>

## 7. CONCLUDING REMARKS AND PERSPECTIVES

In this contribution, Newmark's method has been recalled. It is a widely used integrator in certain fields of the engineering sciences and it is symplectic. This method has been used for benchmarking purposes in our work, where an alternative variational integrator based on Simpson's integration rule has been proposed.
The presented Simpson's numerical scheme applies to the case of multiple degrees of freedom systems with quadratic Lagrangians. It has been formulated linearly with partitioned matrices. The method proves to be symplectic, as demonstrated with a proof that applies to both Newmark's and Simpson's scheme. A sufficient stability condition on the step size is given and it is also proved that the proposed method preserves a certain quadratic form at each time step. Simpson's scheme is therefore conditionally stable.

Numerical trials on a linearized double pendulum have confirmed that the method is fourth-order accurate on both the states and generalized momenta. Numerical evaluations reveal that this convergence order is preserved for long simulations. The proposed method succeeds in predicting the evolution of dynamical systems characterized by quadratic Lagrangians.

An important extension of this work is the treatment of nonlinear multi-degrees-of-freedom systems. In such a configuration, the middle value of the internal interpolation cannot be eliminated. This generalization should enable more applications of the proposed method, relating to Hamiltonian systems. Therefore, this is a natural objective for future developments and is currently under study. An important question relates to noise presence in matrices *M* and *K*. How would this affect the symplectic integrator? This question is relevant in the context of nonlinear dynamical systems. It shall be the object of future developments as well.

A particular subject of interest relating to differential equations is the role of discrete symmetries. The analysis of discrete symmetries has many applications for finding solutions to differential equations. They can simplify a numerical scheme as advocated in [34]. A description of finding discrete symmetries of differential equations has been given in [35]. A discrete symmetry analysis could lead to an improved symplectic integrator and is a future direction for our work.

An improved nonlinear Simpson's variational integrator could find its application to simulate complex nonlinear mechanisms. Some application examples could involve a system of synchronized pendulums [36]; the discrete optimal control of robotic systems [37]; the modal analysis of dynamical systems [38]; the motion analysis of multibody systems evolving in fluid environments [39]; or the motion prediction of sliding rods [24] and soft robots [25].

**Acknowledgments.** This research was conducted while the corresponding author held a one-year visiting position at the Laboratoire des Sciences du Numérique de Nantes (LS2N) of the Institut Mines-Télécom Atlantique (IMTA) of Nantes, France; under the authorization of CONAHCYT. The corresponding author acknowledges the support received from IMTA and CONAHCYT.

## REFERENCES

[1] Herman H. Goldstine. *A History of Numerical Analysis from the 16th through the 19th Century*. Studies in the History of Mathematics and Physical Sciences. Springer, New York, NY, December 1977. doi:10.1007/978-1-4684-9472-3

[2] Hans-Joachim Albinus. The mathematical tourist. *The Mathematical Intelligencer*, 24(3):50–58, 2002. doi:10.1007/BF03024733.

[3] Endre Süli and David F. Mayers. *An Introduction to Numerical Analysis*. Cambridge University Press, Cambridge, UK, 2003.

[4] Albert L. Herman and Bruce A. Conway. Direct optimization using collocation based on high-order gauss-lobatto quadrature rules. *Journal of Guidance, Control, and Dynamics*, 19(3):592–599, 1996. doi:10.2514/3.21662.

[5] Matthew Kelly. An introduction to trajectory optimization: How to do your own direct collocation. *SIAM Review*, 59(4):849–904, 2017. doi:10.1137/16M1062569.

[6] Tomasz Blaszczyk, Jaroslaw Siedlecki, and Mariusz Ciesielski. Numerical algorithms for approximation of fractional integral operators based on quadratic interpolation. *Mathematical Methods in the Applied Sciences*, 41(9):3345–3355, 2018. doi:10.1002/mma.4828.

[7] Aneela Sabir and Mujeeb ur Rehman. A numerical method based on quadrature rules for $\psi$-fractional differential equations. *Journal of Computational and Applied Mathematics*, 419:114684, 2023. doi:10.1016/j.cam.2022.114684.

[8] Zeinab Salamat Mostaghim, Behrouz Parsa Moghaddam, and Hossein Samimi Haghgozar. Numerical simulation of fractional-order dynamical systems in noisy environments. *Computational and Applied Mathematics*, 37(5):6433–6447, 2018. doi:10.1007/s40314-018-0698-z.

[9] M. Fakharany, Mahmoud M. El-Borai, and M. A. Abu Ibrahim. Numerical analysis of finite difference schemes arising from time-memory partial integro-differential equations. *Frontiers in Applied Mathematics and Statistics*, 8, 2022. doi:10.3389/fams.2022.1055071.

[10] Vladimir Igorevich Arnold. *Mathematical Methods of Classical Mechanics*. Springer New York, New York, 2nd edition, 1989.

[11] Richard Phillips Feynman and Albert Roach Hibbs. *Quantum mechanics and path integrals*. International series in pure and applied physics. McGraw-Hill, New York, NY, 1965.
The presented Simpson's numerical scheme applies to the case of multiple degrees of freedom systems with quadratic Lagrangians. It has been formulated linearly with partitioned matrices. The method proves to be symplectic, as demonstrated with a proof that applies to both Newmark's and Simpson's scheme. A sufficient stability condition on the step size is given and it is also proved that the proposed method preserves a certain quadratic form at each time step. Simpson's scheme is therefore conditionally stable.

Numerical trials on a linearized double pendulum have confirmed that the method is fourth-order accurate on both the states and generalized momenta. Numerical evaluations reveal that this convergence order is preserved for long simulations. The proposed method succeeds in predicting the evolution of dynamical systems characterized by quadratic Lagrangians.

An important extension of this work is the treatment of nonlinear multi-degrees-of-freedom systems. In such a configuration, the middle value of the internal interpolation cannot be eliminated. This generalization should enable more applications of the proposed method, relating to Hamiltonian systems. Therefore, this is a natural objective for future developments and is currently under study. An important question relates to noise presence in matrices *M* and *K*. How would this affect the symplectic integrator? This question is relevant in the context of nonlinear dynamical systems. It shall be the object of future developments as well.

A particular subject of interest relating to differential equations is the role of discrete symmetries. The analysis of discrete symmetries has many applications for finding solutions to differential equations. They can simplify a numerical scheme as advocated in [34]. A description of finding discrete symmetries of differential equations has been given in [35]. A discrete symmetry analysis could lead to an improved symplectic integrator and is a future direction for our work.

An improved nonlinear Simpson's variational integrator could find its application to simulate complex nonlinear mechanisms. Some application examples could involve a system of synchronized pendulums [36]; the discrete optimal control of robotic systems [37]; the modal analysis of dynamical systems [38]; the motion analysis of multibody systems evolving in fluid environments [39]; or the motion prediction of sliding rods [24] and soft robots [25].

**Acknowledgments.** This research was conducted while the corresponding author held a one-year visiting position at the Laboratoire des Sciences du Numérique de Nantes (LS2N) of the Institut Mines-Télécom Atlantique (IMTA) of Nantes, France; under the authorization of CONAHCYT. The corresponding author acknowledges the support received from IMTA and CONAHCYT.

## REFERENCES

[1] Herman H. Goldstine. *A History of Numerical Analysis from the 16th through the 19th Century*. Studies in the History of Mathematics and Physical Sciences. Springer, New York, NY, December 1977. doi:10.1007/978-1-4684-9472-3

[2] Hans-Joachim Albinus. The mathematical tourist. *The Mathematical Intelligencer*, 24(3):50–58, 2002. doi:10.1007/BF03024733.

[3] Endre Süli and David F. Mayers. *An Introduction to Numerical Analysis*. Cambridge University Press, Cambridge, UK, 2003.

[4] Albert L. Herman and Bruce A. Conway. Direct optimization using collocation based on high-order gauss-lobatto quadrature rules. *Journal of Guidance, Control, and Dynamics*, 19(3):592–599, 1996. doi:10.2514/3.21662.

[5] Matthew Kelly. An introduction to trajectory optimization: How to do your own direct collocation. *SIAM Review*, 59(4):849–904, 2017. doi:10.1137/16M1062569.

[6] Tomasz Blaszczyk, Jaroslaw Siedlecki, and Mariusz Ciesielski. Numerical algorithms for approximation of fractional integral operators based on quadratic interpolation. *Mathematical Methods in the Applied Sciences*, 41(9):3345–3355, 2018. doi:10.1002/mma.4828.

[7] Aneela Sabir and Mujeeb ur Rehman. A numerical method based on quadrature rules for $\psi$-fractional differential equations. *Journal of Computational and Applied Mathematics*, 419:114684, 2023. doi:10.1016/j.cam.2022.114684.

[8] Zeinab Salamat Mostaghim, Behrouz Parsa Moghaddam, and Hossein Samimi Haghgozar. Numerical simulation of fractional-order dynamical systems in noisy environments. *Computational and Applied Mathematics*, 37(5):6433–6447, 2018. doi:10.1007/s40314-018-0698-z.

[9] M. Fakharany, Mahmoud M. El-Borai, and M. A. Abu Ibrahim. Numerical analysis of finite difference schemes arising from time-memory partial integro-differential equations. *Frontiers in Applied Mathematics and Statistics*, 8, 2022. doi:10.3389/fams.2022.1055071.

[10] Vladimir Igorevich Arnold. *Mathematical Methods of Classical Mechanics*. Springer New York, New York, 2nd edition, 1989.

[11] Richard Phillips Feynman and Albert Roach Hibbs. *Quantum mechanics and path integrals*. International series in pure and applied physics. McGraw-Hill, New York, NY, 1965.
[12] Jean-Maris Souriau. *Structure des systèmes dynamiques*. Maîtrises de mathématiques. Dunod, Paris, France, 1970.

[13] Adi Nusser and Enzo Branchini. On the least action principle in cosmology. *Monthly Notices of the Royal Astronomical Society*, 313(3):587-595, 04 2000. doi:10.1046/j.1365-8711.2000.03261.x.

[14] Alberto Rojo and Anthony Bloch. *The Principle of Least Action: History and Physics*. Cambridge University Press, 2018.

[15] Richard Courant. Variational methods for the solution of problems of equilibrium and vibrations. *Bulletin of the American Mathematical Society*, 49(1):1-23, 1943.

[16] Gregoire Allaire and Alan Craig. *Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, UK, 2007.

[17] Nathan M. Newmark. A Method of Computation for Structural Dynamics. *Journal of the Engineering Mechanics Division*, 85(3):67-94, 1959. doi:10.1061/JMCEA3.0000098.

[18] Jesús María Sanz-Serna. Symplectic integrators for Hamiltonian problems: an overview. *Acta Numerica*, 1:243-286, 1992. doi:10.1017/S0962492900002282.

[19] C. Kane, J. E. Marsden, M. Ortiz, and M. West. Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems. *International Journal for Numerical Methods in Engineering*, 49(10):1295-1325, 2000. doi:10.1002/1097-0207(20001210)49:10<1295>:AID-NME993>3.0.CO;2-W.

[20] Ernst Hairer, Gerhard Wanner, and Christian Lubich. *Geometric Numerical Integration*. Springer Series in Computational Mathematics. Springer-Verlag Berlin Heidelberg, 2006. doi:10.1007/3-540-30666-8.

[21] René de Vogelaere. Methods of integration which preserve the contact transformation property of the hamilton equations. 1956. doi:10.7274/24733422.v1.

[22] Michel Géradin and Daniel J. Rixen. *Mechanical Vibrations: Theory and Application to Structural Dynamics*. John Wiley & Sons, Chichester, UK, 2015.

[23] Anil K. Chopra. *Dynamics of Structures: Theory and Applications to Earthquake Engineering*. Prentice-Hall International Series in Civil Engineering and Engineering Mechanics. Pearson Education Limited, Harlow, UK, 5th edition, 2020.

[24] Frédéric Boyer, Vincent Lebastard, Fabien Candelier, and Federico Renda. Extended hamilton's principle applied to geometrically exact kirchhoff sliding rods. *Journal of Sound and Vibration*, 516:116511, 2022. doi:10.1016/j.jsv.2021.116511.

[25] Frédéric Boyer, Andrea Gotelli, Philipp Tempel, Vincent Lebastard, Federico Renda, and Sébastien Briot. Implicit time-integration simulation of robots with rigid bodies and cosserat rods based on a newton-euler recursive algorithm. *IEEE Transactions on Robotics*, 40:677-696, 2024. doi:10.1109/TR0.2023.3334647.

[26] François Dubois and Juan Antonio Rojas-Quintero. A variational symplectic scheme based on simpson's quadrature. In Frank Nielsen and Frédéric Barbaresco, editors, *Geometric Science of Information*, pages 22-31, Cham, 2023. Springer Nature Switzerland. doi:10.1007/978-3-031-38299-4_3.

[27] François Dubois and Juan Antonio Rojas-Quintero. Simpson's quadrature for a nonlinear variational symplectic scheme. In Emmanuel Franck, Jürgen Fuhrmann, Victor Michel-Dansac, and Laurent Navoret, editors, *Finite Volumes for Complex Applications X—Volume 2, Hyperbolic and Related Problems*, pages 83-92, Cham, 2023. Springer Nature Switzerland. doi:10.1007/978-3-031-40860-1_9.

[28] Pierre-Arnaud Raviart and Jean-Marie Thomas. *Introduction à l'analyse numérique des équations aux dérivées partielles*. Collection Mathématiques appliquées pour la maîtrise. Masson, Paris, France, 1983.

[29] C. Gear. Simultaneous numerical solution of differential-algebraic equations. *IEEE Transactions on Circuit Theory*, 18(1):89-95, 1971. doi:10.1109/TCT.1971.1083221.

[30] Monique Chyba, Ernst Hairer, and Gilles Vilmart. The role of symplectic integrators in optimal control. *Optimal Control Applications and Methods*, 30(4):367-382, 2009. doi:10.1002/oca.855.

[31] Dina Razafindralandy, Vladimir Salnikov, Aziz Hamdouni, and Ahmad Deeb. Some robust integrators for large time dynamics. *Advanced Modeling and Simulation in Engineering Sciences*, 6(1):5, 2019. doi:10.1186/s40323-019-0130-2.

[32] Ge Zhong and Jerrold E. Marsden. Lie-poisson hamilton-jacobi theory and lie-poisson integrators. *Physics Letters A*, 133(3):134-139, 1988. doi:10.1016/0375-9601(88)90773-6.

[33] Wolfram Research, Inc. *Mathematica*, Version 12.3. Champaign, IL, 2021.

[34] Nail Hairullovich Ibragimov. *Elementary Lie Group Analysis and Ordinary Differential Equations*. Wiley, Chichester, UK, 1999.

[35] Peter E. Hydon. *Symmetry Methods for Differential Equations: A Beginner's Guide*. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2000.

[36] Sergey Oktyabrinovich Gladkov and Sofiya Borisovna Bogdanova. About the possibility of synchronization in dynamical systems. *Journal of Physics: Conference Series*, 1479(1):012011, mar 2020. doi:10.1088/1742-6596/1479/1/012011.

[37] Juan Antonio Rojas-Quintero, François Dubois, and Hedy César Ramírez-de Ávila. Riemannian formulation of pontryagin's maximum principle for the optimal control of robotic manipulators. *Mathematics*, 10(7), 2022. doi:10.3390/math10071117.

[38] Nima Shirafkan, Pierre Gosselet, Franz Bamer, Abdelbacet Oueslati, Bernd Markert, and Géry de Saxcé. Constructing the hamiltonian from the behaviour of a dynamical system by proper symplectic decomposition. In Frank Nielsen and Frédéric Barbaresco, editors, *Geometric Science of Information*, pages 439-447, Cham, 2021. Springer International Publishing. doi:10.1007/978-3-030-80209-7_48.

[39] Zdravko Terze, Viktor Pandža, Marijan Andrić, and Dario Zlatar. Lie group dynamics of reduced multibody-fluid systems. *Mathematics and Mechanics of Complex Systems*, 9:167-177, 05 2021. doi:10.2140/memocs.2021.9.167.
[12] Jean-Maris Souriau. *Structure des systèmes dynamiques*. Maîtrises de mathématiques. Dunod, Paris, France, 1970.

[13] Adi Nusser and Enzo Branchini. On the least action principle in cosmology. *Monthly Notices of the Royal Astronomical Society*, 313(3):587-595, 04 2000. doi:10.1046/j.1365-8711.2000.03261.x.

[14] Alberto Rojo and Anthony Bloch. *The Principle of Least Action: History and Physics*. Cambridge University Press, 2018.

[15] Richard Courant. Variational methods for the solution of problems of equilibrium and vibrations. *Bulletin of the American Mathematical Society*, 49(1):1-23, 1943.

[16] Gregoire Allaire and Alan Craig. *Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, UK, 2007.

[17] Nathan M. Newmark. A Method of Computation for Structural Dynamics. *Journal of the Engineering Mechanics Division*, 85(3):67-94, 1959. doi:10.1061/JMCEA3.0000098.

[18] Jesús María Sanz-Serna. Symplectic integrators for Hamiltonian problems: an overview. *Acta Numerica*, 1:243-286, 1992. doi:10.1017/S0962492900002282.

[19] C. Kane, J. E. Marsden, M. Ortiz, and M. West. Variational integrators and the newmark algorithm for conservative and dissipative mechanical systems. *International Journal for Numerical Methods in Engineering*, 49(10):1295-1325, 2000. doi:10.1002/1097-0207(20001210)49:10<1295>:AID-NME993>3.0.CO;2-W.

[20] Ernst Hairer, Gerhard Wanner, and Christian Lubich. *Geometric Numerical Integration*. Springer Series in Computational Mathematics. Springer-Verlag Berlin Heidelberg, 2006. doi:10.1007/3-540-30666-8.

[21] René de Vogelaere. Methods of integration which preserve the contact transformation property of the hamilton equations. 1956. doi:10.7274/24733422.v1.

[22] Michel Géradin and Daniel J. Rixen. *Mechanical Vibrations: Theory and Application to Structural Dynamics*. John Wiley & Sons, Chichester, UK, 2015.

[23] Anil K. Chopra. *Dynamics of Structures: Theory and Applications to Earthquake Engineering*. Prentice-Hall International Series in Civil Engineering and Engineering Mechanics. Pearson Education Limited, Harlow, UK, 5th edition, 2020.

[24] Frédéric Boyer, Vincent Lebastard, Fabien Candelier, and Federico Renda. Extended hamilton's principle applied to geometrically exact kirchhoff sliding rods. *Journal of Sound and Vibration*, 516:116511, 2022. doi:10.1016/j.jsv.2021.116511.

[25] Frédéric Boyer, Andrea Gotelli, Philipp Tempel, Vincent Lebastard, Federico Renda, and Sébastien Briot. Implicit time-integration simulation of robots with rigid bodies and cosserat rods based on a newton-euler recursive algorithm. *IEEE Transactions on Robotics*, 40:677-696, 2024. doi:10.1109/TR0.2023.3334647.

[26] François Dubois and Juan Antonio Rojas-Quintero. A variational symplectic scheme based on simpson's quadrature. In Frank Nielsen and Frédéric Barbaresco, editors, *Geometric Science of Information*, pages 22-31, Cham, 2023. Springer Nature Switzerland. doi:10.1007/978-3-031-38299-4_3.

[27] François Dubois and Juan Antonio Rojas-Quintero. Simpson's quadrature for a nonlinear variational symplectic scheme. In Emmanuel Franck, Jürgen Fuhrmann, Victor Michel-Dansac, and Laurent Navoret, editors, *Finite Volumes for Complex Applications X—Volume 2, Hyperbolic and Related Problems*, pages 83-92, Cham, 2023. Springer Nature Switzerland. doi:10.1007/978-3-031-40860-1_9.

[28] Pierre-Arnaud Raviart and Jean-Marie Thomas. *Introduction à l'analyse numérique des équations aux dérivées partielles*. Collection Mathématiques appliquées pour la maîtrise. Masson, Paris, France, 1983.

[29] C. Gear. Simultaneous numerical solution of differential-algebraic equations. *IEEE Transactions on Circuit Theory*, 18(1):89-95, 1971. doi:10.1109/TCT.1971.1083221.

[30] Monique Chyba, Ernst Hairer, and Gilles Vilmart. The role of symplectic integrators in optimal control. *Optimal Control Applications and Methods*, 30(4):367-382, 2009. doi:10.1002/oca.855.

[31] Dina Razafindralandy, Vladimir Salnikov, Aziz Hamdouni, and Ahmad Deeb. Some robust integrators for large time dynamics. *Advanced Modeling and Simulation in Engineering Sciences*, 6(1):5, 2019. doi:10.1186/s40323-019-0130-2.

[32] Ge Zhong and Jerrold E. Marsden. Lie-poisson hamilton-jacobi theory and lie-poisson integrators. *Physics Letters A*, 133(3):134-139, 1988. doi:10.1016/0375-9601(88)90773-6.

[33] Wolfram Research, Inc. *Mathematica*, Version 12.3. Champaign, IL, 2021.

[34] Nail Hairullovich Ibragimov. *Elementary Lie Group Analysis and Ordinary Differential Equations*. Wiley, Chichester, UK, 1999.

[35] Peter E. Hydon. *Symmetry Methods for Differential Equations: A Beginner's Guide*. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2000.

[36] Sergey Oktyabrinovich Gladkov and Sofiya Borisovna Bogdanova. About the possibility of synchronization in dynamical systems. *Journal of Physics: Conference Series*, 1479(1):012011, mar 2020. doi:10.1088/1742-6596/1479/1/012011.

[37] Juan Antonio Rojas-Quintero, François Dubois, and Hedy César Ramírez-de Ávila. Riemannian formulation of pontryagin's maximum principle for the optimal control of robotic manipulators. *Mathematics*, 10(7), 2022. doi:10.3390/math10071117.

[38] Nima Shirafkan, Pierre Gosselet, Franz Bamer, Abdelbacet Oueslati, Bernd Markert, and Géry de Saxcé. Constructing the hamiltonian from the behaviour of a dynamical system by proper symplectic decomposition. In Frank Nielsen and Frédéric Barbaresco, editors, *Geometric Science of Information*, pages 439-447, Cham, 2021. Springer International Publishing. doi:10.1007/978-3-030-80209-7_48.

[39] Zdravko Terze, Viktor Pandža, Marijan Andrić, and Dario Zlatar. Lie group dynamics of reduced multibody-fluid systems. *Mathematics and Mechanics of Complex Systems*, 9:167-177, 05 2021. doi:10.2140/memocs.2021.9.167.