Skip to main content
SpringerLink
Log in

Implicit Multirate GARK Methods

  • Published:
Journal of Scientific Computing Aims and scope Submit manuscript

Abstract

This work considers multirate generalized-structure additively partitioned Runge–Kutta methods for solving stiff systems of ordinary differential equations with multiple time scales. These methods treat different partitions of the system with different timesteps for a more targeted and efficient solution compared to monolithic single rate approaches. With implicit methods used across all partitions, methods must find a balance between stability and the cost of solving nonlinear equations for the stages. In order to characterize this important trade-off, we explore multirate coupling strategies, problems for assessing linear stability, and techniques to efficiently implement Newton iterations for stage equations. Unlike much of the existing multirate stability analysis which is limited in scope to particular methods, we present general statements on stability and describe fundamental limitations for certain types of multirate schemes. New implicit multirate methods up to fourth order are derived, and their accuracy and efficiency properties are verified with numerical tests.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Log in via an institution

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Fig. 1
Fig. 2
Fig. 3

Similar content being viewed by others

Availability of Data and Material

The implementation of the CUSP problem is available in ODE Test Problems: https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems.

References

  1. Alexander, R.: Diagonally implicit Runge–Kutta methods for stiff ODEs. SIAM J. Numer. Anal. 14(6), 1006–1021 (1977). https://doi.org/10.1137/0714068

    Article  MathSciNet  MATH  Google Scholar 

  2. Andrus, J.: Stability of a multi-rate method for numerical integration of ODEs. Comput. Math. Appl. 25(2), 3–14 (1993)

    Article  MathSciNet  Google Scholar 

  3. Ascher, U.M., Ruuth, S.J., Spiteri, R.J.: Implicit-explicit Runge–Kutta methods for time-dependent partial differential equations. Appl. Numer. Math. 25(2), 151–167 (1997). https://doi.org/10.1016/S0168-9274(97)00056-1

    Article  MathSciNet  MATH  Google Scholar 

  4. Bartel, A., Günther, M.: A multirate W-method for electrical networks in state-space formulation. J. Comput. Appl. Math. 147(2), 411–425 (2002). https://doi.org/10.1016/S0377-0427(02)00476-4

    Article  MathSciNet  MATH  Google Scholar 

  5. Beam, R.M., Warming, R.: An implicit finite-difference algorithm for hyperbolic systems in conservation-law form. J. Comput. Phys. 22(1), 87–110 (1976). https://doi.org/10.1016/0021-9991(76)90110-8

    Article  MathSciNet  MATH  Google Scholar 

  6. Bonaventura, L., Casella, F., Delpopolo, L., Ranade, A.: A self adjusting multirate algorithm based on the TR-BDF2 method. (2018). arXiv preprint arXiv:1801.09118

  7. Constantinescu, E., Sandu, A.: Multirate timestepping methods for hyperbolic conservation laws. J. Sci. Comput. 33(3), 239–278 (2007). https://doi.org/10.1007/s10915-007-9151-y

    Article  MathSciNet  MATH  Google Scholar 

  8. Constantinescu, E.M., Sandu, A.: Extrapolated multirate methods for differential equations with multiple time scales. J. Sci. Comput. 56(1), 28–44 (2013). https://doi.org/10.1007/s10915-012-9662-z

    Article  MathSciNet  MATH  Google Scholar 

  9. Carciopolo, L.D., Bonaventura, L., Scotti, A., Formaggia, L.: A conservative implicit multirate method for hyperbolic problems. Comput. Geosci. (2018). https://doi.org/10.1007/s10596-018-9764-2

  10. Gear, C.: Multirate methods for ordinary differential equations. Tech. rep., Illinois Univ., Urbana (USA). Dept. of Computer Science (1974)

  11. Gear, C.W., Wells, D.: Multirate linear multistep methods. BIT Numer. Math. 24(4), 484–502 (1984)

    Article  MathSciNet  Google Scholar 

  12. Günther, M., Hoschek, M.: ROW methods adapted to electric circuit simulation packages. J. Comput. Appl. Math. 82(1–2), 159–170 (1997). https://doi.org/10.1016/S0377-0427(97)00043-5

    Article  MathSciNet  MATH  Google Scholar 

  13. Günther, M., Rentrop, P.: Multirate ROW methods and latency of electric circuits. Appl. Numer. Math. 13(1), 83–102 (1993). https://doi.org/10.1016/0168-9274(93)90133-C

    Article  MathSciNet  MATH  Google Scholar 

  14. Günther, M., Sandu, A.: Multirate generalized additive Runge–Kutta methods. Numer. Math. 133(3), 497–524 (2016). https://doi.org/10.1007/s00211-015-0756-z

    Article  MathSciNet  MATH  Google Scholar 

  15. Hachtel, C., Bartel, A., Günther, M., Sandu, A.: Multirate implicit Euler schemes for a class of differential-algebraic equations of index-1. J. Comput. Appl. Math. p. 112499 (2019). https://doi.org/10.1016/j.cam.2019.112499

  16. Hairer, E., Wanner, G.: Solving ordinary differential equations II: stiff and differential-algebraic problems, 2 edn. No. 14 in Springer Series in Computational Mathematics. Springer, Berlin (1996)

  17. Hairer, E., Wanner, G., Lubich, C.: Geometric Numerical Integration. Springer, Berlin (2006). https://doi.org/10.1007/3-540-30666-8

  18. Hundsdorfer, W., Savcenco, V.: Analysis of a multirate theta-method for stiff ODEs. Appl. Numer. Math. 59(3–4), 693–706 (2009). https://doi.org/10.1016/j.apnum.2008.03.022

    Article  MathSciNet  MATH  Google Scholar 

  19. Kennedy, C.A., Carpenter, M.H.: Additive Runge–Kutta schemes for convection–diffusion–reaction equations. Appl. Numer. Math. 44(1), 139–181 (2003). https://doi.org/10.1016/S0168-9274(02)00138-1

    Article  MathSciNet  MATH  Google Scholar 

  20. Knoth, O., Wolke, R.: Implicit-explicit Runge-Kutta methods for computing atmospheric reactive flows. Appl. Numer. Math. 28(327–341) (1998)

  21. Kværnø, A.: Stability of multirate Runge-Kutta schemes. Int. J. Differ. Equ. Appl. 1(1) 97–105 (2000)

  22. Kværnø, A., Rentrop, P.: Low order multirate Runge-Kutta methods in electric circuit simulation, Preprint 99/1. University of Karlsruhe, IWRMM, Karlsruhe (1999)

  23. Peaceman, D.W., Rachford Jr., H.H.: The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 3(1), 28–41 (1955). https://doi.org/10.1137/0103003

    Article  MathSciNet  MATH  Google Scholar 

  24. Ralston, A.: Runge–Kutta methods with minimum error bounds. Math. Comput. 16(80), 431–437 (1962). https://doi.org/10.2307/2003133

    Article  MathSciNet  MATH  Google Scholar 

  25. Rice, J.R.: Split Runge–Kutta methods for simultaneous equations. J. Res. Natl. Inst. Stand. Technol. 60 (1960)

  26. Roberts, S., Popov, A.A., Sandu, A.: ODE test problems: a MATLAB suite of initial value problems. arXiv preprint arXiv:1901.04098 (2019). https://github.com/ComputationalScienceLaboratory/ODE-Test-Problems

  27. Roberts, S., Sarshar, A., Sandu, A.: Coupled multirate infinitesimal GARK schemes for stiff systems with multiple time scales. SIAM J. Sci. Comput. 42(3), A1609–A1638 (2020). https://doi.org/10.1137/19M1266952

    Article  MathSciNet  MATH  Google Scholar 

  28. Rodríguez-Gómez, G., González-Casanova, P., Martínez-Carballido, J.: Computing general companion matrices and stability regions of multirate methods. Int. J. Numer. Methods Eng. 61(2), 255–273 (2004)

    Article  MathSciNet  Google Scholar 

  29. Sand, J., Skelboe, S.: Stability of backward Euler multirate methods and convergence of waveform relaxation. BIT Numer. Math. 32(2), 350–366 (1992)

    Article  MathSciNet  Google Scholar 

  30. Sandu, A.: A class of multirate infinitesimal GARK methods. SIAM J. Numer. Anal. 57(5), 2300–2327 (2019). https://doi.org/10.1137/18M1205492

    Article  MathSciNet  MATH  Google Scholar 

  31. Sandu, A., Günther, M.: A generalized-structure approach to additive Runge–Kutta methods. SIAM J. Numer. Anal. 53(1), 17–42 (2015). https://doi.org/10.1137/130943224

    Article  MathSciNet  MATH  Google Scholar 

  32. Sarshar, A., Roberts, S., Sandu, A.: Design of high-order decoupled multirate GARK schemes. SIAM J. Sci. Comput. 41(2), A816–A847 (2019). https://doi.org/10.1137/18M1182875

    Article  MathSciNet  MATH  Google Scholar 

  33. Savcenco, V.: Comparison of the asymptotic stability properties for two multirate strategies. J. Comput. Appl. Math. 220(1–2), 508–524 (2008)

    Article  MathSciNet  Google Scholar 

  34. Savcenco, V.: Construction of a multirate RODAS method for stiff ODEs. J. Comput. Appl. Math. 225(2), 323–337 (2009). https://doi.org/10.1016/j.cam.2008.07.041

    Article  MathSciNet  MATH  Google Scholar 

  35. Savcenco, V., Hundsdorfer, W., Verwer, J.: A multirate time stepping strategy for stiff ordinary differential equations. BIT Numer. Math. 47(1), 137–155 (2007)

    Article  MathSciNet  Google Scholar 

  36. Scheidemann, V.: Introduction to Complex Analysis in Several Variables. Springer, Berlin (2005)

    MATH  Google Scholar 

  37. Schlegel, M., Knoth, O., Arnold, M., Wolke, R.: Multirate Runge–Kutta schemes for advection equations. J. Comput. Appl. Math. 226, 345–357 (2009)

    Article  MathSciNet  Google Scholar 

  38. Sexton, J.M., Reynolds, D.R.: Relaxed multirate infinitesimal step methods. arXiv preprint arXiv:1808.03718 (2018). Submitted to J. Comput. Appl. Math

  39. Skelboe, S., Andersen, P.U.: Stability properties of backward Euler multirate formulas. SIAM J. Sci. Stat. Comput. 10(5), 1000–1009 (1989)

    Article  MathSciNet  Google Scholar 

  40. Strang, G.: On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5(3), 506–517 (1968)

    Article  MathSciNet  Google Scholar 

  41. Verhoeven, A., Beelen, T., El Guennouni, A., Maten, ter, E., Mattheij, R., Tasic, B.: Error analysis of BDF compound-fast multirate method for differential-algebraic equations, CASA-report, vol. 0610. Technische Universiteit Eindhoven (2006)

  42. Verhoeven, A., El Guennouni, A., Ter Maten, E., Mattheij, R.: A general compound multirate method for circuit simulation problems. In: Scientific Computing in Electrical Engineering, pp. 143–149. Springer, Berlin (2006)

  43. Verhoeven, A., Ter Maten, E.J.W., Mattheij, R.M., Tasić, B.: Stability analysis of the BDF slowest-first multirate methods. Int. J. Comput. Math. 84(6), 895–923 (2007)

    Article  MathSciNet  Google Scholar 

  44. Wensch, J., Knoth, O., Galant, A.: Multirate infinitesimal step methods for atmospheric flow simulation. BIT Numer. Math. 49(2), 449–473 (2009)

    Article  MathSciNet  Google Scholar 

  45. Zanna, A.: Discrete variational methods and symplectic generalized additive Runge–Kutta methods. arXiv preprint arXiv:2001.07185 (2020)

  46. Zhao, D., Li, Y., Barbič, J.: Asynchronous implicit backward Euler integration. In: Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’16, pp. 1–9. Eurographics Association, Goslar Germany, Germany (2016). https://doi.org/10.2312/sca.20161217

Download references

Acknowledgements

The authors acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper. We also thank Jeffrey Hittinger and Valentin Dallerit for helpful discussions.

Funding

The work of S. Roberts (in part), J. Loffeld, and C.S Woodward was supported by the Lawrence Livermore Laboratory Directed Research and Development Program under tracking number 17-ERD-035. Their work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. Lawrence Livermore National Security, LLC. LLNL-JRNL-795454. The work of S. Roberts (in part), A. Sarshar, and A. Sandu was supported by the National Science Foundation (awards CCF–1613905 and ACI–1709727), Air Force Office of Scientific Research (award DDDAS FA9550-17-1-0015), and by the Computational Science Laboratory at Virginia Tech. The work of S. Roberts (in part) was supported by the Virginia Space Grant Consortium.

This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.

Author information

Authors and Affiliations

  1. Computational Science Laboratory, Department of Computer Science, Virginia Tech, Blacksburg, USA

    Steven Roberts, Arash Sarshar & Adrian Sandu

  2. Center for Applied Scientific Computing, Lawrence Livermore National Laboratory, Livermore, USA

    John Loffeld & Carol S. Woodward

Authors
  1. Steven Roberts
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. John Loffeld
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Arash Sarshar
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Carol S. Woodward
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Adrian Sandu
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Steven Roberts.

Ethics declarations

Conflict of interest

The authors declare that they have no conflict of interest.

Code Availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendices

Appendix A: Third Order Compound-Fast MrGARK

A third order compound-fast method, which we will refer to as compound-fast MrGARK SDIRK3, is built on the following base method of Alexander [1]:

(A.1)

Here, \(\gamma \approx 0.44\) is the middle root of \(6 \gamma ^3-18 \gamma ^2+9 \gamma -1 = 0\). The coupling coefficients are

$$\begin{aligned} a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,1}&= \frac{\left( 6 \gamma ^2-24 \gamma +5\right) \left( 2 \gamma ^2+2 \gamma (\lambda -1)+(\lambda -1)^2\right) -8 \gamma ^3 M^2}{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad -\frac{\left( 6 \gamma ^3-30 \gamma ^2-15 \gamma +5\right) M (\gamma +\lambda -1)}{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,2}&= \frac{-5 (\lambda -1)^2+12 \gamma ^4 (M-1)+4 \gamma ^3 (M-1) (3 \lambda +4 M-17)}{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad + \frac{\gamma ^2 \left( -6 \lambda ^2+68 \lambda +(82-72 \lambda ) M-72\right) +2 \gamma (\lambda -1) (14 \lambda +5 M-19)}{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,3}&= -\frac{4 \gamma \left( (\lambda -1)^2+2 \gamma ^2 (M-1)^2-2 \gamma (\lambda -1) (2 M-1)\right) }{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,1}&= -\frac{-2 \left( 6 \gamma ^2-24 \gamma +5\right) (\gamma (\lambda +1)+(\lambda -1) \lambda )+16 \gamma ^3 M^2}{2 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad - \frac{\left( 6 \gamma ^3-30 \gamma ^2-15 \gamma +5\right) M (\gamma +2 \lambda -1)}{2 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,2}&= \frac{\gamma \left( 28 \lambda ^2-33 \lambda +5 (2 \lambda -1) M-5\right) + 2 \gamma ^3 \left( 8 M^2-3 (\lambda +1)+3 (2 \lambda -7) M\right) }{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad + \frac{6 \gamma ^4 M -5 (\lambda -1) \lambda +\gamma ^2 \left( -6 \lambda ^2+34 \lambda +(41-72 \lambda ) M+28\right) }{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,3}&= -\frac{4 \gamma \left( (\lambda -1) \lambda +2 \gamma ^2 (M-1) M+\gamma (\lambda +(2-4 \lambda ) M+1)\right) }{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,1}&= \frac{-36 \gamma ^5+252 \gamma ^4+20 \lambda ^2-4 \gamma ^3 \left( 8 M^2+6 \lambda M+129\right) }{4 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad + \frac{24 \gamma ^2 \left( \lambda ^2+5 \lambda M+13\right) +\gamma \left( -96 \lambda ^2+60 \lambda M-69\right) -20 \lambda M+5}{4 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,2}&= \frac{36 \gamma ^5-276 \gamma ^4-5 \left( 4 \lambda ^2+1\right) +\gamma ^3 \left( 64 M^2+48 \lambda M+588\right) }{4 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2} \nonumber \\&\quad + \frac{-12 \gamma ^2 \left( 2 \lambda ^2+24 \lambda M+29\right) +\gamma \left( 112 \lambda ^2+40 \lambda M+73\right) }{4 (\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,3}&= \frac{\gamma \left( 6 \gamma ^3-4 \lambda ^2-2 \gamma ^2 \left( 4 M^2+9\right) +\gamma (16 \lambda M+9)-1\right) }{(\gamma -1) \left( 6 \gamma ^2-20 \gamma +5\right) M^2}. \end{aligned}$$
(A.2)

Appendix B: Fourth Order Compound-Fast MrGARK

For the compound-fast method of order four, we start with a new base method, solving the coupling and base conditions together. This allows more flexibility to keep the coupling coefficients bounded functions of \(\lambda \) and M. The following L-stable base method was derived:

(B.1)

When paired with the following coupling coefficients, we have the compound-fast MrGARK SDIRK4 scheme:

$$\begin{aligned} a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,1}&{=} \frac{{-}165 \left( 64 \lambda ^4{-}192 \lambda ^3{+}240 \lambda ^2{-}148 \lambda {+}37\right) {+}2688 (4 \lambda {-}3) M^3{-}3408 \left( 8 \lambda ^2{-}12 \lambda {+}5\right) M^2{+}448 \left( 64 \lambda ^3{-}144 \lambda ^2{+}120 \lambda {-}37\right) M}{1836 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,2}&{\!=\!} \frac{1110 \left( 64 \lambda ^4{-}192 \lambda ^3{\!+\!}240 \lambda ^2{\!-\!}148 \lambda {\!+\!}37\right) {\!-\!}1296 M^4{-}2292 (4 \lambda {-}3) M^3{\!+\!}8781 \left( 8 \lambda ^2{\!-\!}12 \lambda {\!+\!}5\right) M^2{\!-\!}2101 \left( 64 \lambda ^3{\!-\!}144 \lambda ^2{\!+\!}120 \lambda {\!-\!}37\right) M}{22464 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,3}&{=} {-}\frac{125 \left( {-}33 \left( 64 \lambda ^4{-}192 \lambda ^3{+}240 \lambda ^2{-}148 \lambda {+}37\right) {+}336 (4 \lambda {-}3) M^3{-}552 \left( 8 \lambda ^2{-}12 \lambda {+}5\right) M^2{+}83 \left( 64 \lambda ^3{-}144 \lambda ^2{+}120 \lambda {-}37\right) M\right) }{22464 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,4}&{=} \frac{1331 \left( {-}5 \left( 64 \lambda ^4{-}192 \lambda ^3{+}240 \lambda ^2{-}148 \lambda {+}37\right) {+}32 (4 \lambda {-}3) M^3{-}60 \left( 8 \lambda ^2{-}12 \lambda {+}5\right) M^2{+}11 \left( 64 \lambda ^3{-}144 \lambda ^2{+}120 \lambda {-}37\right) M\right) }{56576 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{1,5}&\!=\! \frac{\!-\!85 \left( 64 \lambda ^4\!-\!192 \lambda ^3\!+\!240 \lambda ^2\!-\!148 \lambda \!+\!37\right) \!+\!192 M^4\!+\!16 (4 \lambda -3) M^3\!-\!648 \left( 8 \lambda ^2-12 \lambda \!+\!5\right) M^2\!\!+\!\!175 \left( 64 \lambda ^3\!-\!144 \lambda ^2\!+\!120 \lambda \!-\!37\right) M}{3328 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,1}&= \frac{-165 \left( 64 \lambda ^4-48 \lambda ^2+68 \lambda -17\right) +10752 \lambda M^3-3408 \left( 8 \lambda ^2-1\right) M^2+448 \left( 64 \lambda ^3-24 \lambda +17\right) M}{1836 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,2}&= \frac{1110 \left( 64 \lambda ^4-48 \lambda ^2+68 \lambda -17\right) -1296 M^4-9168 \lambda M^3+8781 \left( 8 \lambda ^2-1\right) M^2-2101 \left( 64 \lambda ^3-24 \lambda +17\right) M}{22464 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,3}&= -\frac{125 \left( -33 \left( 64 \lambda ^4-48 \lambda ^2+68 \lambda -17\right) +1344 \lambda M^3-552 \left( 8 \lambda ^2-1\right) M^2+83 \left( 64 \lambda ^3-24 \lambda +17\right) M\right) }{22464 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,4}&= \frac{1331 \left( 5 \left( -64 \lambda ^4+48 \lambda ^2-68 \lambda +17\right) +128 \lambda M^3-60 \left( 8 \lambda ^2-1\right) M^2+11 \left( 64 \lambda ^3-24 \lambda +17\right) M\right) }{56576 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{2,5}&= \frac{-85 \left( 64 \lambda ^4-48 \lambda ^2+68 \lambda -17\right) +192 M^4+64 \lambda M^3-648 \left( 8 \lambda ^2-1\right) M^2+175 \left( 64 \lambda ^3-24 \lambda +17\right) M}{3328 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,1}&= \frac{-33 \left( 32000 \lambda ^4-76800 \lambda ^3+84720 \lambda ^2-56180 \lambda +15773\right) +215040 (5 \lambda -3) M^3}{183600 M^4} \nonumber \\&\quad + \frac{-3408 \left( 800 \lambda ^2-960 \lambda +353\right) M^2+448 \left( 6400 \lambda ^3-11520 \lambda ^2+8472 \lambda -2809\right) M}{183600 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,2}&= \frac{222 \left( 32000 \lambda ^4-76800 \lambda ^3+84720 \lambda ^2-56180 \lambda +15773\right) -129600 M^4-183360 (5 \lambda -3) M^3}{2246400 M^4} \nonumber \\&\quad + \frac{8781 \left( 800 \lambda ^2-960 \lambda +353\right) M^2-2101 \left( 6400 \lambda ^3-11520 \lambda ^2+8472 \lambda -2809\right) M}{2246400 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,3}&= \frac{33 \left( 32000 \lambda ^4-76800 \lambda ^3+84720 \lambda ^2-56180 \lambda +15773\right) -134400 (5 \lambda -3) M^3}{89856 M^4} \nonumber \\&\quad + \frac{2760 \left( 800 \lambda ^2-960 \lambda +353\right) M^2-415 \left( 6400 \lambda ^3-11520 \lambda ^2+8472 \lambda -2809\right) M}{89856 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,4}&= \frac{1331 \left( -32000 \lambda ^4+76800 \lambda ^3-84720 \lambda ^2+56180 \lambda +2560 (5 \lambda -3) M^3\right) }{5657600 M^4} \nonumber \\&\quad + \frac{1331 \left( -60 \left( 800 \lambda ^2-960 \lambda +353\right) M^2+11 \left( 6400 \lambda ^3-11520 \lambda ^2+8472 \lambda -2809\right) M-15773\right) }{5657600 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{3,5}&= \frac{-17 \left( 32000 \lambda ^4-76800 \lambda ^3+84720 \lambda ^2-56180 \lambda +15773\right) +19200 M^4+1280 (5 \lambda -3) M^3}{332800 M^4} \nonumber \\&\quad + \frac{-648 \left( 800 \lambda ^2-960 \lambda +353\right) M^2+175 \left( 6400 \lambda ^3-11520 \lambda ^2+8472 \lambda -2809\right) M}{332800 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{4,1}&= \frac{-33 \left( 425920 \lambda ^4-619520 \lambda ^3+280720 \lambda ^2-35660 \lambda +2387\right) +1300992 (11 \lambda -4) M^3}{2443716 M^4} \nonumber \\&\quad + \frac{-137456 \left( 264 \lambda ^2-192 \lambda +29\right) M^2+448 \left( 85184 \lambda ^3-92928 \lambda ^2+28072 \lambda -1783\right) M}{2443716 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{4,2}&= \frac{222 \left( 425920 \lambda ^4-619520 \lambda ^3+280720 \lambda ^2-35660 \lambda +2387\right) -1724976 M^4-1109328 (11 \lambda -4) M^3}{29899584 M^4} \nonumber \\&\quad + \frac{354167 \left( 264 \lambda ^2-192 \lambda +29\right) M^2-2101 \left( 85184 \lambda ^3-92928 \lambda ^2+28072 \lambda -1783\right) M}{29899584 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{4,3}&= -\frac{25 \left( -33 \left( 425920 \lambda ^4-619520 \lambda ^3+280720 \lambda ^2-35660 \lambda +2387\right) +813120 (11 \lambda -4) M^3\right) }{29899584 M^4} \nonumber \\&\quad - \frac{25 \left( -111320 \left( 264 \lambda ^2-192 \lambda +29\right) M^2+415 \left( 85184 \lambda ^3-92928 \lambda ^2+28072 \lambda -1783\right) M\right) }{29899584 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{4,4}&= \frac{-425920 \lambda ^4+619520 \lambda ^3-280720 \lambda ^2+35660 \lambda +15488 (11 \lambda -4) M^3}{56576 M^4} \nonumber \\&\quad + \frac{-2420 \left( 264 \lambda ^2-192 \lambda +29\right) M^2+11 \left( 85184 \lambda ^3-92928 \lambda ^2+28072 \lambda -1783\right) M-2387}{56576 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{4,5}&= \frac{-17 \left( 425920 \lambda ^4-619520 \lambda ^3+280720 \lambda ^2-35660 \lambda +2387\right) +255552 M^4+7744 (11 \lambda -4) M^3}{4429568 M^4} \nonumber \\&\quad + \frac{-26136 \left( 264 \lambda ^2-192 \lambda +29\right) M^2+175 \left( 85184 \lambda ^3-92928 \lambda ^2+28072 \lambda -1783\right) M}{4429568 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{5,1}&= \frac{16 \lambda \left( -165 \lambda ^3+168 M^3-426 \lambda M^2+448 \lambda ^2 M\right) }{459 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{5,2}&= -\frac{-8880 \lambda ^4+162 M^4+1146 \lambda M^3-8781 \lambda ^2 M^2+16808 \lambda ^3 M}{2808 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{5,3}&= \frac{125 \lambda \left( 33 \lambda ^3-21 M^3+69 \lambda M^2-83 \lambda ^2 M\right) }{351 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{5,4}&= \frac{1331 \lambda \left( -10 \lambda ^3+4 M^3-15 \lambda M^2+22 \lambda ^2 M\right) }{1768 M^4}, \nonumber \\ a^{\left\{ {\mathfrak {f}, \mathfrak {s}}, \lambda \right\} }_{5,5}&= \frac{-85 \lambda ^4+3 M^4+\lambda M^3-81 \lambda ^2 M^2+175 \lambda ^3 M}{52 M^4}. \end{aligned}$$
(B.2)

Appendix C: Conservation of Linear Invariants

For some applications, it is desirable that a numerical integrator uphold invariant properties of the system, such as conservation of mass and energy. These invariants are generally expressed as first integrals of the system. It is well known that non-partitioned explicit and implicit Runge–Kutta methods conserve linear first integrals but must meet certain restrictions to conserve quadratic ones, e.g. as with symplectic methods. Partitioned Runge–Kutta methods must obey additional restrictions to conserve even linear invariants. A detailed discussion about preservation of first integrals by Runge–Kutta methods can be found in [17].

GARK schemes are partitioned methods, and here we briefly describe conditions for them to uphold linear first integrals. Consider an ODE Eq. (1.1) satisfying the following linear invariant:

$$\begin{aligned} \zeta ^T \left( f^{\left\{ {\mathfrak {f}}\right\} }(y) + f^{\left\{ {\mathfrak {s}}\right\} }(y) \right) = 0 \quad \Rightarrow \quad \zeta ^T y(t) = \text {const} \end{aligned}$$
(C.1)

When a GARK method applied to this system, the step update Eq. (2.1) satisfies

$$\begin{aligned} \zeta ^T y_{n+1} = \zeta ^T y_{n} + H \, \sum _{j=1}^{s^{\left\{ {\mathfrak {f}}\right\} }} b^{\left\{ {\mathfrak {f}}\right\} }_j \, \zeta ^T \, f^{\left\{ {\mathfrak {f}}\right\} }\left( Y^{\left\{ {\mathfrak {f}}\right\} }_j \right) + H \, \sum _{j=1}^{s^{\left\{ {\mathfrak {s}}\right\} }} b^{\left\{ {\mathfrak {s}}\right\} }_j \, \zeta ^T \, f^{\left\{ {\mathfrak {s}}\right\} }\left( Y^{\left\{ {\mathfrak {s}}\right\} }_j \right) \end{aligned}$$
(C.2)

In order to apply Eq. (C.1), the arguments of \(f^{\left\{ {\mathfrak {f}}\right\} }\) and \(f^{\left\{ {\mathfrak {s}}\right\} }\) must be identical. In general, the arguments in Eq. (C.2) are different: \(Y^{\left\{ {\mathfrak {f}}\right\} }_j \ne Y^{\left\{ {\mathfrak {s}}\right\} }_j\). Moreover, the number of slow stages can be different than the number of fast stages, which prevents the pairing of terms as in Eq. (C.1). Therefore, a general GARK method cannot be expected to preserve linear invariants.

There are special cases, however, where is is possible. If the subsystems individually satisfy

$$\begin{aligned} \zeta ^T f^{\left\{ {\mathfrak {f}}\right\} }(y) = 0, \quad \text {and} \quad \zeta ^T f^{\left\{ {\mathfrak {s}}\right\} }(y) = 0, \end{aligned}$$

one can see Eq. (C.2) simplifies to \(\zeta ^T y_{n+1} = \zeta ^T y_{n}\). Also, when

$$\begin{aligned} \mathbf {A}^{\left\{ {\mathfrak {f}, \mathfrak {f}}\right\} }= \mathbf {A}^{\left\{ {\mathfrak {s}, \mathfrak {f}}\right\} }, \quad \mathbf {A}^{\left\{ {\mathfrak {f}, \mathfrak {s}}\right\} }= \mathbf {A}^{\left\{ {\mathfrak {s}, \mathfrak {s}}\right\} }, \quad \mathbf {b}^{\left\{ {\mathfrak {f}}\right\} }= \mathbf {b}^{\left\{ {\mathfrak {s}}\right\} }, \end{aligned}$$
(C.3)

we can achieve \(s^{\left\{ {\mathfrak {f}}\right\} }\) = \(s^{\left\{ {\mathfrak {s}}\right\} }\) and \(Y^{\left\{ {\mathfrak {f}}\right\} }_j = Y^{\left\{ {\mathfrak {s}}\right\} }_j\). In fact, this GARK method degenerates into a special class of partitioned Runge–Kutta schemes known to preserve linear invariants [17]. The multirate methods in [7], for example, satisfy Eq. (C.3).

Rights and permissions

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Roberts, S., Loffeld, J., Sarshar, A. et al. Implicit Multirate GARK Methods. J Sci Comput 87, 4 (2021). https://doi.org/10.1007/s10915-020-01400-z

Download citation

  • Received:

  • Revised:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1007/s10915-020-01400-z

Keywords

Mathematics Subject Classification

Search

Navigation