# American Institute of Mathematical Sciences

April  2022, 9(2): 299-328. doi: 10.3934/jcd.2021030

## Symplectic P-stable additive Runge—Kutta methods

 Department of Mathematics, University of Bergen, Postbox 7803, 5020 Bergen, Norway

Received  December 2020 Revised  July 2021 Published  April 2022 Early access  December 2021

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: Antonella Zanna. Symplectic P-stable additive Runge—Kutta methods. Journal of Computational Dynamics, 2022, 9 (2) : 299-328. doi: 10.3934/jcd.2021030
##### References:
 [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.

show all references

##### References:
 [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 $(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
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$
Left: Individual oscillatory energies $I_i$ and total oscillatory energy $I = \sum_i I_i$. Right: Energy error $|H-H_0|$. The simulations are performed in $[0, 200]$ with $\omega = 50$ and $h = 2/\omega = 0.04$. See text for initial conditions
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
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$
Hamiltonian max error ratio for $h = 0.04$ and $h = 0.02$, $T = 200$. We observe peaks in correspondence of the resonances
Errors at $T = 3$ in the slow positions (left) and slow momenta (right) for the Lobatto-Gauss method of order 4 against the step size $h$ for $\omega = 10, \ldots, 10^4$. The lines for $h^2$ and $h^4$ are plotted for convenience
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, $h = 0.005$. It is known that the Yoshida technique can lead to a large error constant in the leading error term
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
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
 [1] Sihong Shao, Huazhong Tang. Higher-order accurate Runge-Kutta discontinuous Galerkin methods for a nonlinear Dirac model. Discrete and Continuous Dynamical Systems - B, 2006, 6 (3) : 623-640. doi: 10.3934/dcdsb.2006.6.623 [2] Elisa Giesecke, Axel Kröner. Classification with Runge-Kutta networks and feature space augmentation. Journal of Computational Dynamics, 2021, 8 (4) : 495-520. doi: 10.3934/jcd.2021018 [3] Da Xu. Numerical solutions of viscoelastic bending wave equations with two term time kernels by Runge-Kutta convolution quadrature. Discrete and Continuous Dynamical Systems - B, 2017, 22 (6) : 2389-2416. doi: 10.3934/dcdsb.2017122 [4] Wenjuan Zhai, Bingzhen Chen. A fourth order implicit symmetric and symplectic exponentially fitted Runge-Kutta-Nyström method for solving oscillatory problems. Numerical Algebra, Control and Optimization, 2019, 9 (1) : 71-84. doi: 10.3934/naco.2019006 [5] Lijin Wang, Jialin Hong. Generating functions for stochastic symplectic methods. Discrete and Continuous Dynamical Systems, 2014, 34 (3) : 1211-1228. doi: 10.3934/dcds.2014.34.1211 [6] Juan Carlos Marrero, David Martín de Diego, Ari Stern. Symplectic groupoids and discrete constrained Lagrangian mechanics. Discrete and Continuous Dynamical Systems, 2015, 35 (1) : 367-397. doi: 10.3934/dcds.2015.35.367 [7] Per Christian Moan, Jitse Niesen. On an asymptotic method for computing the modified energy for symplectic methods. Discrete and Continuous Dynamical Systems, 2014, 34 (3) : 1105-1120. doi: 10.3934/dcds.2014.34.1105 [8] Yanyan Shi, Yajuan Sun, Yulei Wang, Jian Liu. Study of adaptive symplectic methods for simulating charged particle dynamics. Journal of Computational Dynamics, 2019, 6 (2) : 429-448. doi: 10.3934/jcd.2019022 [9] Richard A. Norton, David I. McLaren, G. R. W. Quispel, Ari Stern, Antonella Zanna. Projection methods and discrete gradient methods for preserving first integrals of ODEs. Discrete and Continuous Dynamical Systems, 2015, 35 (5) : 2079-2098. doi: 10.3934/dcds.2015.35.2079 [10] Holger Heumann, Ralf Hiptmair. Eulerian and semi-Lagrangian methods for convection-diffusion for differential forms. Discrete and Continuous Dynamical Systems, 2011, 29 (4) : 1471-1495. doi: 10.3934/dcds.2011.29.1471 [11] Andrew J. Steyer, Erik S. Van Vleck. Underlying one-step methods and nonautonomous stability of general linear methods. Discrete and Continuous Dynamical Systems - B, 2018, 23 (7) : 2859-2877. doi: 10.3934/dcdsb.2018108 [12] Robert I. McLachlan, G. R. W. Quispel. Discrete gradient methods have an energy conservation law. Discrete and Continuous Dynamical Systems, 2014, 34 (3) : 1099-1104. doi: 10.3934/dcds.2014.34.1099 [13] Zimo Zhu, Gang Chen, Xiaoping Xie. Semi-discrete and fully discrete HDG methods for Burgers' equation. Communications on Pure and Applied Analysis, , () : -. doi: 10.3934/cpaa.2021132 [14] Xuefeng Shen, Khoa Tran, Melvin Leok. High-order symplectic Lie group methods on $SO(n)$ using the polar decomposition. Journal of Computational Dynamics, 2022  doi: 10.3934/jcd.2022003 [15] Antonia Katzouraki, Tania Stathaki. Intelligent traffic control on internet-like topologies - integration of graph principles to the classic Runge--Kutta method. Conference Publications, 2009, 2009 (Special) : 404-415. doi: 10.3934/proc.2009.2009.404 [16] Richard A. Norton, G. R. W. Quispel. Discrete gradient methods for preserving a first integral of an ordinary differential equation. Discrete and Continuous Dynamical Systems, 2014, 34 (3) : 1147-1170. doi: 10.3934/dcds.2014.34.1147 [17] Diogo A. Gomes. Viscosity solution methods and the discrete Aubry-Mather problem. Discrete and Continuous Dynamical Systems, 2005, 13 (1) : 103-116. doi: 10.3934/dcds.2005.13.103 [18] Li-Xia Liu, Sanyang Liu, Chun-Feng Wang. Smoothing Newton methods for symmetric cone linear complementarity problem with the Cartesian $P$/$P_0$-property. Journal of Industrial and Management Optimization, 2011, 7 (1) : 53-66. doi: 10.3934/jimo.2011.7.53 [19] Xinlong Feng, Huailing Song, Tao Tang, Jiang Yang. Nonlinear stability of the implicit-explicit methods for the Allen-Cahn equation. Inverse Problems and Imaging, 2013, 7 (3) : 679-695. doi: 10.3934/ipi.2013.7.679 [20] Syvert P. Nørsett, Andreas Asheim. Regarding the absolute stability of Størmer-Cowell methods. Discrete and Continuous Dynamical Systems, 2014, 34 (3) : 1131-1146. doi: 10.3934/dcds.2014.34.1131

Impact Factor:

## Tools

Article outline

Figures and Tables
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 1, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.5466118454933167, "perplexity": 2294.9417158716224}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2022-27/segments/1656104660626.98/warc/CC-MAIN-20220706030209-20220706060209-00627.warc.gz"}