Classical 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 implicitness as the underlying Lobatto IIIA-B pair, and, in addition, they are P-stable, therefore suitable for application to highly oscillatory problems.
Citation: |
Figure 1. Left: Plot of the stability functions for the the Lobatto IIIA-B and Gauss–Legendre combinations of order two (IMEX), order four and order six with coefficients constructed by interpolation $ (20) $. These methods are P-stable. Right: Stability function plot for the methods with coefficients constructed by collocation $ (21) $. For P-stability, the function must have values between $ -1 $ and $ 1 $ for all $ \mu $. These methods are not P-stable. See text for their interval of stability
Figure 2. Modified frequency for the methods in the Lobatto IIIA and Gauss—Legendre family of order two (IMEX), four and six. The straight line is the identity function, for which $ \tilde \lambda = \lambda $. The IMEX retains the correct frequency of oscillations up to $ h\lambda\approx 1 $. The order two method retains the correct frequency up to $ h\lambda\approx 2 $, while the order six method up to $ h\lambda \approx 3 $
Figure 4. Left: Maximum deviation in the Hamiltonian (total) energy error. Right: Maximum deviation in scaled oscillatory energy $ \omega I $ error. The computation is performed in [0, 200] for $ h\omega/\pi \in ( 0, 4.5] $, $ h = 0.02 $. The peaks correspond to the resonances of the methods. These occur when $ \cos(h\tilde \omega) = \pm1 $, namely when $ h\omega/\pi \approx 1.1 $ for the order four method and when $ h\omega/\pi \approx 1, 2.47 $ for the order six method. See text for details
Figure 5. Left: Individual oscillatory energies $ I_i $ and total oscillatory energy $ I = \sum_i I_i $, as in Figure 3, with step size $ h = 0.1 $, $ \omega = 50 $. Right: Same oscillatory energies as in the left plot. Here the step size is kept fixed to $ h = 0.1 $ but the frequency $ \omega $ as well as the length of the interval is scaled by a factor of $ 20 $
Figure 13. Same parameters and initial conditions as in 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
Figure 14. Comparison of the Hamiltonian error for the Lobatto-Gauss methods (left) and Yoshida methods of the same order (right), same parameters and initial conditions as in Figure 3, except for the step-size, $ h = 0.005 $. It is known that the Yoshida technique can lead to a large error constant in the leading error term
Figure 15. Left: Maximum error in Hamiltonian (total) energy. Right: Maximum deviation in scaled oscillatory energy $ \omega I $. Same parameters and initial conditions as in Figure 4, this time with the IMEX and its order four and six implementations using the Yoshida technique
Figure 16. Same parameters and initial conditions as in Figure 5, this time with the IMEX and its order four and six implementations using the Yoshida technique
Figure 17. Hamiltonian max error ratio for the IMEX and its higher order implementations with the Yoshida technique. Same initial conditions and parameters as in Figure 6
[1] | E. Celledoni and E. H. Høyseth, The averaged Lagrangian method, J. Comput. Appl. Math., 316 (2017), 161-174. doi: 10.1016/j.cam.2016.09.047. |
[2] | G. J. Cooper and A. Sayfy, Additive Runge–Kutta methods for stiff ordinary differential equations, Math. Comp., 40 (1983), 207-218. doi: 10.1090/S0025-5718-1983-0679441-1. |
[3] | E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration, Structure-Preserving Algorithms for Ordinary Differential Equations, Springer Series in Computational Mathematics, Springer, 2010. |
[4] | L. Jay, Specialized partitioned additive Runge–Kutta methods for systems of overdetermined DAEs with holonomic constraints, SIAM J. Numer. Anal., 45 (2007), 1814-1842. doi: 10.1137/060667475. |
[5] | L. Jay, Structure preservation for constrained dynamics with super partitioned additive Runge–Kutta methods, SIAM J. Sci. Comput., 20 (1998), 416-446. doi: 10.1137/S1064827595293223. |
[6] | L. O. Jay and L. R. Petzold, Highly oscillatory systems and periodic stability, Technical Report, Army High Performance Computing Research Center, Stanford, CA, (1995), 95–105. |
[7] | J. E. Marsden and M. West, Discrete mechanics and variational integrators, Acta Numerica, 10 (2001), 357-514. doi: 10.1017/S096249290100006X. |
[8] | R. I. McLachlan and A. Stern, Modified trigonometric integrators, SIAM J. Numer. Anal., 52 (2014), 1378-1397. doi: 10.1137/130921118. |
[9] | R. I. McLachlan, Y. Sun and P. S. P. Tse, Linear stability of partitioned Runge-Kutta methods, SIAM J. Numer. Anal., 49 (2011), 232-263. doi: 10.1137/100787234. |
[10] | F. Pfeil, A Higher Order IMEX Method for Solving Highly Oscillatory Problems, Master's thesis, University of Bergen, Norway, June 2019. |
[11] | A. Sandu and M. Günther, A generalized-structure approach to additive Runge–Kutta methods, SIAM J. Numer. Anal., 53 (2015), 17-42. doi: 10.1137/130943224. |
[12] | A. Stern and E. Grinspun, Implicit-explicit variational integration of highly oscillatory problems, Multiscale Model. Simul., 7 (2009), 1779-1794. doi: 10.1137/080732936. |
[13] | G. M. Tanner, Generalized Additive Runge–Kutta Methods for Stiff Odes, PhD thesis, University of Iowa, 2018. |
[14] | T. Wenger, S. Ober-Blöbaum and S. Leyendecker, Variational integrators of mixed order for dynamical systems with multiple time scales and split potentials, ECCOMAS Congress 2016, 2016. doi: 10.7712/100016.1920.10163. |
[15] | T. Wenger, S. Ober-Blöbaum and S. Leyendecker, Construction and analysis of higher order variational integrators for dynamical systems with holonomic constraints, Adv. Comput. Math., 43 (2017), 1163-1195. doi: 10.1007/s10444-017-9520-5. |
[16] | H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A, 150 (1990), 262-268. doi: 10.1016/0375-9601(90)90092-3. |
[17] | A. Zanna, A family of modified trigonometric integrators for highly oscillatory problems, FoCM, Barcelona, 2017. |
[18] | M. Zhang and R. D. Skeel, Cheap implicit symplectic integrators, Appl. Numer. Math., 25 (1997), 297-302. doi: 10.1016/S0168-9274(97)00066-4. |
Left: Plot of the stability functions for the the Lobatto IIIA-B and Gauss–Legendre combinations of order two (IMEX), order four and order six with coefficients constructed by interpolation
Modified frequency for the methods in the Lobatto IIIA and Gauss—Legendre family of order two (IMEX), four and six. The straight line is the identity function, for which
Left: Individual oscillatory energies
Left: Maximum deviation in the Hamiltonian (total) energy error. Right: Maximum deviation in scaled oscillatory energy
Left: Individual oscillatory energies
Hamiltonian max error ratio for
Errors at
Error in the slow positions (left) and slow momenta (right) for the Lobatto-Gauss method of order 6
Error in the slow positions (left) and slow momenta (right) for the IMEX method with a Yoshida time stepping for a method of order 4. There is a order two reduction in the momenta, but only an order one reduction for the error in the slow positions
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
Comparison of the error in the slow positions (left) and slow momenta for the interpolation Lobatto—Gauss (solid line) and the IMEX-Yoshida method (dashed line) of order four
Comparison of the error in the slow positions (left) and slow momenta for the interpolation Lobatto—Gauss (solid line) and the IMEX-Yoshida method (dashed line) of order six
Same parameters and initial conditions as in 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
Comparison of the Hamiltonian error for the Lobatto-Gauss methods (left) and Yoshida methods of the same order (right), same parameters and initial conditions as in Figure 3, except for the step-size,
Left: Maximum error in Hamiltonian (total) energy. Right: Maximum deviation in scaled oscillatory energy
Same parameters and initial conditions as in Figure 5, this time with the IMEX and its order four and six implementations using the Yoshida technique
Hamiltonian max error ratio for the IMEX and its higher order implementations with the Yoshida technique. Same initial conditions and parameters as in Figure 6