Simpson’s Variational Integrator for Systems with Quadratic Lagrangians

: This contribution proposes a variational symplectic integrator aimed at linear systems issued from the least action principle. An internal quadratic ﬁnite-element 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.


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: 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 the fractional form) have been developed [6].Some applications involve solving initialvalue 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's evolution in the form of differential equations.It is valid across disciplines such as classical and quantum mechanics, cosmology, electromagnetism, optics, and relativity [10][11][12][13][14].As such, this variational principle is closely involved in the development of the finite-element method [15], which is used for the space and time integration of differential equations [16].
Numerical schemes for dynamical systems issued from the principle of least action are typically referred to as variational [17][18][19][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][19][20][21].One remarkable example of such methods is Newmark's integrator [17,19], which is very popular for solving problems in the dynamics of structures [22,23] and has recently been geometrized to solve the 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-element 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.The 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-degree 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-degree 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's convergence is revealed to be of the fourth order.The manuscript ends with a brief discussion and concluding remarks in Section 7.

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 where L is the system Lagrangian.We focus on dynamical systems for which the Lagrangian can be expressed quadratically as where M and K are symmetric, positive-definite n-dimensional matrices with constant coefficients; q ∈ R n .
We can discretize Equation ( 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 Equation (2): where L d (q , q r ) is the discrete form of the Lagrangian (3).Subscripts and r stand for "left" and "right" values, respectively.Let us consider a centered finite-difference approximation: dq dt q r − q h , and a midpoint quadrature: The discrete Lagrangian becomes

Discrete Euler-Lagrange Equations
The discrete action (4) being a sum, only two terms contain the variables q j : So, when the discrete action is stationary (δS d = 0 for arbitrary variations δq j of the states q j ), only two terms remain.Necessarily, ∂L d ∂q r q j−1 , q j + ∂L d ∂q q j , q j+1 = 0. ( The generalized momenta p j ∈ R n are defined, on the right, as Therefore, the first term of Equation ( 5) is identified as p j , so applying Equation (6) in Equation ( 5) leads to Then, p j+1 is constructed following Equation (6): Using Equations ( 7) and ( 8), it can be established that Equations ( 9) are consistent with dp dt = −Kq and p = M dq dt , respectively.

Newmark's Scheme
System (9) can then be arranged in linear form as where η ∈ R 2n , η = (p, q) T and I n is the n-dimensional identity matrix.Newmark's symplectic scheme is obtained by matrix inversion of Equation (10).We can establish that where matrices A and B are defined in Equation ( 11) above.It has been observed that this particular variant of Newmark's method is unconditionally stable and second-order convergent [19].

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 Equation ( 2) and the Lagrangian has the structure of Equation (3).
Note that q(0) = q , q h 2 = q m and q(h) = 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, where derivatives g , g r ∈ R n are given by Gear's scheme [29].Gear's scheme is used as the differentiation approximation for q(t) ∈ P 2 as where g m ∈ R n .The above confirms that a first-order centered finite difference is recovered by 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 q j q(t j ), ∀ 0 j N; Taking Equation ( 14), q(t) is a quadratic polynomial vector function within the interval [t j , t j+1 ] with t = t j + θh, q = q j , q m = q j+ 1 /2 , q r = q j+1 .

Discrete Lagrangian
Let us recall that the continuous action is defined by Equation ( 2) and that the Lagrangian is defined by Equation (3).In the present case, the discrete action sum Σ d for a motion t → q(t) between the initial time and a given final time T > 0 is discretized with N regular intervals as where L h (q , q m , q r ) is the discrete form L h (q , q m , q r ) h 0 L dt, of the Lagrangian (3).Using Simpson's rule (1), the polynomial approximation (14) of the states, and derivatives (15), the discrete Lagrangian of a linear system is expressed as

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 q j and q j+ 1 /2 : Maupertuis's stationary action principle [10] implies that δΣ d = 0 for an arbitrary variation of the internal degree of freedom δq j+ 1 /2 ∈ [t j , t j+1 ].Considering Gear's scheme (15), we have where g i is the i-th component of g and q k is the k-th component of q.
When δΣ d = 0, ∂L h ∂q m q j , q j+ 1 /2 , q j+1 = 0 by necessity.This conforms to the discrete Euler-Lagrange equations at the middle of the interpolation interval: However, g j − g j+1 = 1 h −4q j + 8q j+ 1 /2 − 4q j+1 , so Equation ( 17) becomes This last equation is consistent with M dt 2 + Kq = 0. Additionally, for an arbitrary variation δq j , the Euler-Lagrange equations are given by the necessary condition that The generalized momenta p j are defined, on the right, as Therefore, the first term of Equation ( 19) is identified as p j , and it can established that because −3g j − 4g j+ 1 /2 + g j+1 = 1 h 14q j − 16q j+ 1 /2 + 2q j+1 .Equation ( 18) is then multi- plied by h /3, and the result is added to Equation (21).This eliminates q j+ 1 /2 from the first term of the right-hand side: Then, p j+1 is calculated according to Equation ( 20) 18) is then multiplied by − h /3, and the result is added to Equation (23).This eliminates q j+ 1 /2 from the first term of the right-hand side: Using Equations ( 22) and ( 24), we can establish that Equations ( 25) are consistent with dp dt = −Kq and p = M dq dt , respectively.Note that the term h 2 12 K in the second equation above vanishes as h → 0.

First Variant of Simpson's Scheme
The system composed of Equations ( 18) and ( 25) can be rearranged as where System (26) can then be arranged in linear form: where η ∈ R 2n , η = (p, q) T ; The first variant of Simpson's scheme is obtained by matrix inversion of Equation ( 27).We can establish that q j+ 1 /2 where matrices A σ and B σ are defined in Equation (28) above.

Second Variant of Simpson's Scheme
Simpson's scheme's internal degree of freedom can be eliminated using the first equation of System (26): This equation approximates the middle point when h → 0, because then L → I n .Substituting this value into the third equation of System (26) leads to and the second equation of System (26) remains unchanged.Therefore, the internal degree of freedom is successfully eliminated so that, now, where η ∈ R 2n ; η = (p, q) T ; The second variant of Simpson's symplectic scheme is obtained by matrix inversion of Equation (30).We can establish that where matrices A s and B s are defined in Equation ( 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 Equation (11) and Equation (31).The symplectic property of both schemes is analyzed in 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.

Symplectic Property
A symplecticity proof is obtained by verifying that Φ corresponds to the scheme transformation matrix and characterizes a discrete time evolution of the system.J is sometimes referred to as the canonical matrix for Hamiltonian systems [30] and has the property that J −1 = J T = −J.When Equation (33) holds, it means that Φ is an area-preserving transformation and that the scheme (12) is symplectic (see, e.g., [18][19][20]31] for more details on this demonstration).Proposition 1.An implicit scheme of the type are square, partitioned, invertible matrices and blocks X and Y are symmetric and positive-definite.
Proof of Proposition 1.Let us first make explicit the transformation A −1 .Since A is square and partitioned, its inversion is performed using auxiliary variables α and β.Let us establish that Subtracting both equations above gives Since Z is the sum of two symmetric, positive-definite matrices, it is invertible.Equation ( 36) is then substituted into the first equation of System (35): Matrix A from Equation ( 34) is inverted in Equations ( 36) and (37) as: Then, it suffices to verify Equation (33) with Φ = A −1 B. Thus,
Proof of Proposition 2. In Newmark's scheme's formulation (12), matrices A n and B n from Equation ( 11) are of the form of Equation ( 34), because X n and Y n (11) are symmetric and positive-definite.By Proposition 1, Newmark's scheme is symplectic.

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 Equation (34).For this, blocks X s and Y s are required to be symmetric and positive-definite.
Proof of Proposition 3. Since M and K are symmetric, Y s is symmetric if and only if its first term: is symmetric, Y s is also symmetric.
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: where (µ, κ) are the smallest and (m, k) are the largest eigenvalues of matrices M and K, respectively; ϕ = 0 is an eigenvector.Taking M 1 /2 ϕ = ψ and then and so, The above expression is positive for This inequality is a sufficient stability condition for Simpson's scheme.Let us remark that k /µ 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 The stability condition (40) can also be stated as This condition is similar to the stability condition characterizing the mono-dimensional case for Simpson's scheme [26].
Proof of Proposition 4. When the condition (40) is met, Equation (39) becomes and LK −1 is positive-definite.Therefore, KL −1 is also positive-definite, and recalling Proposition 3, it is symmetric.Consequently, and Y s is positive-definite.Now, only the positive-definiteness of block X s from Equation (31) remains to be proven.
Proof of Proposition 5. From inequality (38), and substituting the condition (40) for the first term of the right-hand side of the above inequality, Therefore, matrix X s is positive-definite.
Proof of Proposition 6.For the second variant of Simpson's scheme (32), matrices A s and B s from Equation ( 31) are of the form of Equation ( 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.
These results prove that the proposed Simpson's scheme is symplectic.

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 φ(p, q) verifies φ(p j+1 , q j+1 ) = φ(p j , q j ), where (p j+1 , q j+1 ) and (p j , q j ) satisfy the dynamics of Simpson's scheme, Equations ( 31) and (32).Proposition 7. Given an implicit scheme of the type where are square, partitioned, invertible matrices and blocks X and Y are symmetric and positive-definite, there exists a quadratic form: Proof of Proposition 7. Let us expand Equation (41): The above can also be written as Therefore, by multiplying ξ by the first equation above on the left and by the second equation above on the right, p j+1 + p j T ξ p j+1 − p j = X q j+1 − q j T ξ(−Y) q j+1 + q j = − q j+1 − q j T XξY q j+1 + q j .( 43) is symmetric and positive-definite, it is deduced that ζ = XξY is symmetric and positivedefinite.Following from Equation (43), p j+1 T ξp j+1 + q j+1 T ζq j+1 = p j T ξp j + q j T ζq j , and the property is proven since φ(p j+1 , q j+1 ) = φ(p j , q j ).

Linear Double Pendulum Model and Exact Solution
This section presents a case study for subsequent numerical experiments.

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.The 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 where an overdot indicates time differentiation.Potential energy is calculated as and finally, the system Lagrangian L = T − V can be explicated as 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) 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 q = (q 1 , q 2 ).
Using Equation (45), the linear form L L of Lagrangian (44) becomes where the second term of the cos(q 1 − q 2 ) approximation in Equation ( 45) vanishes when multiplying the product q1 q2 .Generalized momenta are defined as According to the Lagrangian (46), we have ).
Motion equations are then obtained by applying Euler-Lagrange equations

Exact Solution
Equation (47) can also be established as a linear system of the form where The general solution of Equation ( 48) is of the form where x 1 and x 2 are eigenvectors and ω denotes the oscillation frequency.Two characteristic frequencies (ω 1 , ω 2 ) are determined by the solution of the auxiliary equation det Let us focus on the case where In this particular case, the oscillation frequencies are given by with a mass ratio µ r = m 2/m 1 and frequency ω 0 = √ g /l.
Eigenvectors x 1 and x 2 are then obtained by solving (K − ω i 2 M)x i = 0 for i = 1 and i = 2: Solving the above system gives Finally, using Equation (50), the general solution of Equation (48) (or Equation ( 47)) can be established as where constants (c 1 , c 2 , ϕ 1 , ϕ 2 ) are given by the chosen initial conditions on the positions and velocities.

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.The results are provided for reference since RK4 is among the most popular methods available.For benchmarking purposes, we applied 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.The computations were carried out using Wolfram's Mathematica software (version 12.3) [33].The figure plots were then created using exported data with the pgfplots package from L A T E X. Table 1 specifies the constants and initial conditions used for all of our simulations.Using these values and following Equation (51) with null initial phases (ϕ 1 , ϕ 2 ), the exact solution that serves as the main reference in our comparisons is where ω 1 and ω 2 are given by Equation (49).
Table 1.Constants and initial conditions used for numerical simulations.

Constants Initial Conditions
Frequency ω 0 is used to show the results in terms of an oscillation period t such that t = 1 ω 0 .
Therefore, both the total simulation duration T and step size h are given in terms of 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.

Configuration Parameters and Generalized Momenta
We begin by comparing the configuration parameter solutions 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 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.
Configuration parameters', q, evolution for the linear double pendulum.Initial conditions are specified in Table 1.The step size is fixed as h = 0.1 t. Simpson's integrator tracks the exact solution with more precision than Newmark's method and Runge-Kutta's integrator.
Figures 2 and 3 show that Simpson's integrator is more precise than both Newmark's and Runge-Kutta's integrators on a short simulation (T = 1 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 4.
Newmark's integrator precision can be increased by refining the step size.With h = 0.01 t, the 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.
Configuration parameters', q, 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 t. Simpson's solutions correctly follow the exact solution for longer simulations.

Phase Portraits
With a step size of h = 0.1 t, Newmark's solutions' deviations are particularly visible when tracing the motion phase portrait.Figure 5 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.The total simulation time was limited to T = 3 t for visualization purposes.

Exact Newmark Simpson
Figure 5. Phase portraits' evolution for the linear double pendulum.Initial conditions are specified in Table 1.The step size is fixed as h = 0.1 t, and three periods are shown (T = 3 t).Simpson's phase portrait clearly follows the exact one.

Energy Conservation
The following function gives the system energy: 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 System (25) introduces the small and vanishing quantity − h 2 12 K q j+1 −q j h 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 6 shows that Simpson's solutions lead to a non-conserved energy H(p, q).Nevertheless, the maximum relative error with respect to the initial value is extremely small even for h = 0.1 t, as the 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 6).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.The 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.
Proposition 7 shows that Simpson's scheme preserves a quadratic form given by the function of Equation (42).Matrices ξ s and ζ s (where subscript s stands for Simpson) are according to Proposition 7 as Figure 7 plots the absolute error on function φ(p, q) of Equation (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.

Convergence
The error e(t) = q(t) − q ex (t), and its convergence rate is measured following the prescriptions found in [16].The schemes' precision was evaluated using an ∞ error norm e ∞ = q − q ex ∞ = sup n |q n − q ex n |.Several simulations were performed for decreasing values of h between h = 0.1 t and h = 0.001 t.The e ∞ norm was calculated for each case.These errors are plotted in Figure 8, 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 8).These trials confirmed 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 performed 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].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 t.It can be observed that Newmark's scheme's 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 ∞ 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 2. Convergence order with respect to simulation length for motion simulations performed for a linearized double pendulum (see Figure 1).The initial conditions are specified in Table 1.Error norms e ∞ for momenta p and states q increase with simulation length T. Newmark's scheme's convergence decays 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.

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.Simpson's numerical scheme presented 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 was given, and it was also proven 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 fourthorder accurate on both the states and generalized momenta.Numerical evaluations revealed 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 non-linear 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 non-linear 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 in simulating complex non-linear 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]. 0 Figure 6.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.The 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.

Figure 7 .
Figure 7. Simpson's scheme preserves the quadratic form φ(p, q) of Equation (42).The absolute errors are minimal and may come from accumulated rounding errors.
Generalized momenta's, p, evolution for the linear double pendulum during ten periods.Initial conditions are specified in Table1.The step size is fixed as h = 0.1 t. Simpson's solutions correctly follow the exact solution for longer simulations.