{'fraction_non_alphanumeric': 0.04591579301452038, 'fraction_numerical': 0.01813645792453888, 'mean_word_length': 4.130140926085706, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Symplectic P-stable Additive Runge--Kutta Methods

Symplectic partitioned Runge--Kutta methods can be obtained from a variational formulation where all the terms in the discrete Lagrangian are treated with the same quadrature formula. We construct a family of symplectic methods allowing the use of different quadrature formulas (primary and secondary) for different terms of the Lagrangian. In particular, we study a family of methods using Lobatto quadrature (with corresponding Lobatto IIIA-B symplectic pair) as a primary method and Gauss--Legendre quadrature as a secondary method. The methods have the same favourable implicitness as the underlying Lobatto IIIA-B pair, and, in addition, they are \emph{P-stable}, therefore suitable for application to highly oscillatory problems.

1. Introduction. In this paper we introduce a family of Runge-Kutta methods of additive type particularly suited to highly oscillatory problems. Our method are derived from a variational formulation, using different quadrature formulas for different parts of the Lagrangian. We will consider mainly a formulation where we use a primary method (giving rise to a symplectic PRK) and a secondary method, that is based on different quadrature weight and nodes. The final formulation of the hybrid method can be classified as special subclass of symplectic additive Runge-Kutta (ARK) methods. ARK were introduced already in the 80's [2] to deal with stiff ODEs. In addition, we mention the more recent GARK methods [11], which were also proposed to add flexibility in treating different force terms in the differential equation by different sets of coefficients. By opportune repetitions of coefficients or introduction of zero patterns, ARK can be put into a GARK form and vice versa. In the context of algebraic differential equations, similar approaches have been followed by [5] with special attention to structure preservation in Hamiltonian systems. Recently, GARK methods for stiff ODEs and DAEs were considered in [13] with focus on the combination Gauss/Radau IIA and Gauss/Lobatto IIIC.
Our motivation comes from the study of highly oscillatory problems, trigonometric integrators (see [3] and references therein), in particular, the intriguing properties of the second-order implicit-explicit (IMEX) method originally proposed by [18] and further analyzed in a variational setting in [12] and as a modified trigonometric integrator in [8]. This method is equivalent to applying the "midpoint rule" 1 to the fast, linear part of the system, and the leapfrog (Störmer/Verlet) method to the slow, nonlinear part. It has the following properties: (i) it is symplectic; (ii) it is free of artificial resonances; (iii) it is the unique method that correctly captures slow energy exchange to leading order; (iv) it conserves the total energy and a modified oscillatory energy up to to second order; (v) it is uniformly second-order accurate in the slow components; and (vi) it has the correct magnitude of deviations of the fast oscillatory energy, which is an adiabatic invariant [8].
The Störmer/Verlet method belongs to the family of Lobatto IIIA-B partitioned Runge-Kutta methods (PRK). In an unpublished report from 1995, Jay and Petzold studied the linear stability of Lobatto PRK and proved that none of the methods in this family is P-stable, as they are not unconditionally stable when applied to the harmonic oscillator. They concluded that these methods were not suitable for highly oscillatory systems [6]. Further stability properties of these Lobatto PRK were also studied in the context of multisymplectic integration and the wave equation in [9].
Being the lack of P-stability well established for Lobatto PRK, it is therefore quite a surprise that the combination implicit midpoint and Störmer/Verlet is unconditionally stable. Intrigued by the properties of the IMEX, [17] introduced a family of symplectic, unconditionally stable modified trigonometric integrators of second order, which included the IMEX as a special case.
In this paper we construct higher order integrators pursuing the variational approach of the Lagrangian formalism. The technique used is very close to the one described in [3] for the derivation of symplectic PRK methods. The main idea is similar to that described above for the second order IMEX: to use a Lobatto method for the kinetic energy and slow potential (the latter being costly to compute) and Gauss-Legendre of the same order for the linear highly oscillatory part (easy to compute). To avoid the introduction of further function evaluation of the potential, we approximate the internal stages values by two techniques, interpolation and collocation. Although we focus especially on the Lobatto and Gauss-Legendre combination as primary and secondary method respectively, the derivation presented is general and applies to different combinations of primary and secondary methods.
Other variational approaches exist, especially in the community of computational mechanics, see for instance [7]. Recently, the latter approach has been used, together to a splitting of the Lagrangian, in the context of higher order variational integrators for dynamical systems with holonomic constraints [15] and in order to devise mixed order integrators for systems with multiple scales [14]. The approach in [15] and the one presented in this paper have several common features but also diversities, like for instance the way the Lagrange multipliers are handled. Having said this, it is not unlikely that some of the methods derived by the two approaches will coincide for some similar choices of coefficients and some problems, but a thorough comparison is outside the scope of the present paper.
The paper is organized as follows. In Section 2 we show the general theory for the derivation of the methods and how to construct the coefficients by either interpolation or collocation. In Section 3 we prove some results on the order of the proposed methods. In Section 4 we study the P-stability of the methods and in Section 5 we show how the methods can be put in the framework of modified trigonometric integrators. In Section 6 we show several numerical tests on the Fermi-Pasta-Ulam-Tsingou problem and compare with higher order construction of the IMEX method using the Yoshida time-stepping technique. Finally, we have some concluding remarks and in the Appendix we present explicitly tables with the coefficients for the methods of the Lobatto-Gauss-Legendre family based on interpolation and collocation for order two, four and six.
2. Variational derivation. It is well known that symplectic Partitioned Runge-Kutta methods (PRK) can be obtained by a variational approach, doing discrete variations on a discrete Lagrangian approximating the continuous Lagrangian L(q,q) [3].
Consider a Lagrangian L(q,q) and assume that it can be written as sum of two (or more) terms, Whereas the derivation of symplectic PRK uses the same quadrature for all the terms, we consider the case when one would like to use a different quadrature for one or more terms in the sum. A motivating example is the case of highly oscillatory problems in molecular dynamics, with a Lagrangian of the form where V 1 is a slow potential while V 2 is a fast oscillating potential, for instance of the form V 2 = − 1 2 q T Ω 2 q, Ω being a diagonal matrix with elements ω i 1. A natural splitting in this context would be In this paper, we restrict the discussion to the case when the Lagrangian is split in two terms as above, but the generalization to several terms is straightforward. We focus on a discrete Lagrangian of the form with where the coefficients (A, b, c) are the coefficients of a standard RK method with s 1 stages (primary method), while (b,c) are the weights and nodes of the secondary quadrature with s 2 weights and nodes respectively. To avoid the introduction of extra internal stages due to the secondary method, we assume that theQ k can be written asQ for some coefficientsã k,j to be determined, k = 1, . . . , s 2 and j = 1, . . . , s 1 . We consider the discrete action A = N −1 n=0 L h (q n , q n+1 ), and require the discrete action to be stationary, assuming no variation at the endpoints q 0 and q N , leading to the usual discrete Euler-Lagrange equations, which, by introduction of the conjugate variables, lead to a one-step method, (q n , p n ) → (q n+1 , p n+1 ). Because of (2) and (4), the action also introduces additional variables, theQ i s, subject to the constraint (3). We thus require the action A to be stationary also with respect to these variables, and handle the constraint (3) by the method of Lagrangian multipliers by augmenting the discrete Lagrangian (1) as In a manner very similar to the derivation of symplectic PRK described in [3], derivation of (6) with respect to λ imposes the constraint (3), while derivation with respect to theQ j gives the relation between the multiplier λ and the other variables, We set With this notation, and using the relations ∂Qi ∂Qj = ha ij I, ∂Qi ∂Qj = hã ij I, the constraint conditions in equation (7) read The final formulation of the method is obtained from the conjugate variables p n and p n+1 as in (5) and thereafter eliminating λ using (10). By direct computation, we have where in (12) we have used (10) and in (13) we have used s1 l=1 b l ∂Q l ∂qn = −I which comes from the derivation of (3).
By a similar token, we find so that, eliminating λ, we get the relation To obtain the definition of the P j , we use (13) and substitute in (10) to obtain We recognize that the first set of coefficients is a i,j = b j − b j a j,i /b i , so that the L 1 part is treated with a classical symplectic pair of PRK [3]. The second set of coefficients imposes the symplecticity condition for the use of the secondary method in the treatment of the L 2 .
2.1. General format of the methods. The generalization to L = L 1 + L 2 + · · · is straightforward and leads to a symplectic subclass of ARK methods. In this paper we describe in detail the case with two terms, that is L 1 = 1 2q Tq − V 1 (q), L 2 = −V 2 (q). Let −∇V 1 = F 1 and −∇V 2 = F 2 be the forces corresponding to the potentials V 1 , V 2 . Denote by (A, b, c) the primary method so that, with ( A, b, c), it forms symplectic PRK pair. Let (b,c) be the secondary method (only quadrature weights and nodes are necessary). We have The first two equations are in terms of the Q i andQ k variables, the third equation defines the P i s, and the last equation means that theP k variables play an auxiliary role only, see (9). The constraint relation (7) kãk,j F 2 (Q k ) = λ j b j allows us to find the derivativesQ i (= P i ) at the intermediate stages. The full method reads where ã i,k =b k −b kãk,i bi . In matrix notation, where 1 m1,m2 is the m 1 × m 2 matrix with all entries equal to one.
The method requires the evaluation of extra internal stages for the secondary method (theQ i ) but only as many function evaluations of F 1 as for the underlying (symplectic) PRK method, allowing for different number of function evaluations for F 2 . This can be particularly interesting when F 1 is expensive, while F 2 is cheap to compute.
The methods (17) can be applied to an Hamiltonian systeṁ with Hamiltonian energy H(q, p) = p Tq − L(q,q) where L(q,q) = 1 2q Theorem 2.1. The methods (17) with Ã and A as in (18) and (19) are symplectic.
Proof. Follows immediately from the fact that the variational derivation of the methods uses essentially generating forms of the first kind, that is, L h (q n , q n+1 ) is a type-I symplectic generating function for the method.
While the symplectic requirements for the matrix A are well known in the the context of PRK, those for Ã are the same as those derived by algebraic arguments in [11,5].
One of the most common approaches to derive coefficients for the coupling matrix A in methods of the type (17), is to impose order condition in addition to the linear constraints (18) and (19) which guarantee symplecticity. The solution of the overall nonlinear system needs not be unique, as we shall see below. In this paper, we follow a different approach which leads to at least two solutions for the matrixÃ and consequently Ã .
For simplicity, we restrict our analysis to the interval [t 0 , t 0 + h].

2.2.
Construction of the matrixÃ. There are two natural choices to construct the approximationsQ j in (4). We assume that the primary method is described by the RK tableau .
For the secondary method, we assumec k andb k , the nodes and weights of the secondary method, to be given.
Interpolation: Given the primary nodes be the ith cardinal Lagrange polynomial and construct the interpolating poly-nomialQ where the Q l = q 0 + h s1 j=1 a l,jQj are obtained by the primary method. Substituting Q l in (20), recalling that j L j (τ ) = 1 and computing in the nodes t 0 +c i h of the secondary method, we recover (4), with coefficients Let L(c) the s 2 × s 1 matrix with elements L(c) i,j = L j (c i ) of the primary Lagrange cardinal polynomials evaluated in the secondary nodes. Theñ where A is the coefficient matrix of the primary method. Collocation: Another natural choice is to use the interpolation of theQ j s: we use the cardinal interpolating polynomials L j (τ ) constructed with the nodes of the primary method to constructQ ≈ s1 j=1 L j Thereafter, evaluating in the secondary nodes t 0 +c i h, we recover (4) with coefficients 3. Order of the methods. A general treatment of the order conditions for these ARK methods can be developed using the algebraic tree theory in a manner very similar to the order analysis of the ARK, GARK methods [11,13] using the formalism of colored trees [3]. The order conditions are used to derive the coefficients of the methods. In our setting, the primary method, leading to a PRK pair (A, A, b, c), and the secondary method, (Ã, Ã ,b,c), are given by the choices (21)-(22), but the order of the resulting method (17) is not obvious.
Lemma 3.1. Assume that for the primary method, Ac k−1 = 1 k c k , where the power is intended component-wise on the vector elements. With the same notation as above, if s 1 ≥ s 2 ≥ 1, we havẽ In particular, s1 j=1ã i,j =c i , i = 1, . . . , s 2 . Moreover, if: 1. the quadrature formula based on the nodesb i is exact for polynomials of degree at least s 1 − 1 for the interpolation (21) and the primary method satisfies for all j; and 2. the quadrature formula based onb i and b i are of order at least s 1 + 1 for the collocation (22), then we have Proof. We first prove (23) in the case k = 1 (c 0 = 1 s1 ).
For the interpolative scheme (20), we have where the second last passage holds provided that j a l,j = c l , which is true as long as the primary method has order at least one. The function  The proof for k > 1 for the interpolative methods follows by a similar argument as for k = 1, using the property of the primary RK method that Ac k−1 = c k k and the fact that the L i are interpolating polynomials on s 1 nodes interpolating exactly up to degree s 1 − 1.
For the collocative stages (22), we have since the L j interpolate exactly polynomials up do degree s 1 − 1 as above.
For the proof of (24), we observe that where in the last passage we have used the fact that b i = 1. In the interpolative setting (21), we look at the term In the collocative setting (22), since the f i s are polynomials of degree s 1 and the quadrature formula based on the nodesb i has order at least s 1 .
The second last passage follows provided that the integration formula with weights and nodes (b, c) is exact for polynomials of degree completes the proof.
Proof. To prove the theorem it is sufficient to show that that quadrature formula interpolating the nodesc using the nodes c, as well as the quadrature formula collocating the nodesc, have also order r when f (x) = x k , k = 0, . . . , r, for which 1 0 x k dx = 1 k+1 . We start with proving (27) for the interpolative formulas. For k = 0, . . . , s 1 − 1 the statement is immediate as the function j L j (x)f (c j ) exactly interpolates polynomials of degree up to degree k = s 1 − 1, hence j L j (c i )c k j =c k i . Hence by virtue of (26) the statement follows.
For the collocative formulas and (28), when f (x) = x k , we have f (x) = kx k−1 so that the interpolation j L j (x)c k−1 j = x k−1 is exact for polynomials of degree k = 0, . . . , s 1 . Consequently,f i =c k i and the statement follows. When k = s 1 +1, . . . , r− 1, we refer again to the computations in Lemma 3.1: which completes the proof.
We are especially interested on the family of methods generated by the Lobatto IIIA-B (primary method) and Gauss-Legendre (secondary method) of the same order (r = 2s 2 = 2(s 1 − 1)). These quadrature formulas are superconvergent and the proof of superconvergence is heavily based on the roots and weights of the corresponding orthogonal polynomials, so that, in principle, the interpolation might destroy the superconvergence. Fortunately, this does not happen because the methods satisfy the hypotheses of Theorem 3.2, and the order is preserved. This statement is summarized in the Corollary below. (17) with primary method Lobatto IIIA-B with s 1 stages and secondary method Gauss-Legendre with s 2 = s 1 − 1 stages has order r = 2(s 1 − 1) = 2s 2 both for coefficients based on interpolation and collocation. 4. P-stability. P-stability is a desirable property when applying a numerical method to highly oscillatory systems. The test model is the harmonic oscillator

Corollary 1. The methods
whose exact solution can be written as where D ω = diag(1, ω). It is well known that the application of a s-stages PRK pair with coefficient (A, b) and ( A, b) yields a numerical approximation We are interested in methods that preserve the unit modulus of the eigenvalues of the rotation matrix Θ(µ). It is well known that symmetric RK methods are P-stable, and, as a consequence, the methods Lobatto IIIA and Lobatto IIIB, taken individually, are P-stable. However, the PRK combination Lobatto IIIA-B, which include the Verlet scheme for order 2, is not P-stable [6,9]. Motivated by the positive results of the IMEX method, that was proven to be P-stable (unconditionally stable, [8]), we study the methods (17) and, the same spirit of the IMEX methods, the oscillatory part is treated by the secondary method (i.e. we set F 1 = 0).  (17) is given as Moreover, as the methods are symplectic, Proof. For the test equation (29) (F 1 = 0), the method (17) in [t 0 , t 0 + h] can be written as To ease notation, we denote by capital lettersQ the vector of the internal stages Q i , by P the vector of the internal momenta P i , and abuse notation, to avoid the use of tensor products. So, for instance,ÃP has, as an i-component, the vector a i,1 P 1 + · · · +ã i,s1 P s1 . In block form, we have which we use to solve for theQ and P . From q 1 = q 0 +hb T P and p 1 = p 0 −hω 2bTQ , we get where the last passage follows in a manner very similar as corresponding proof for PRK methods with M (µ) as in (33). As for (34), if the method is symplectic, then it must be volume preserving for Hamiltonian systems, which, in this case implies that det M (µ) = 1.
Since the eigenvalues of the matrix M (µ) are because of the determinant condition (34) one has λ 1 λ 2 = 1. Hence the eigenvalues lie on the unit circle if and only if In addition, when trM = 2, the eigenvalues are both equal to 1, while for trM = −2, the eigenvalues are both equal to −1. When studying P-stability, we will refer to the function 1 2 |trM (µ)| as stability function of the method. Provided that the method is P-stable, it can be interpreted as an oscillator with a modified frequency. Comparing with the matrix Θ(µ) in (30), we have 1 2 trM (µ) = cos(ωh) = cos(μ),μ =ωh (37) corresponding to a modified frequencyω satisfying Corollary 2. The IMEX method is P-stable.
Because of the symplecticity of the methods (17), it is obvious that in order to study the P-stability it is sufficient to look at the diagonal elements M 1,1 and M 2,2 of the matrix M (µ) in (33). By direct computation, one has that then M 1,1 = M 2,2 .
Proof. We use the formal series (I + G) −1 = k (−1) k G k . The first part of the statement says that we can push Ã andÃ on the other right hand side using (23) and (24) from Lemma 3.1, that isÃ1 s1 =c and Ã 1 s2 = c.
For the second part of the statement, if all the infinite terms of the series in M 1,1 and M 2,2 are equal for k = 1, 2, . . ., then the series are also equal, even if the series do not converge. To prove this, note that the matrices (Ã Ã ) and ( ÃÃ ) have the same n nonzero eigenvalues λ 1 , . . . , λ n , n ≤ min{s 1 , s 2 }. By the Cayley-Hamilton theorem, G m can obtained as a linear combination of I, . . . , G n−1 for m ≥ n. Therefore only the terms in (43) need be checked.
Remark. Note that for k = 0, we have b T c = 1 2 =b Tc is always verified for methods of order at least one.
The combination Lobatto IIIA-B and Gauss-Legendre of the same order satisfies the requirements of Lemma 4.3, hence it is sufficient to check (43) only up to k = 1 (method of order 4) and k = 2 for the method of order six.
where the first and third passages follow from the exactness of the quadrature, the second passage by integration by parts. Proof. Because of Lemma 4.3, it is sufficient to prove (43). Although in general Ãc k−1 = c k k for k > 1 (see (24)), it is possible to show that provided D(k) holds. To see this, note that, through the definition of Ã , The first term in the second passage follows because of Theorem 3.2, while the second term is a consequence of Lemma 4.4 and D(k). Using Lemma 3.1, it is easy to see that (44) follows.
Next, consider the termb T (Ã Ã ) kc . We have, We continue in a similar way by applying (44), so that b T (Ã Ã ) kc = 1 (2k + 1)!b Tc2k+1 = 1 (2k + 2)! by virtue of Theorem 3.2. A similar procedure is applied to the term b T ( ÃÃ ) k c: where we have used b T Ã =b T (I − diag(c)) and Lemma 3.1 onÃc. Treating recursively the term (Ã Ã ) k−1c2 as before, one arrives to the desired result.
Because of the superconvergence of the Lobatto and Gauss-Legendre quadrature, we have that all the assumptions of the previous theorem are satisfied for combinations of the method of the same order r, therefore, we can state the following result.

5.
The methods as modified trigonometric integrators. We consider the application to the test equation Theorem 5.1 (Modified trigonometric integrator). Consider the symplectic methods (20) applied to the test oscillatory problem (47). Assume that the primary method has symmetric stages and that | 1 2 trM (µ)| ≤ 1, with matrix M (µ) as in (33) having two independent eigenvectors in case of equality. Then the method can be considered as a symplectic modified trigonometric integrator with modified frequency satisfying the implicit relation and can be written in the form for s 1 implicitly defined filter functions where A i is the ith column of A. The p-variables are reconstructed from the formula where the ψ i are the same as in (50).
Proof. As in the proof of P-stability, we ease notation and denote by capital letters Q the vector of the internal stagesQ i , by P the vector of the internal momenta P i , by F (Q) the vector of the f (Q i ) and abuse notation, to avoid the use of tensor products. Thus, the expression b T F (Q) means Similarly, for matrix products, the expression AF (Q) has, as the i-component, the vector a i,1 f (Q 1 ) + · · · + a i,s1 f (Q s1 ), etc.
Proceeding as for P-stability, we see that which we use to solve for theQ and P . From q 1 = q 0 + hb T P and . One has Thus q 1 is given by where, as above, µ = ωh andμ =ωh is the modified frequency and F (Q + ) indicates that the internal stages are in [0, h]. The cos(μ) and sinc(μ) terms come form D ω M (µ)D ω −1 in the usual way, provided that | 1 2 M (µ)| ≤ 1. By replacing h with −h, we have where, as above, F (Q − ) indicates that indicates that the internal stages are in [0, −h] Taking the sum of q 1 and q −1 , we obtain while subtracting the two expressions, we obtain With some simple algebraic manipulations, it is easy to recover the filter functions. The theorem statement follows by assuming that the primary method has symmetric stages.
Remark. The above theorem is also valid for all the methods described in the paper in the region where the step size h is such that | 1 2 trM | ≤ 1.
When the first node c 1 = 0 then Q 1 = Q −1 = q 0 so the first term on the right hand side of (49) becomes 2ψ 1 (μ)f (q 0 ) while it cancels in (51). Moreover, in the case of the Lobatto primary method, c s1 = 1 hence Q s1 = q 1 and Q −s1 = q −1 . However, the last column of the matrix A is zero, and so is the last filter function ψ s1 .
For the IMEX method, we have c 1 = 0, c 2 = 1 (s 1 = 2), hence (51) gives We have ψ 2 = 0 and and we recover its expression as a modified trigonometric integrator with filter functions φ = 1, ψ(ξ) = cos ξ satisfying the implicit relation cos(μ) = (1 + µ 2 4 ) −1 , as derived in [8]. Similarly, for the order four Lobatto-Gauss-Legendre method, we have ) and with filter functions ψ i , i = 1, 2, 3, satisfying the implicit relations ψ 1 (μ) = 2(−µ 2 + 12) µ 4 + 12µ 2 + 144 , ψ 2 (μ) = 2(µ 2 + 24) µ 4 + 12µ 2 + 144 , ψ 3 = 0 (53) (φ i = 1, i = 1, 2, 3). The modified frequency is given by (45). Finally, for the of order 6 Lobatto-Gauss-Legendre method, we have similar expressions, with filters implicitly defined by and modified frequency given by (46). 6. Numerical experiments. As a test bed, we consider the Fermi-Pasta-Ulam-Tsingou (FPUT, formerly FPU) problem of alternating soft and stiff springs, that has been extensively used in literature to study methods for oscillatory problems. Because of the oscillatory nature of the problem, among all the methods proposed, we test only those that are P-stable, as methods that are not P-stable are likely to produce diverging solution as soon as the step size leaves the region of P-stability. Therefore, in what follows, all the numerical experiments are performed with the Lobatto-Gauss-Legendre family (17) with coefficients by interpolation (21). We will compare these methods also with higher order integrators obtained using the IMEX (coinciding with the Lobatto-Gauss-Legendre method by interpolation of order 2) and the Yoshida time stepping technique.
6.1. The Fermi-Pasta-Ulam-Tsingou problem. For comparison with [3,8], we consider the same setup with 2 points of unit mass representing alternating soft nonlinear springs and stiff linear springs. Setting q to be the concatenation of slow (index s) and fast (index f ) position variables, q = [q s,1 , . . . , q s, , q f,1 , . . . , q f, ] T , and p the corresponding momenta, the Hamiltonian reads In our setup, the nonlinear potential and kinetic energy are treated with the Lobatto IIIA-B pair, while the linear stiff energy ω 2 2 i=1 q 2 f,i is treated with the Gauss-Legendre methods based on interpolation.
The total oscillatory energy I,  (3) shows the oscillatory energies for each spring and the total oscillatory energy, comparing the IMEX method (which is the lowest method in the class) and the higher order proposed method based interpolation (Lobatto-Gauss-Legendre of order 4 and 6). The right plot shows the corresponding error in the Hamiltonian energy.
When the modified frequencyω is such that cos(hω) = ±1, see Figure 1, left plot, we expect to observe resonances. This happens for hω/π = 2 √ 3/π ≈ 1.1 for order four method and for hω/π = √ 10/π ≈ 1 and hω/π = 2 √ 15/π ≈ 2.47 for the order six method. Resonances can be observed in the preservation of the Hamiltonian (total) energy of the system and in the scaled total oscillatory energy ωI in the range (0, 4.5π], the latter being more uniform in dealing with the frequencies. It is clear that the width of the resonance region is inversely proportional to the curvature at the resonance point. The flatter the stability function is at the resonance points in Figure 1, the wider the region of resonance. The left plot in Figure 5 displays the solution obtained by the methods by taking a relatively large step size, with hω/π ≈ 1.59. The approximations to the solutions are still fairly acceptable and the methods do not display excessive oscillations as other trigonometric integrators do.  The right plot in Figure 5 depicts the behavior of the methods as they approach their high-frequency limit, in a similar experiment as in [8]. We keep h = 0.1 but take ω = 1000 with a ratio hω/π ≈ 31.8. In this experiment ω is scaled by a factor of 20 (compared to 200 in [8]) and the time interval must also be scaled correspondingly to [0,4000].
A rough analysis of the slow energy exchange can be performed by using the modified trigonometric integrator form of the method and the expansion of the exact and numerical solution using modulate Fourier expansions, see [3,8]. One difficulty with respect to the standard analysis using modified trigonometric integrators is the presence of more filter functions ψ in (49) and of the internal stages of the methods. However, performing a Taylor expansion of the internal stages, one can  put the methods in the form and apply the standard analysis as for trigonometric integrators.
Using the same setup as in [3,8], one finds that α = 1 + 1 6 µ 2 + O(µ 4 ), β = 1, and γ = 1 − 1 1728 µ 6 + O(µ 8 ). In order to preserve the slow energy exchange at a correct rate, it is required that α = β = γ = 1, a property that is satisfied only by the IMEX, as proven in [8]. It is in particular the value of α that has the strongest effect on the slow energy exchange. Nevertheless, the methods perform way better than classical trigonometric integrators.  In the convergence region we would expect that the ratio would be 16 for the method of order 4 and 64 for the method of order 6, however the plot does not cover well the convergence region. Overall, we see that the methods have a conservation of O(h 2 ), except from the regions corresponding to resonances. This behavior seems to indicate that the methods suffer of order reduction, a phenomenon that is not uncommon for higher order methods in prescribed regions of the step-size. This effect will be discussed more thoroughly below. 6.2. Order reduction. Ultimately, it is the error in the slow variables one of the most relevant quantities in the numerical simulations of these kind of problems, because the fast variables will be in any case poorly resolved. In figures (7-8) we show the errors in the slow variables for the FPUT problem for different values of the step-size and different ω. The errors are evaluated at T = 3 and the exact solution is computed using Matlab's ode45 to about machine precision (setting AbsTol, RelTol = 1e-14). It is observed that the methods suffer of order reduction both in the positions and the momenta, manifested as a platou in the error plots.
In figures (9)-(10) we repeat the same experiments by methods of order 4 and 6 obtained from the IMEX and using the Yoshida technique [16]. Also in this case one can observe an order reduction, from order 4 to order 3 for the positions and from order 4 to 2 for the momenta for the method of order 4. Similarly, one observes a reduction from order 6 to order 3 for the positions and from order 6 to order 2 for the momenta for the method of order 6. In summary, the order reduction is similar to that of the Lobatto-Gauss-Legendre on the momenta, but is one order less on the positions.
It is not clear why the Yoshida technique gives a lesser order reduction for the positions and marginally also for the momenta. We conjecture that it might be due to the fact that the method uses step sizes αh and βh, rather than just h, and the  use of these two step sizes might reduce the resonance effects of the single step size. (12) show a comparison of the errors for methods of the same order. It is observed that for larger step-sizes, the Lobatto-Gauss-Legendre methods have smaller error (about two orders of magnitude) than IMEX with Yoshida timestepping. For smaller step-sizes, there is no obvious answer and the choice of the method will most likely depend on the application under consideration.

Figures (11)-
The Lobatto-Gauss-Legendre of order 4 and 6 have implicit stages, which require one and two functions evaluations respectively. In our numerical experiments, we have solved the implicit stages by fixed point iteration. The number of function evaluations will then depend on the number of fixed point iterations. For small step sizes, we observed 1-2 fixed point iterations. For larger step sizes (but still in the convergence region) we never observed more than 10 iterations, a typical number was ≈5-6. In comparison, the order 4 IMEX with the Yoshida technique would require 3 function evaluations and 9 function evaluation for order 6. However, the Yoshida techniques have larger error in the regions of convergence, especially in the larger step size regions. This error is about two-three orders of magnitude larger Figure 10. Error in the slow positions (left) and slow momenta (right) for the IMEX method with a Yoshida time stepping yielding a method of order 6. Also in this case there is an observable reduction in the order. We have four orders loss for the momenta and three order loss in the positions. than the Lobatto-Gauss-Legendre methods, indicating that these can be used with a larger step size, resulting in an overall cheaper method.
For completeness, the numerical experiments in Figures 3-6 are repeated for IMEX and its order four and six implementations using the Yoshida time stepping technique.
7. Conclusions and further remarks. We have introduced a family of symplectic methods based on a variational derivation. The main idea is to use different integration quadrature formulas for different terms of the Lagrangian. The introduction of extra internal stages is solved either by interpolation or by collocation. In particular, we have derived a higher order generalization of the IMEX method (using the Verlet method and an interpolated form of the Implicit Midpoint Rule), namely the LobattoIIIA-B-Gauss-Legendre family of arbitrary order, and present the coefficients explicitly for the methods of order 4 and 6. We have proved that these method possess the expected order and shown that the methods with internal   Figure 3, this time using the IMEX and its order four and six implementation using the Yoshida technique. The Hamiltonian error (right) does not separate as clearly as the Lobatto-Gauss methods of corresponding order.
stages solved by interpolation are P-stable, making these particularly interesting in the context of oscillatory problem. We have also observed that these higher order methods might suffer from resonance and from order reduction. The methods are thoroughly tested on the FPUT problem and their behavior is compared to higher order IMEX implementations using the Yoshida time-stepping technique.
The proposed methods might be considered as special subclass of additive Runge-Kutta methods (ARK). The advantage of the variational derivation is that the methods are automatically symplectic, therefore particularly suited to geometric integration. It will be interesting to explore further this mixed technique for other choices of primary/secondary methods and the use other techniques, like treating some of the terms by averaged Lagrangian methods in the spirit of [1]. Possibly, this mixed approach might lead to further interesting numerical method that might not be easily discovered using the classical algebraic theory of RK and ARK methods.