Symplectic P-stable additive Runge—Kutta methods

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 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

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 6.  Hamiltonian max error ratio for $ h = 0.04 $ and $ h = 0.02 $, $ T = 200 $. We observe peaks in correspondence of the resonances

Figure 7.  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

Figure 8.  Error in the slow positions (left) and slow momenta (right) for the Lobatto-Gauss method of order 6

Figure 9.  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

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

Figure 11.  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

Figure 12.  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

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