Integrators

Create a vector of S solution matrices of type DT to store the solution of S internal stages for a problem with DxM dimensions.

Create a vector of S solution vectors of type DT to store the solution of S internal stages for a problem with D dimensions.

Create a vector of (S,M+1) solution vectors of type DT to store the solution of S internal stages and M random processes for a problem with D dimensions.

Create a vector of S+1 solution vectors of type DT to store the solution of S internal stages and the solution of the previous timestep for a problem with D dimensions.

Explicit Runge-Kutta Method

ERK(tableau)

Explicit Runge-Kutta integrator cache.

Explicit Runge-Kutta method solving the system

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} . \end{aligned}\]

Implicit Runge-Kutta Method

IRK(tableau)

Implicit Runge-Kutta integrator cache.

Fields

• x: nonlinear solver solution vector
• q̄: solution at previous timestep
• Q: internal stages of solution
• V: internal stages of vector field
• Y: vector field of internal stages
• J: Jacobi matrices for all internal stages

Implicit Runge-Kutta method solving the explicit system of equations for ODEs

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} . \end{aligned}\]

or the implicit systems of equations for IODEs and LODEs,

\[\begin{aligned} P_{n,i} &= \vartheta (Q_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , \\ F_{n,i} &= f (Q_{n,i}, V_{n,i}) , & P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} . \end{aligned}\]

If implicit_update is set to true, the update is computed by solving

\[\vartheta(q_{n+1}) = \vartheta(q_{n}) + h \sum \limits_{i=1}^{s} b_{i} \, f (Q_{n,j}, V_{n,j}) ,\]

otherwise it is computed explicitly by

\[q_{n+1} = q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} .\]

Usually we are interested in Lagrangian systems, where

\[\begin{aligned} P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , \end{aligned}\]

and tableaus satisfying the symplecticity conditions

\[\begin{aligned} b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , & \bar{b}_i &= b_i . \end{aligned}\]

Implicit Runge-Kutta integrator cache.

Fields

• q̄: solution at previous timestep
• p̄: momentum at previous timestep
• Q: internal stages of solution
• V: internal stages of vector field
• Θ: internal stages of one-form $\vartheta$
• F: internal stages of force field

Diagonally Implicit Runge-Kutta Method

DIRK(tableau)

Diagonally implicit Runge-Kutta integrator cache.

Fields

• x: nonlinear solver solution vector
• Q: internal stages of solution q
• V: internal stages of vector field v = q̇
• Y: summed vector field of internal stages Q

Diagonally implicit Runge-Kutta integrator solving the system

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} . \end{aligned}\]

Explicit Partitioned Runge-Kutta Method

EPRK(tableau)

Explicit Runge-Kutta integrator cache.

Explicit partitioned Runge-Kutta method solving the system

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}, P_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\ F_{n,i} &= f (t_i, Q_{n,i}, P_{n,i}) , & P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , & p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} . \end{aligned}\]

Usually we are interested in Hamiltonian systems, where

\[\begin{aligned} V_{n,i} &= \dfrac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & F_{n,i} &= - \dfrac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , \end{aligned}\]

and tableaus satisfying the symplecticity conditions

\[\begin{aligned} b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , & \bar{b}_i &= b_i . \end{aligned}\]

Implicit Partitioned Runge-Kutta Method

IPRK(tableau)

Implicit partitioned Runge-Kutta integrator cache.

Fields

• x: nonlinear solver solution vector
• q̄: solution at previous timestep
• p̄: momentum at previous timestep
• Q: internal stages of solution q
• P: internal stages of momentum p
• V: internal stages of vector field v = q̇
• F: internal stages of vector field f = ṗ
• Y: summed vector field of internal stages Q
• Z: summed vector field of internal stages P

Implicit partitioned Runge-Kutta method solving the system

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}, P_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\ F_{n,i} &= f (t_i, Q_{n,i}, P_{n,i}) , & P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , & p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} , \end{aligned}\]

Usually we are interested in Hamiltonian systems, where

\[\begin{aligned} V_{n,i} &= \dfrac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & F_{n,i} &= - \dfrac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , \end{aligned}\]

and tableaus satisfying the symplecticity conditions

\[\begin{aligned} b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , & \bar{b}_i &= b_i . \end{aligned}\]

For implicit equations like IODEs and LODEs this method is solving the system

\[\begin{aligned} P_{n,i} &= \vartheta (Q_{n,i}, V_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\ F_{n,i} &= f (Q_{n,i}, V_{n,i}) , & P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , & p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} , \end{aligned}\]

Usually we are interested in Lagrangian systems, where

\[\begin{aligned} P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , \end{aligned}\]

and tableaus satisfying the symplecticity conditions

\[\begin{aligned} b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , & \bar{b}_i &= b_i . \end{aligned}\]

Fully implicit Runge-Kutta method with TableauBackwardEuler.

Reference:

Augustin-Louis Cauchy.
Équations différentielles ordinaires. Cours inédit (fragment). Douzième leçon.
Ed. Christian Gilain, Etudes Vivantes, 1981.
Page 102, Equation (5), Θ=1.

Diagonally implicit Runge-Kutta method with TableauCrankNicolson.

Reference:

J. Crank and P. Nicolson.
A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type.
Mathematical Proceedings of the Cambridge Philosophical Society, Volume 43, Issue 1, Pages 50-67, 1947.
doi: 10.1017/S0305004100023197

Diagonally implicit Runge-Kutta method with TableauCrouzeix.

Reference:

M.Crouzeix.
Sur L'approximation des équations différentielles opérationelles linéaires par des méthodes de Runge-Kutta.
Thesis. Université de Paris, 1975.

Explicit Runge-Kutta method with TableauExplicitEuler.

Reference:

Leonhard Euler.
Institutiones calculi differentialis cum eius vsu in analysi finitorum ac doctrina serierum.
Imp. Acad. Imper. Scient. Petropolitanae, Opera Omnia, Vol.X, [I.6], 1755.
In: Opera Omnia, 1st Series, Volume 11, Institutiones Calculi Integralis. Teubner, Leipzig, Pages 424-434, 1913.
Sectio secunda. Caput VII. De integratione aequationum differentialium per approximationem. Problema 85.

Explicit Runge-Kutta method with TableauExplicitMidpoint.

Reference:

Carl Runge.
Über die numerische Auflösung von Differentialgleichungen.
Mathematische Annalen, Volume 46, Pages 167-178, 1895.
doi: 10.1007/BF01446807.
Equation (2)

Explicit Runge-Kutta method with TableauForwardEuler.

Reference:

Leonhard Euler.
Institutiones calculi differentialis cum eius vsu in analysi finitorum ac doctrina serierum.
Imp. Acad. Imper. Scient. Petropolitanae, Opera Omnia, Vol.X, [I.6], 1755.
In: Opera Omnia, 1st Series, Volume 11, Institutiones Calculi Integralis. Teubner, Leipzig, Pages 424-434, 1913.
Sectio secunda. Caput VII. De integratione aequationum differentialium per approximationem. Problema 85.

Fully implicit Runge-Kutta method with TableauGauss.

References:

John C. Butcher.
Implicit Runge-Kutta processes.
Mathematics of Computation, Volume 18, Pages 50-64, 1964.
doi: 10.1090/S0025-5718-1964-0159424-9.

John C. Butcher.
Gauss Methods. 
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_115.

Explicit Runge-Kutta method with TableauHeun2.

Reference:

Karl Heun.
Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen.
Zeitschrift für Mathematik und Physik, Volume 45, Pages 23-38, 1900.
Algorithm II.

Explicit Runge-Kutta method with TableauHeun3.

Reference:

Karl Heun.
Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen.
Zeitschrift für Mathematik und Physik, Volume 45, Pages 23-38, 1900.
Algorithm VI.

Fully implicit Runge-Kutta method with TableauIRK3.

Reference:

Ernst Hairer and Gerhard Wanner.
Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems.
Springer, 1996.
Section IV.5, W-transformation.
The s=2, γ=1/2 member of the family of algebraically stable methods of order 2s-1,
obtained from the Gauss method by X = Xₛ + γ eₛ eₛᵀ.

Fully implicit Runge-Kutta method with TableauImplicitEuler.

Reference:

Augustin-Louis Cauchy.
Équations différentielles ordinaires. Cours inédit (fragment). Douzième leçon.
Ed. Christian Gilain, Etudes Vivantes, 1981.
Page 102, Equation (5), Θ=1.

Fully implicit Runge-Kutta method with TableauImplicitMidpoint.

Reference:

Augustin-Louis Cauchy.
Équations différentielles ordinaires. Cours inédit (fragment). Douzième leçon.
Ed. Christian Gilain, Etudes Vivantes, 1981.
Page 102, Equation (5), Θ=1/2.

Diagonally implicit Runge-Kutta method with TableauKraaijevangerSpijker.

Reference:

J. F. B. M. Kraaijevanger and M. N. Spijker.
Algebraic stability and error propagation in Runge-Kutta methods.
Applied Numerical Mathematics, Volume 5, Issues 1-2, Pages 71-87, 1989.
doi: 10.1016/0168-9274(89)90025-1

Explicit Runge-Kutta method with TableauKutta3.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 440

Runge-Kutta method with TableauLobattoIII.

References:

John C. Butcher.
Integration processes based on Radau quadrature formulas
Mathematics of Computation, Volume 18, Pages 233-244, 1964.
doi: 10.1090/S0025-5718-1964-0165693-1.

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123.

Runge-Kutta method with TableauLobattoIIIA.

References:

Byron Leonard Ehle
On Padé approximations to the exponential function and a-stable methods for the numerical solution of initial value problems.
Research Report CSRR 2010, Dept. AACS, University of Waterloo, 1969.

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123

Partitioned Runge-Kutta method with TableauLobattoIIIA for $q$ and TableauLobattoIIIB for $p$.

Partitioned Runge-Kutta method with TableauLobattoIIIA for $q$ and TableauLobattoIIIĀ for $p$.

Runge-Kutta method with TableauLobattoIIIB.

References:

Byron Leonard Ehle.
On Padé approximations to the exponential function and a-stable methods for the numerical solution of initial value problems.
Research Report CSRR 2010, Dept. AACS, University of Waterloo, 1969.

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123.

Partitioned Runge-Kutta method with TableauLobattoIIIB for $q$ and TableauLobattoIIIA for $p$.

Partitioned Runge-Kutta method with TableauLobattoIIIB for $q$ and TableauLobattoIIIB̄ for $p$.

Runge-Kutta method with TableauLobattoIIIC.

References:

F. H. Chipman.
A-stable Runge-Kutta processes.
BIT, Volume 11, Pages 384-388, 1971.
doi: 10.1007/BF01939406.

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123.

Partitioned Runge-Kutta method with TableauLobattoIIIC for $q$ and TableauLobattoIIIC̄ for $p$.

Partitioned Runge-Kutta method with TableauLobattoIIIC̄ for $q$ and TableauLobattoIIIC for $p$.

Runge-Kutta method with TableauLobattoIIID.

References:

R.P.K. Chan.
On symmetric Runge-Kutta methods of high order.
Computing, Volume 45, Pages 301-309, 1990.
doi: 10.1007/BF02238798

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123.

Partitioned Runge-Kutta method with TableauLobattoIIID for $q$ and TableauLobattoIIID̄ for $p$.

Runge-Kutta method with TableauLobattoIIIE.

References:

R.P.K. Chan.
On symmetric Runge-Kutta methods of high order.
Computing, Volume 45, Pages 301-309, 1990.
doi: 10.1007/BF02238798

Laurent O. Jay.
Lobatto Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_123.

Partitioned Runge-Kutta method with TableauLobattoIIIE for $q$ and TableauLobattoIIIĒ for $p$.

Runge-Kutta method with TableauLobattoIIIF.

References:

Wang Fangzong and Liao Xiaobing.
A Class of Lobatto Methods of Order 2s.
Journal of Applied Mathematics, Volume 46, Pages 6-10, 2016.

Partitioned Runge-Kutta method with TableauLobattoIIIF for $q$ and TableauLobattoIIIF̄ for $p$.

Runge-Kutta method with TableauLobattoIIIF̄.

Partitioned Runge-Kutta method with TableauLobattoIIIF̄ for $q$ and TableauLobattoIIIF for $p$.

Runge-Kutta method with TableauLobattoIIIG.

Partitioned Runge-Kutta method with TableauLobattoIIIG for $q$ and TableauLobattoIIIḠ for $p$.

Partitioned Runge-Kutta method TableauGauss for both $q$ and $p$.

References:

John C. Butcher.
Implicit Runge-Kutta processes.
Mathematics of Computation, Volume 18, Pages 50-64, 1964.
doi: 10.1090/S0025-5718-1964-0159424-9.

John C. Butcher.
Gauss Methods. 
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_115.

Diagonally implicit Runge-Kutta method with TableauQinZhang.

Reference:

M.-Z. Qin and M.-Q. Zhang.
Symplectic Runge-Kutta algorithms for Hamilton systems.
Journal of Computational Mathematics, Supplementary Issue, Pages 205-215, 1992.

Explicit Runge-Kutta method with TableauRK21.

Reference:

Karl Heun.
Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen.
Zeitschrift für Mathematik und Physik, Volume 45, Pages 23-38, 1900.
Algorithm II.

Explicit Runge-Kutta method with TableauRK22.

Reference:

Carl Runge
Über die numerische Auflösung von Differentialgleichungen.
Mathematische Annalen, Volume 46, Pages 167-178, 1895.
doi: 10.1007/BF01446807.
Equation (3)

Explicit Runge-Kutta method with TableauRK31.

Reference:

John C. Butcher
Numerical Methods for Ordinary Differential Equations. Wiley, 2016.
Page 99

Explicit Runge-Kutta method with TableauRK32.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 440

Explicit Runge-Kutta method with TableauRK4.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 443

Explicit Runge-Kutta method with TableauRK41.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 443

Explicit Runge-Kutta method with TableauRK416.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 443

Explicit Runge-Kutta method with TableauRK42.

Reference:

John C. Butcher
Numerical Methods for Ordinary Differential Equations. Wiley, 2016.
Page 102

Explicit Runge-Kutta method with TableauRK438.

Reference:

Wilhelm Kutta
Beitrag zur Näherungsweisen Integration totaler Differentialgleichungen
Zeitschrift für Mathematik und Physik, Volume 46, Pages 435-453, 1901.
Page 441

Explicit Runge-Kutta method with TableauRK5.

Reference:

John C. Butcher
Numerical Methods for Ordinary Differential Equations. Wiley, 2016.
Page 103

Runge-Kutta method with TableauRadauIA.

References:

Byron Leonard Ehle
On Padé approximations to the exponential function and a-stable methods for the numerical solution of initial value problems.
Research Report CSRR 2010, Dept. AACS, University of Waterloo, 1969.

Runge-Kutta method with TableauRadauIB.

Reference:

Sun Geng
Construction of high order symplectic Runge-Kutta methods
Journal of Computational Mathematics, Volume 11, Pages 250-260, 1993.

Runge-Kutta method with TableauRadauIIA.

References:

Byron Leonard Ehle
On Padé approximations to the exponential function and a-stable methods for the numerical solution of initial value problems.
Research Report CSRR 2010, Dept. AACS, University of Waterloo, 1969.

Owe Axelsson.
A class of A-stable methods.
BIT, Volume 9, Pages 185-199, 1969.
doi: 10.1007/BF01946812.

Ernst Hairer and Gerhard Wanner.
Radau Methods.
In: Engquist B. (eds). Encyclopedia of Applied and Computational Mathematics. Springer, Berlin, Heidelberg. 2015.
doi: 10.1007/978-3-540-70529-1_139.

Runge-Kutta method with TableauRadauIIB.

Reference:

Sun Geng
Construction of high order symplectic Runge-Kutta methods
Journal of Computational Mathematics, Volume 11, Pages 250-260, 1993.

Explicit Runge-Kutta method with TableauRalston2.

Reference:

Anthony Ralston.
Runge-Kutta Methods with Minimum Error Bounds.
Mathematics of Computation, Volume 16, Pages 431-437, 1962.
doi: 10.1090/S0025-5718-1962-0150954-0.
Equation (3.5)

Explicit Runge-Kutta method with TableauRalston3.

Reference:

Anthony Ralston.
Runge-Kutta Methods with Minimum Error Bounds.
Mathematics of Computation, Volume 16, Pages 431-437, 1962.
doi: 10.1090/S0025-5718-1962-0150954-0.
Equation (4.10)

Explicit Runge-Kutta method with TableauRunge2.

Reference:

Carl Runge
Über die numerische Auflösung von Differentialgleichungen.
Mathematische Annalen, Volume 46, Pages 167-178, 1895.
doi: 10.1007/BF01446807.
Equation (3)

Fully implicit Runge-Kutta method with TableauSRK3.

Reference:

Shan Zhao and Guo-Wei Wei.
A unified discontinuous Galerkin framework for time integration.
Mathematical Methods in the Applied Sciences, Volume 37, Issue 7, Pages 1042-1071, 2014.
doi: 10.1002/mma.2863.

Explicit Runge-Kutta method with TableauSSPRK2.

Reference:

Chi-Wang Shu, Stanley Osher.
Efficient implementation of essentially non-oscillatory shock-capturing schemes.
Journal of Computational Physics, Volume 77, Issue 2, Pages 439-471, 1988.
doi: 10.1016/0021-9991(88)90177-5.
Equation (2.16)

Explicit Runge-Kutta method with TableauSSPRK3.

Reference:

Chi-Wang Shu, Stanley Osher.
Efficient implementation of essentially non-oscillatory shock-capturing schemes.
Journal of Computational Physics, Volume 77, Issue 2, Pages 439-471, 1988.
doi: 10.1016/0021-9991(88)90177-5.
Equation (2.18)

SymplecticEulerA

Symplectic Euler method using explicit Euler for q and implicit Euler for p.

SymplecticEulerB

Symplectic Euler method using implicit Euler for q and explicit Euler for p.

Partitioned Runge-Kutta Method

PRK(tableau)

Returns an explicit or implicit partitioned Runge-Kutta method depending on the tableau.

Runge-Kutta Method

RK(tableau)

Returns an explicit, implicit or diagonally implicit Runge-Kutta method depending on the tableau.

Holds the coefficients of a projected Gauss-Legendre Runge-Kutta method.

The construction follows the W-transformation: with the Legendre basis $\hat{P}_{j}$ orthonormal on $[0,1]$ evaluated at the $s$ Gauss-Legendre nodes $c_i$,

\[P_{ij} = \hat{P}_{j-1} (c_i) , \qquad Q = P^{-1} ,\]

the Gauss tableau itself is recovered as $a = P X Q$ from the tridiagonal

\[X_{11} = \tfrac{1}{2} , \qquad X_{i+1,i} = + \xi_i , \qquad X_{i,i+1} = - \xi_i , \qquad \xi_i = \frac{1}{2 \sqrt{4 i^2 - 1}} ,\]

and a skew perturbation $A = P W Q$ is built from the skew $W$ with the single nonzero pair $W_{s,s-1} = +1$, $W_{s-1,s} = -1$. This yields the one-parameter family of tableaus

\[a(\lambda) = a + \lambda A ,\]

used by PGLRK (with $\lambda$ chosen to conserve energy) and by VPRKpTableau (with $\lambda$ chosen to enforce the Dirac constraint).

Two structural properties of $A$ matter, and both are asserted in the tests:

• $B A$ is skew for $B = \mathrm{diag}(b)$, because $P^T B P = I$. Hence $b_i A_{ij} + b_j A_{ji} = 0$ and $a(\lambda)$ satisfies the symplecticity condition for every $\lambda$ and every choice of nonzero slot in $W$.
• $b^T A = 0$ and $A \mathbb{1} = 0$, so $a(\lambda)$ retains the quadrature and row-sum conditions — and therefore the order — but this holds only for $s \geq 3$. At $s = 2$ the pair $W_{2,1} / W_{1,2}$ coincides with the order-determining $\xi_1$ entries of $X$, and the perturbation destroys consistency. The $(2,1)$ slot is the only offending slot, at every $s$.

For that reason $s \geq 3$ is required. A two-stage method would remain symplectic but lose consistency, i.e. it would be a symplectic method of reduced order.

Reference

Luigi Brugnano, Felice Iavernaro and Donato Trigiante. Energy- and quadratic invariants-preserving integrators based upon Gauss collocation formulae. SIAM Journal on Numerical Analysis 50(6), pp. 2897-2916, 2012. doi: 10.1137/110856617.

getTableauPGLRK(coeff::CoefficientsPGLRK, λ, a::AbstractMatrix)
getTableauPGLRK(coeff::CoefficientsPGLRK, λ)

Evaluate the one-parameter family of tableaus $a(\lambda) = a + \lambda A$ of a CoefficientsPGLRK. The three-argument form writes into a in place and accepts any AbstractMatrix destination, so that λ may be a dual number during automatic differentiation.

Formal Lagrangian Runge-Kutta Method.

Applies a Runge-Kutta method to the formal Lagrangian of a noncanonical Hamiltonian system. Given an LODE or IODE whose vector field $v$ is determined by

\[\dot{q} = v(q) = \bar{\Omega}^{-T} (q) \, \nabla H (q) , \qquad \bar{\Omega}_{ij} (q) = \vartheta_{i,j} (q) - \vartheta_{j,i} (q) ,\]

the formal Lagrangian introduces an adjoint variable $p$ and reads

\[\bar{L} (q, \dot{q}, p, \dot{p}) = p^{T} \dot{q} - p^{T} v(q) + \vartheta^{T} (q) \, v(q) - H(q) ,\]

where the last two (momentum-modified) terms are chosen so that $p = \vartheta(q)$ solves the adjoint equations. Variation with respect to $p$ returns the physical equations of motion, and variation with respect to $q$ gives the adjoint equations

\[\dot{p} = \dfrac{\partial \bar{L}}{\partial q} = f(q, v) + \left( \dfrac{\partial v}{\partial q} \right)^{T} \big( \vartheta(q) - p \big) , \qquad f(q, v) = \big( \nabla \vartheta (q) \big)^{T} v - \nabla H (q) .\]

Discretisation

The position is advanced by an ordinary implicit Runge-Kutta method applied to $\dot{q} = v(q)$,

\[\begin{aligned} V_{n,i} &= v (t_i, Q_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \end{aligned}\]

and the adjoint variable by the corresponding Runge-Kutta method applied to the adjoint equations with $\bar{a}_{ij} = a_{ij}$, $\bar{b}_{i} = b_{i}$,

\[\begin{aligned} P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, \dot{P}_{n,j} , & \dot{P}_{n,i} &= F_{n,i} + J_{n,i}^{T} \big( \Theta_{n,i} - P_{n,i} \big) , & p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} \, \dot{P}_{n,i} , \end{aligned}\]

where $\Theta_{n,i} = \vartheta (Q_{n,i})$, $F_{n,i} = f (Q_{n,i}, V_{n,i})$ and $J_{n,i} = \partial v / \partial q (Q_{n,i})$. Because the adjoint equations are linear in $P$, these stage equations are not solved by Newton iteration but by one linear solve of size $d s \times d s$,

\[\big( I + h \, (a \otimes J^{T}) \big) \, P = \mathbb{1} \otimes p_{n} + h \, (a \otimes I) \, \big( F + J^{T} \Theta \big) ,\]

performed after the position stages have converged. The manuscript notes that the adjoint variables need not be solved for in practice, which is why this is a post-hoc computation rather than part of the nonlinear residual.

Requirements and caveats

• The tableau must be Gauss-Legendre (or at least have an invertible $a$); the elimination of the discrete adjoint momenta relies on it.
• The problem's $\bar{v}$ must depend on $q$ only. For example GeometricProblems.HarmonicOscillator.lodeproblem has $\bar{v}(t,q,p) = p_1$, for which $\partial \bar{v} / \partial q$ vanishes identically and the method is meaningless; use degenerate_lodeproblem or GeometricProblems.LotkaVolterra2d.lodeproblem instead.
• Symplecticity is not established. The reference states explicitly that variational Runge-Kutta methods are not symplectic with respect to the noncanonical symplectic form, promises "pseudo-symplecticity" and then proves nothing; conjugate symplecticity is recorded as an open conjecture in one draft and as an unproven theorem in another. There is likewise no order theorem — the order is inherited from the underlying tableau and is verified numerically in test/verification/flrk_convergence_tests.jl. Accordingly issymplectic returns missing, and the energy is expected to drift slowly rather than be conserved.

Reference

Michael Kraus. Variational Integrators for Noncanonical Hamiltonian Systems. Working manuscript.

Constructors

FLRK(tableau::Tableau)
FLRK(method::RKMethod)

Formal Lagrangian Runge-Kutta integrator cache.

Fields

• x: nonlinear solver solution vector
• q̃, p̃, ṽ, f̃: temporaries for the initial guess
• Q, V, Y: internal stages of the position, its vector field and Y = h A V
• Θ: one-form $\vartheta$ at the internal stages
• P: adjoint momentum at the internal stages
• F: force $f = (\nabla \vartheta)^{T} v - \nabla H$ at the internal stages
• G: $J^{T} (\Theta - P)$ at the internal stages
• Ṗ: total adjoint force $F + G$
• J: Jacobians $\partial v / \partial q$ at the internal stages
• A, δP: system matrix and right-hand side of the linear adjoint solve
• jac: automatic-differentiation Jacobian of $\bar{v}$

Solve the linear stage equations of the adjoint (momentum) variable and assemble the total adjoint force $\dot{P} = F + J^{T} (\Theta - P)$.

Runs on the DT cache only, after the position stages have converged.

Projected Gauss-Legendre Runge-Kutta method: an energy-preserving Gauss collocation method for ODEs carrying an invariant h.

The method uses the one-parameter family of tableaus of CoefficientsPGLRK,

\[a(\lambda) = a + \lambda A , \qquad A = P W Q ,\]

where $a$ is the $s$-stage Gauss tableau in W-transformed form and $W$ is skew. The stages

\[\begin{aligned} Y_{n,i} &= \sum \limits_{j=1}^{s} a_{ij} (\lambda) \, V_{n,j} , & Q_{n,i} &= q_{n} + h \, Y_{n,i} , & V_{n,i} &= v (t_{n,i}, Q_{n,i}) , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \end{aligned}\]

are solved for each trial $\lambda$, and $\lambda$ is then determined by a bisection on the scalar energy-conservation condition

\[H (q_{n+1} (\lambda)) = H (q_{0}) .\]

Because $B A$ is skew for $B = \mathrm{diag}(b)$, the perturbed tableau satisfies the symplecticity condition $b_{i} a_{ij} + b_{j} a_{ji} = b_{i} b_{j}$ for every fixed $\lambda$, and because $b^{T} A = 0$ and $A \mathbb{1} = 0$ it also retains the quadrature and row-sum conditions, hence the order $2s$. Both identities are asserted in test/methods/pglrk_coefficients_tests.jl.

That does not make the method symplectic. $\lambda$ is determined from $q_{n}$ at every step, so the step map is a member of the $\lambda$-family composed with a state-dependent choice of parameter, and its Jacobian carries an additional rank-one term $(\partial \Psi / \partial \lambda) (\partial \lambda / \partial q)^{T}$. Exact energy conservation and symplecticity are in any case mutually exclusive for a general Hamiltonian system unless the method reproduces the exact flow (Ge & Marsden, 1988). Measured on a nonlinear pendulum at $\Delta t = 0.8$, the defect $\vert J^{T} \Omega J - \Omega \vert$ tends to $3.1 \times 10^{-7}$ as the finite-difference step is refined, against $7.3 \times 10^{-12}$ for Gauss(3) — the latter being finite-difference noise rather than a defect. Accordingly issymplectic returns missing.

The reference energy is that of the initial condition of the problem, not of the previous step, so the absolute energy error stays bounded instead of performing a random walk.

Note that $\lambda$ is sought in the bracket $[-\lambda_{\max}, +\lambda_{\max}]$ with $\lambda_{\max} = 2 \cdot 10^{-s}$ by default. If the energy residual does not change sign across that bracket there is no root, and the method falls back to $\lambda = 0$, i.e. to the plain Gauss method, rather than applying a large perturbation; the number of steps on which that happens is counted in the cache's nfallback.

Reference

Luigi Brugnano, Felice Iavernaro and Donato Trigiante. Energy- and quadratic invariants-preserving integrators based upon Gauss collocation formulae. SIAM Journal on Numerical Analysis 50(6), pp. 2897-2916, 2012. doi: 10.1137/110856617.

Constructors

PGLRK(coefficients::CoefficientsPGLRK; λmax = 2 / 10^s)
PGLRK(s::Int, T = Float64; λmax = 2 / 10^s)

Projected Gauss-Legendre Runge-Kutta integrator cache.

Mutable, because the projection parameter λ, the reference energy h₀ and the working tableau ā = a + λ A are updated in the course of a step.

Fields

• x: nonlinear solver solution vector, holding Y
• x₀: λ-independent starting point for the stage solves
• q: candidate solution $q_{n+1}$ for the current λ
• q̃, ṽ: temporaries for the initial guess
• Q, V, Y: internal stages
• ā: working tableau $a + \lambda A$
• λ: current value of the projection parameter
• h, h₀: current and reference value of the invariant
• nfallback: number of steps on which no root was found in $[-\lambda_{\max}, +\lambda_{\max}]$ and the method fell back to $\lambda = 0$, i.e. to the plain Gauss method

Solve the stage equations for a trial λ and return the energy residual h₀ - h(q_{n+1}(λ)).

The stage solve is restarted from the λ-independent initial guess x₀ on every call. Warm-starting from the previous trial λ instead makes the returned residual depend on the sequence of trials rather than on λ alone, which breaks the bisection: the bracket then collapses to a denormal instead of converging on the root.

Determine the projection parameter λ by bisection on the energy residual.

Falls back to λ = 0 (the plain Gauss method) when the residual already vanishes, or when it does not change sign across the bracket — in the latter case SimpleSolvers.bisection returns an endpoint, i.e. the largest admissible perturbation rather than a root. Such steps are counted in cache(int).nfallback.

Set the projection parameter λ and recompute the working tableau ā = a + λ A. Both live in the DT cache, since λ is real-valued and constant with respect to the nonlinear solver unknowns.

Discrete Euler-Lagrange Method.

We consider discrete Euler-Lagrange equations of the form

\[D_1 L_d (q_{n}, q_{n+1}) + D_2 L_d (q_{n-1}, q_{n}) = 0,\]

where $q_{n}$ approximates the solution $q(t_{n})$ and $L_d$ is any discrete Lagrangian.

Discrete Euler-Lagrange integrator cache.

Midpoint Variational Integrator in position-momentum form.

We consider a discrete Lagrangian of the form

\[L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) ,\]

where $q_{n}$ approximates the solution $q(t_{n})$. The Euler-Lagrange equations are computed as:

\[\begin{aligned} p_{n } &= - D_{1} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) , \\ p_{n+1} &= \hphantom{-} D_{2} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, \frac{q_{n+1} - q_{n}}{h} \bigg) . \end{aligned}\]

The first equation can be solved implicitly for $q_{n+1}$ given $(q_{n}, p_{n})$. The second equation can be used to explicitly compute $p_{n+1}$.

Trapezoidal Variational Integrator in position-momentum form.

We consider a discrete Lagrangian of the form

\[L_d (q_{n}, q_{n+1}) = \frac{h}{2} \bigg[ L \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + L \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) \bigg] ,\]

where $q_{n}$ approximates the solution $q(t_{n})$. The Euler-Lagrange equations are computed as:

\[\begin{aligned} p_{n } &= - D_{1} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) , \\ p_{n+1} &= \hphantom{-} D_{2} L_{d} (q_{n}, q_{n+1}) = \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n}, \frac{q_{n+1} - q_{n}}{h} \bigg) + \frac{1}{2} \frac{\partial L}{\partial v} \bigg( q_{n+1}, \frac{q_{n+1} - q_{n}}{h} \bigg) . \end{aligned}\]

The first equation can be solved implicitly for $q_{n+1}$ given $(q_{n}, p_{n})$. The second equation can be used to explicitly compute $p_{n+1}$.

Variational Partitioned Runge-Kutta Method

\[\begin{aligned} P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\ F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} - d_i \lambda , & p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} , \\ && 0 &= \sum \limits_{i=1}^{s} d_i V_i , && \end{aligned}\]

satisfying the symplecticity conditions

\[\begin{aligned} b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , & \bar{b}_i &= b_i . \end{aligned}\]

Constructors

VPRK(tableau::PartitionedTableau, d=nothing)
VPRK(tableau::Tableau, args...; kwargs...) = VPRK(PartitionedTableau(tableau), args...; kwargs...)
VPRK(tableau1::Tableau, tableau2::Tableau, args...; kwargs...) = VPRK(PartitionedTableau(tableau1, tableau2), args...; kwargs...)

Variational Partitioned Runge-Kutta Method that uses Gauss-Legendre coefficients for both $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIII for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIA for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIA for the coefficients $a_{ij}$ and TableauLobattoIIIB for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIA for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIB for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIB for the coefficients $a_{ij}$ and TableauLobattoIIIA for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIB for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIC for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIC for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIC̄ for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIID for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIID for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIE for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIE for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIF for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIF for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIF̄ for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIF̄ for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIG for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauLobattoIIIG for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauRadauIIA for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses TableauRadauIIB for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta Method that uses the symmetric TableauSRK3 coefficients with 3 stages for both $a_{ij}$ and $\bar{a}_{ij}$.

Variational Partitioned Runge-Kutta method with projection in the tableau, for degenerate Lagrangian systems. EXPERIMENTAL.

Where the other projected VPRK methods append a projection step to an unmodified Runge-Kutta method, this one perturbs the tableau itself. It uses the family of CoefficientsPGLRK,

\[a(\lambda) = a + A(\lambda) , \qquad A(\lambda) = P \, W(\lambda) \, Q ,\]

where $a$ is the $s$-stage Gauss tableau and $W(\lambda)$ is skew with the $d$ free parameters $\lambda$ placed in its high-index corner,

\[W_{s-k+1,\,s-k} = + \lambda_{k} , \qquad W_{s-k,\,s-k+1} = - \lambda_{k} , \qquad k = 1, \dots, d .\]

The stage equations are those of a VPRK method with the perturbed tableau,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \big( a_{ij} + A_{ij} (\lambda) \big) \, V_{n,j} , & P_{n,i} &= \vartheta (Q_{n,i}, V_{n,i}) , \\ Z_{n,i} &= \sum \limits_{j=1}^{s} \big( a_{ij} + A_{ij} (\lambda) \big) \, F_{n,j} , & P_{n,i} &= p_{n} + h \, Z_{n,i} , \end{aligned}\]

and $\lambda$ is fixed by requiring that the Dirac constraint hold at the end of the step,

\[\vartheta (q_{n+1}) - p_{n+1} = 0 , \qquad q_{n+1} = q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \qquad p_{n+1} = p_{n} + h \sum \limits_{i=1}^{s} b_{i} \, F_{n,i} .\]

That is $d$ scalar conditions for the $d$ components of $\lambda$, so the stage equations and the constraint are solved together as one nonlinear system of size $d (s+1)$.

Because $B A(\lambda)$ is skew for $B = \mathrm{diag}(b)$, the perturbed tableau satisfies the symplecticity condition for every $\lambda$.

Stage-count requirement

The $k$-th multiplier occupies the skew slot $(s-k+1, s-k)$, so $d$ multipliers need $s \geq d + 1$ merely to fit into $W$ — this is asserted. Full order additionally requires $s \geq d + 2$, because the slot $(2,1)$ is the one pair that coincides with the order-determining $\xi_1$ entries of $X$ and destroys the conditions $b^{T} A = 0$ and $A \mathbb{1} = 0$. At $s = d + 1$ the method therefore remains symplectic but loses order. For a two-dimensional problem that means $s \geq 4$ for full order.

Status

This method has no written-up reference. The nearest documented relative, sec:glrk of Notes on Projected VIs for Degenerate Lagrangian Systems, places the multiplier in the stage equations rather than in the tableau, and its symplecticity computation is truncated; the Internal Projection section of Projected Variational Integrators for Degenerate Lagrangian Systems that would cover this construction is an empty stub. The tableau family itself is that of Brugnano, Iavernaro & Trigiante (see CoefficientsPGLRK), with the Dirac constraint substituted for their energy constraint. Accordingly the order is measured rather than claimed, and issymplectic reflects only the algebraic property above.

Constructors

VPRKpTableau(coefficients::CoefficientsPGLRK)
VPRKpTableau(s::Int, T = Float64)

Cache of VPRKpTableau.

W, T and A are ST-typed because $\lambda$ is one of the nonlinear solver unknowns, so the perturbed tableau carries dual numbers during the Jacobian computation.

Fields

• x: nonlinear solver solution vector, holding (V₁, …, V_S, λ)
• λ: the Dirac-constraint multipliers
• q, p, θ: solution and one-form at the end of the step
• z: permanent zeros, for the unused velocity slot of ϑ
• q̃, p̃, ṽ, f̃: temporaries for the initial guess
• Q, P, V, F, Y, Z: internal stages
• W, T, A: skew generator, workspace and the perturbation $A = P W Q$

Recompute the tableau perturbation A = P W(λ) Q for the current multipliers.

The k-th multiplier is placed in the skew slot (s-k+1, s-k), i.e. as far into the high-index corner of W as possible, which is where the perturbation leaves the order conditions intact.

Symplectic Euler-A Degenerate Variational Integrator.

Symplectic Euler-B Degenerate Variational Integrator.

Midpoint Degenerate Variational Integrator.

Trapezoidal Degenerate Variational Integrator.

Degenerate variational integrator cache.

Fields

• q: internal stages of solution
• v: internal stages of vector field
• Θ: implicit function evaluated on solution
• f: vector field of implicit function

Degenerate Variational Runge-Kutta (DVRK) method for degenerate Lagrangian systems of the form

\[L (q, \dot{q}) = \vartheta (q) \cdot \dot{q} - H (q) , \qquad \vartheta_{\mu} = 0 \quad \text{for} \quad \mu = d/2+1, \, ..., \, d ,\]

that is, Lagrangians linear in the velocities whose symplectic potential has $d/2$ identically vanishing components. The internal stages read

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , & P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) = \vartheta (Q_{n,i}) , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, F_{n,j} , & F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , \end{aligned}\]

and the update rules are

\[\begin{aligned} q^{\mu}_{n+1} &= q^{\mu}_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V^{\mu}_{n,i} , & p^{\mu}_{n+1} &= p^{\mu}_{n} + h \sum \limits_{i=1}^{s} b_{i} \, F^{\mu}_{n,i} , && \mu = 1, \, ..., \, d/2 , \\ && p^{\mu}_{n+1} &= \vartheta^{\mu} (q_{n+1}) , && \mu = 1, \, ..., \, d . \end{aligned}\]

Note that the quadrature updates are imposed only on the first $d/2$ components. The remaining components of $q_{n+1}$ carry no quadrature update and are determined implicitly by the constraint $p_{n+1} = \vartheta (q_{n+1})$, which is also what defines $p_{n+1}$ on output.

The method preserves the noncanonical symplectic two-form $\omega = \mathrm{d} \vartheta$ provided that

• the symplectic potential is given in a gauge in which $\vartheta_{\mu} = 0$ for $\mu > d/2$, as assumed above,
• the $d/2 \times d/2$ matrix $\partial \vartheta_{\mu} / \partial q^{\nu}$, for $\mu \le d/2 < \nu$, is invertible,
• the coefficient matrix $A = (a_{ij})$ is invertible,
• the coefficients satisfy the symplecticity condition $b_{i} a_{ij} + b_{j} a_{ji} = b_{i} b_{j}$, and
• the momentum is initialised consistently with $p_{0} = \vartheta (q_{0})$.

The first two conditions are independent: the second says that the constraint $p_{n+1} = \vartheta (q_{n+1})$ can be solved for the components of $q_{n+1}$ that carry no quadrature update, and it is also what makes $\omega$ nondegenerate, but it does not by itself imply that the remaining components of $\vartheta$ vanish.

The Gauss-Legendre tableaus satisfy the two conditions on the coefficients for any number of stages, and DVRK(Gauss(s)) attains the full order $2s$ on Lagrangians of the above form. Applied to a Lagrangian outside this class — for instance the same system written in a gauge in which no component of $\vartheta$ vanishes, which leaves the dynamics unchanged — the method remains convergent but the order drops to $s$.

A Degenerate Variational Runge-Kutta method is instantiated by either passing a Runge-Kutta tableau or a Runge-Kutta method:

DVRK(tableau::Tableau; check_conditions = true)
DVRK(method::RKMethod; check_conditions = true)

The constructor checks the two tableau conditions and warns if either is violated; pass check_conditions = false to suppress the warnings.

Degenerate Variational Runge-Kutta integrator cache.

Fields

• q: initial guess of solution
• v: initial guess of vector field
• θ: initial guess of symplectic potential
• f: initial guess of force field
• Q: internal stages of solution
• V: internal stages of vector field
• Θ: implicit function of internal stages
• F: vector field of implicit function

check_dvrk_initial_conditions(problem)

Warn if the initial momentum of problem is not consistent with the primary constraint p₀ = ϑ(t₀, q₀, v₀). Symplecticity of DVRK is conditional on this consistency, and an inconsistent p₀ silently degrades the method.

Continuous Galerkin Variational Integrator.

• b: weights of the quadrature rule
• c: nodes of the quadrature rule
• x: nodes of the basis
• m: mass matrix
• a: derivative matrix
• r₀: reconstruction coefficients at the beginning of the interval
• r₁: reconstruction coefficients at the end of the interval

Cache shared by all five DGVI variants.

The field set is the union of what the variants actually read; the legacy per-variant caches allocated roughly twice as many vectors, most of which were computed and never used. Φ, Λ, Θ, Θ̄, G and Ḡ hold the per-node values of the jump path integral and are empty (F = 0) for every variant but DGVIPI.

Abstract supertype of the Discontinuous Galerkin Variational Integrators.

DGVIs discretise the action of a fully degenerate Lagrangian

\[L(q, \dot{q}) = \vartheta (q) \cdot \dot{q} - H(q)\]

by a piecewise-polynomial trajectory that is not required to be continuous across the interval boundaries, plus a numerical flux that integrates the Lagrangian across the resulting jumps. Within each interval,

\[q_h(t) \vert_{(t_{n}, t_{n+1})} = \sum \limits_{i=1}^{S} x_{n,i} \, \varphi_{i} \bigg( \frac{t - t_{n}}{t_{n+1} - t_{n}} \bigg) ,\]

and with a quadrature rule of $R$ nodes $c_i$ and weights $b_i$ the mass and derivative matrices and the boundary reconstruction vectors are

\[m_{ij} = \varphi_j (c_i) , \qquad a_{ij} = \varphi_j' (c_i) , \qquad r^{+}_{j} = \varphi_j (0) , \qquad r^{-}_{j} = \varphi_j (1) ,\]

so that

\[Q_{n,i} = m_{ij} x_{n,j} , \qquad V_{n,i} = \frac{a_{ij}}{h} x_{n,j} , \qquad q_{n}^{+} = r^{+}_{j} x_{n,j} , \qquad q_{n+1}^{-} = r^{-}_{j} x_{n,j} .\]

All variants share the interior contribution to the discrete Euler-Lagrange equations,

\[0 = \sum \limits_{i=1}^{R} b_i \big[ h \, m_{ij} \, F_{n,i} + a_{ij} \, P_{n,i} \big] + (\text{flux terms}) ,\]

with $P_{n,i} = \vartheta (Q_{n,i})$ and $F_{n,i} = f (Q_{n,i}, V_{n,i})$, and differ only in the flux and in the equation that closes the step. See the individual method docstrings and docs/src/integrators/dgvi.md.

The five variants are DGVI, DGVIPI, DGVIP0, DGVIP1 and DGVIEXP.

State that a DGVI variant carries from one step to the next.

DGVI is a genuine $(q_n, p_n)$ map and uses none of it; the other four propagate one or two of the jump values, which are not part of the SolutionStep. This lives in the DT cache only — it is constant with respect to the nonlinear solver unknowns and must never be dual-typed — and is therefore always reached through dgvi_state rather than by naming the field.

The state is seeded once, on the first step, from the initial condition. A DGVI GeometricIntegrator is therefore single-run: restarting one from different initial conditions would carry over the jump values of the previous run. Build a fresh integrator instead.

Per-variant precondition on the (basis, quadrature) pair, checked by the generated constructors below. The default imposes nothing; DGVIP0 overrides it.

Compute the one-form and the forces at the quadrature nodes.

Compute the solution at the quadrature nodes and the two boundary reconstructions qₙ⁺ = r⁺·X and qₙ₊₁⁻ = r⁻·X.

Compute the velocities at the quadrature nodes, V = a X / h.

Compute the shared DGVI coefficient block (b, c, x, m, a, r⁻, r⁺) from a basis and a quadrature rule.

Compute the jump path-integral coefficient block (β, γ, μ⁻, μ⁺, α⁻, α⁺, ρ⁻, ρ⁺) of a Discontinuity.

μ∓ and α∓ are the "mass" and "derivative" vectors of the connecting path, μ^∓_i = φ^∓(γ_i) and α^∓_i = dφ^∓/dτ(γ_i), and ρ∓ = φ^∓(1/2) reconstruct the nodal value q_n = Φ(1/2; q_n^-, q_n^+).

The jump values a DGVI variant carries from one step to the next, as a DGVIState.

Always reach the state through this function. DGVIState is a field of the ST-parameterised DGVICache, so cache(int, ST).state also exists — but as a separate, dual-typed object, so writing to it would be silently lost. Going through dgvi_state keeps that impossible to get wrong.

Seed the carried-over jump values on the first step. The generic version starts from a continuous trajectory, qₙ⁻ = qₙ⁺ = qₙ, and ϑ(qₙ⁻) = pₙ; variants override it.

Number of degrees of freedom the step adds on top of the S basis coefficients.

Number of quadrature nodes of the jump path integral; zero for every variant except DGVIPI.

Interior contribution to the discrete Euler-Lagrange equations, shared by all variants: b[i] = Σⱼ bⱼ (h mⱼᵢ Fⱼ + aⱼᵢ Pⱼ). The variant's residual! adds its flux terms on top and writes the trailing block(s).

Carry the jump values of this step over to the next one. Overridden per variant.

Discontinuous Galerkin Variational Integrator with a trapezoidal numerical flux.

The discrete action augments the interior quadrature with the flux evaluated as a trapezoidal rule between the nodal value $q_n$ and the two one-sided limits,

\[\mathcal{A}_d = h \sum \limits_{n} \bigg[ \sum \limits_{i=1}^{R} b_i \big[ \vartheta (Q_{n,i}) \cdot V_{n,i} - H(Q_{n,i}) \big] + \frac{\vartheta (q_n) + \vartheta (q_n^+)}{2} \cdot (q_n^+ - q_n) + \frac{\vartheta (q_{n+1}^-) + \vartheta (q_{n+1})}{2} \cdot (q_{n+1} - q_{n+1}^-) \bigg] ,\]

whose stationarity gives, for all $j$,

\[0 = \sum \limits_{i=1}^{R} b_i \big[ h \, m_{ij} \, F_{n,i} + a_{ij} \, P_{n,i} \big] + r^{+}_{j} \, \frac{\vartheta (q_n) + \vartheta (q_n^+)}{2} - r^{-}_{j} \, \frac{\vartheta (q_{n+1}^-) + \vartheta (q_{n+1})}{2} + \frac{r^{+}_{j}}{2} \, \nabla \vartheta^{T} (q_n^+) \cdot (q_n^+ - q_n) + \frac{r^{-}_{j}}{2} \, \nabla \vartheta^{T} (q_{n+1}^-) \cdot (q_{n+1} - q_{n+1}^-) ,\]

together with the variation with respect to the nodal value,

\[\vartheta (q_n^+) = \vartheta (q_n^-) + \nabla \vartheta^{T} (q_n) \cdot (q_n^+ - q_n^-) .\]

Introducing

\[p_{n} = \vartheta (q_{n}^-) + \nabla \vartheta^{T} (q_{n}) \cdot (q_{n} - q_{n}^-)\]

turns the latter into

\[\vartheta (q_n^+) = p_n + \nabla \vartheta^{T} (q_n) \cdot (q_n^+ - q_n) ,\]

so the scheme is a genuine map $(q_n, p_n) \mapsto (q_{n+1}, p_{n+1})$ and needs no state beyond the solution itself — unlike the other four variants. The step is closed by

\[p_{n+1} = \vartheta (q_{n+1}^-) + \nabla \vartheta^{T} (q_{n+1}) \cdot (q_{n+1} - q_{n+1}^-) .\]

Initial conditions $q_0$ and $p_0 = \vartheta(q_0)$ must be prescribed.

The unknowns are the $S$ basis coefficients plus $q_{n+1}$.

Constructor

DGVI(basis::Basis, quadrature::QuadratureRule)

Discontinuous Galerkin Variational Integrator with a path-integral numerical flux.

This is the general variant, and the one derived in docs/src/integrators/dgvi.md. The non-conservative product of the one-form and the jump is discretised by integrating along a path $\Phi(\tau; q^-, q^+)$ connecting the two one-sided limits,

\[\int \limits_0^1 \vartheta \big( \Phi(\tau; q^-, q^+) \big) \cdot \frac{d \Phi}{d\tau} \, d\tau \approx \sum \limits_{i=1}^{F} \beta_i \, \vartheta \big( \Phi (\gamma_i) \big) \cdot \Phi' (\gamma_i) ,\]

with a quadrature rule of $F$ nodes $\gamma_i$ and weights $\beta_i$. Writing the path as $\Phi(\tau) = q^- \varphi^-(\tau) + q^+ \varphi^+(\tau)$ and setting

\[\mu^{\mp}_i = \varphi^{\mp} (\gamma_i) , \qquad \alpha^{\mp}_i = \frac{d\varphi^{\mp}}{d\tau} (\gamma_i) ,\]

the path values and their derivatives are

\[\Phi_i = \mu^-_i q^- + \mu^+_i q^+ , \qquad \Lambda_i = \alpha^-_i q^- + \alpha^+_i q^+ ,\]

and the discrete Euler-Lagrange equations read

\[0 = \sum \limits_{i=1}^{R} b_i \big[ h \, m_{ij} \, F_{n,i} + a_{ij} \, P_{n,i} \big] + \sum \limits_{i=1}^{F} \beta_i \Big[ r^{+}_{j} \alpha^{+}_i \, \vartheta (\Phi_i) + r^{-}_{j} \alpha^{-}_i \, \vartheta (\bar{\Phi}_i) + r^{+}_{j} \mu^{+}_i \, \nabla \vartheta^{T} (\Phi_i) \cdot \Lambda_i + r^{-}_{j} \mu^{-}_i \, \nabla \vartheta^{T} (\bar{\Phi}_i) \cdot \bar{\Lambda}_i \Big] ,\]

where barred quantities are those of the jump at $t_{n+1}$. Note that the full path is integrated at both boundaries, but only the $\partial / \partial q_n^+$ terms are taken at $t_n$ and only the $\partial / \partial q_{n+1}^-$ terms at $t_{n+1}$, so there is no double counting.

The step is closed by reconstructing the nodal value at the midpoint of the path,

\[q_{n} = \Phi (1/2; q_n^-, q_n^+) = \rho^- q_n^- + \rho^+ q_n^+ ,\]

and $q_n^-$ is carried over from the previous step. With Discontinuity(PathIntegralLinear(), LobattoLegendreQuadrature(2)) this reduces to the trapezoidal flux on the full jump with $\rho^- = \rho^+ = 1/2$.

Constructor

DGVIPI(basis::Basis, quadrature::QuadratureRule, jump::Discontinuity)

Seed the carried-over jump values. q₀⁻ = q₀ and q₀⁺ follows from the midpoint reconstruction q₀ = ρ⁻ q₀⁻ + ρ⁺ q₀⁺.

Discontinuous Galerkin Variational Integrator with a projection at the initial time of each interval. EXPERIMENTAL.

The interior contribution and the trapezoidal flux are those of DGVI, but the step is closed by a projection condition imposed at $t_n$ rather than by the definition of $p_n$,

\[\vartheta (q_n^+) - \vartheta (q_n^-) - \nabla \vartheta^{T} (q_n) \cdot (q_n^+ - q_n^-) = 0 ,\]

in which the one-sided one-forms are reconstructed from the stage values of the one-form,

\[\vartheta (q_n^+) = \sum \limits_{i} r^{+}_{i} P_{n,i} , \qquad \vartheta (q_{n+1}^-) = \sum \limits_{i} r^{-}_{i} P_{n,i} .\]

Consequently $\vartheta (q_n^-)$ and $q_n^-$ are carried over from the previous step.

Constructor

DGVIP0(basis::Basis, quadrature::QuadratureRule)

Discontinuous Galerkin Variational Integrator with a projection at the final time of each interval. EXPERIMENTAL.

The interior contribution and the trapezoidal flux are those of DGVI. The step is closed by two conditions: continuity of the carried-over right limit,

\[q_{n}^{+} - r^{+} \cdot x_{n} = 0 ,\]

and a projection imposed at $t_{n+1}$,

\[\vartheta (q_{n+1}^{+}) - \vartheta (q_{n+1}^{-}) - \nabla \vartheta^{T} (q_{n+1}) \cdot (q_{n+1}^{+} - q_{n+1}^{-}) = 0 .\]

The unknowns are therefore the $S$ basis coefficients, $q_{n+1}$ and $q_{n+1}^{+}$, and $q_{n}^{+}$ is carried over from the previous step.

Constructor

DGVIP1(basis::Basis, quadrature::QuadratureRule)

Discontinuous Galerkin Variational Integrator with a midpoint-average numerical flux. EXPERIMENTAL.

Instead of the trapezoidal flux of DGVI, the one-form is evaluated once, at the average of the two one-sided limits,

\[\langle q \rangle_{n} = \tfrac{1}{2} \big( q_{n}^{-} + q_{n}^{+} \big) , \qquad [\![ q ]\!]_{n} = q_{n}^{+} - q_{n}^{-} ,\]

giving the flux $\vartheta ( \langle q \rangle_{n} ) \cdot [\![ q ]\!]_{n}$ — the jump/average form of the manuscript with averaging parameter $\alpha = 1/2$. The discrete equations become

\[0 = \sum \limits_{i=1}^{R} b_i \big[ h \, m_{ij} \, F_{n,i} + a_{ij} \, P_{n,i} \big] + r^{+}_{j} \, \vartheta ( \langle q \rangle_{n} ) - r^{-}_{j} \, \vartheta ( \langle q \rangle_{n+1} ) + \tfrac{1}{2} r^{+}_{j} \, \nabla \vartheta^{T} ( \langle q \rangle_{n} ) \cdot [\![ q ]\!]_{n} + \tfrac{1}{2} r^{-}_{j} \, \nabla \vartheta^{T} ( \langle q \rangle_{n+1} ) \cdot [\![ q ]\!]_{n+1} ,\]

closed by identifying the nodal value with the average,

\[q_{n+1} = \langle q \rangle_{n+1} ,\]

and by continuity with the carried-over right limit, $q_{n}^{+} - r^{+} \cdot x_{n} = 0$. Both $q_{n}^{-}$ and $q_{n}^{+}$ are propagated from step to step.

Constructor

DGVIEXP(basis::Basis, quadrature::QuadratureRule)

Hamilton-Pontryagin Integrator using midpoint quadrature.

We consider a discrete Lagrangian of the form

\[L_d (q_{n}, q_{n+1}) = h \, L \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) ,\]

where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. The discrete Hamilton-Pontryagin action reads

\[A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_{d}^{\alpha,\beta} (q_{n}, q_{n+1}, v_{n+1/2}) + \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] ,\]

where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is

\[\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} .\]

In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$.

The equations of motion, that are solved by this integrator, is computed as:

\[\begin{aligned} 0 &= \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) \\ &+ h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) . \end{aligned}\]

Upon defining the momentum

\[p_{n} = h \, D_2 \phi_h (q_{n-1}, q_{n}; a_{n-1,n}) \cdot p_{n-1/2} + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n-1} + q_{n}}{2}, v_{n-1/2} \bigg) ,\]

we can rewrite the equations of motion as

\[\begin{aligned} 0 &= p_{n} + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ p_{n+1/2} &= \frac{\partial L}{\partial v} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) , \\ p_{n+1} &= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + \frac{h}{2} \, \frac{\partial L}{\partial q} \bigg( \frac{q_{n} + q_{n+1}}{2}, v_{n+1/2} \bigg) . \end{aligned}\]

Given $(q_{n}, p_{n})$, the first four equations are solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method. The last equation provides an explicit update for $p_{n+1}$.

Hamilton-Pontryagin Integrator using trapezoidal quadrature.

We consider a discrete Lagrangian of the form

\[L_d (q_{n}, q_{n+1}) = \frac{h}{2} \big[ L (q_{n}, v_{n+1/2}) + L (q_{n+1}, v_{n+1/2}) \big] ,\]

where $q_{n}$ approximates the solution $q(t_n)$ and $v_{n+1/2}$ is the velocity, which is assumed to be constant in the interval $[t_{n}, t_{n+1}]$. The discrete Hamilton-Pontryagin action reads

\[A_d [q_d] = \sum \limits_{n=0}^{N-1} \bigg[ L_d (q_{n}, q_{n+1}) + h \left< p_{n+1/2} , \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) - v_{n+1/2} \right> \bigg] ,\]

where $\phi_h$ is a map that computes the velocity $v_{n+1/2}$ as a function of $q_{n}$, $q_{n+1}$ and a set of parameters $a_{n,n+1}$. A trivial example of such a map that does not depend on any parameters $a_{n,n+1}$ is

\[\phi_h (q_{n}, q_{n+1}; a_{n,n+1}) = \frac{q_{n+1} - q_{n}}{h} .\]

In order to solve the discrete Euler-Lagrange equations, the user needs to specify the map $\phi_h$, its gradients with respect to $q_{n}$ and $q_{n+1}$, denoted by $D_1 \phi_h$ and $D_2 \phi_h$, respectively, the gradient with respect to the parameters, denoted by $D_a \phi_h$, and an initial set of parameters $a_{0}$.

The equations of motion, that are solved by this integrator, are:

\[\begin{aligned} 0 &= p_{n} + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n}, v_{n+1/2}) + h \, D_1 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ 0 &= D_a \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} , \\ v_{n+1/2} &= \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) , \\ p_{n+1/2} &= \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n}, v_{n+1/2}) + \frac{1}{2} \, \frac{\partial L}{\partial v} (q_{n+1}, v_{n+1/2}) , \\ p_{n+1} &= h \, D_2 \phi_h (q_{n}, q_{n+1}; a_{n,n+1}) \cdot p_{n+1/2} + \frac{h}{2} \, \frac{\partial L}{\partial q} (q_{n+1}, v_{n+1/2}) . \end{aligned}\]

Given $(q_{n}, p_{n})$, the first four equations are solved for $q_{n+1}$, where $v_{n+1/2}$, $p_{n+1/2}$, and $a_{n,n+1}$ are treated as internal variables similar to the internal stages of a Runge-Kutta method. The last equation provides an explicit update for $p_{n+1}$.

Implements the GeometricIntegrator interface for the exact solution of a substep in SODE.

Symmetric splitting method with general stages. See McLachlan, Quispel, 2003, Equ. (4.11).

Basic method: Lie composition

\[\begin{aligned} \varphi_{\tau}^{A} &= \varphi_{\tau}^{v_1} \circ \varphi_{\tau}^{v_2} \circ \dotsc \circ \varphi_{\tau}^{v_{r-1}} \circ \varphi_{\tau}^{v_r} \\ \varphi_{\tau}^{B} &= \varphi_{\tau}^{v_r} \circ \varphi_{\tau}^{v_{r-1}} \circ \dotsc \circ \varphi_{\tau}^{v_2} \circ \varphi_{\tau}^{v_1} \end{aligned}\]

Integrator:

\[\varphi_{\tau}^{GS} = \varphi_{a_1 \tau}^{A} \circ \varphi_{b_1 \tau}^{B} \circ \dotsc \circ \varphi_{b_1 \tau}^{B} \circ \varphi_{a_1 \tau}^{A}\]

General splitting method for vector fields with two components A and B.

Integrator:

\[\varphi_{\tau} = \varphi_{b_s \tau}^{B} \circ \varphi_{a_s \tau}^{A} \circ \dotsc \circ \varphi_{b_1 \tau}^{B} \circ \varphi_{a_1 \tau}^{A}\]

Non-symmetric splitting method. See McLachlan, Quispel, 2003, Equ. (4.10). The methods A and B are the composition of all vector fields in the SODE and its adjoint, respectively.

Basic method: Lie composition

\[\begin{aligned} \varphi_{\tau}^{A} &= \varphi_{\tau}^{v_1} \circ \varphi_{\tau}^{v_2} \circ \dotsc \circ \varphi_{\tau}^{v_{r-1}} \circ \varphi_{\tau}^{v_r} \\ \varphi_{\tau}^{B} &= \varphi_{\tau}^{v_r} \circ \varphi_{\tau}^{v_{r-1}} \circ \dotsc \circ \varphi_{\tau}^{v_2} \circ \varphi_{\tau}^{v_1} \end{aligned}\]

Integrator:

\[\varphi_{\tau}^{NS} = \varphi_{b_s \tau}^{B} \circ \varphi_{a_s \tau}^{A} \circ \dotsc \circ \varphi_{b_1 \tau}^{B} \circ \varphi_{a_1 \tau}^{A}\]

Symmetric splitting method with symmetric stages. See McLachlan, Quispel, 2003, Equ. (4.6).

Basic method: symmetric Strang composition

\[\varphi_{\tau}^{A} = \varphi_{\tau/2}^{v_1} \circ \varphi_{\tau/2}^{v_2} \circ \dotsc \circ \varphi_{\tau/2}^{v_{r-1}} \circ \varphi_{\tau/2}^{v_r} \circ \varphi_{\tau/2}^{v_r} \circ \varphi_{\tau/2}^{v_{r-1}} \circ \dotsc \circ \varphi_{\tau/2}^{v_2} \circ \varphi_{\tau/2}^{v_1}\]

Integrator:

\[\varphi_{\tau}^{SS} = \varphi_{a_1 \tau}^{A} \circ \varphi_{a_2 \tau}^{A} \circ \dotsc \circ \varphi_{a_s \tau}^{A} \circ \dotsc \circ \varphi_{a_2 \tau}^{A} \circ \varphi_{a_1 \tau}^{A}\]

Lie-Trotter Splitting A

For a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$, the splitting reads

\[\Phi_{\Delta t} = \varphi^{2}_{\Delta t} \circ \varphi^{1}_{\Delta t}\]

Reference:

H. F. Trotter.
On the product of semi-groups of operators.
Proceedings of the American Mathematical Society, Volume 10, Pages 545-551, 1959.
doi: 10.1090/S0002-9939-1959-0108732-6.

Lie-Trotter Splitting B

For a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$, the splitting reads

\[\Phi_{\Delta t} = \varphi^{1}_{\Delta t} \circ \varphi^{2}_{\Delta t}\]

Reference:

H. F. Trotter.
On the product of semi-groups of operators.
Proceedings of the American Mathematical Society, Volume 10, Pages 545-551, 1959.
doi: 10.1090/S0002-9939-1959-0108732-6.

Alias for Strang

McLachlan's 2nd order symmetric, minimum error composition method

The composition reads

\[\Phi_{\Delta t} = \varphi_{\alpha \Delta t} \circ \varphi^{*}_{(1/2 - \alpha) \Delta t} \circ \varphi_{(1/2 - \alpha) \Delta t} \circ \varphi^{*}_{\alpha \Delta t} ,\]

where the parameter $\alpha$ can be optimized, e.g., to minimize the solution error. McLachlan arrives at $\alpha = 0.1932$ as a generally useful value.

Reference:

Robert I. McLachlan.
On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
doi: 10.1137/0916010.

McLachlan's 4th order symmetric, minimum error composition method

The composition reads

\[\Phi_{\Delta t} = \varphi_{\alpha_5 \Delta t} \circ \varphi^{*}_{\beta_5 \Delta t} \circ \dotsc \circ \varphi_{\alpha_2 \Delta t} \circ \varphi^{*}_{\beta_2 \Delta t} \circ \varphi_{\alpha_1 \Delta t} \circ \varphi^{*}_{\beta_1 \Delta t} ,\]

with

\[\begin{aligned} \beta_1 &= \alpha_5 = \frac{14 - \sqrt{19}}{108} , & \alpha_1 &= \beta_5 = \frac{146 + 5 \sqrt{19}}{540} , & & \\ \beta_2 &= \alpha_4 = \frac{- 23 - 20 \sqrt{19}}{270} , & \alpha_2 &= \beta_4 = \frac{-2 + 10 \sqrt{19}}{135} , & \beta_3 &= \alpha_3 = \frac{1}{5} . \end{aligned}\]

The coefficients are optimised to provide an integrator with minimal solution error.

Reference:

Robert I. McLachlan.
On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
doi: 10.1137/0916010.

Strang Splitting

For a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$, the splitting reads

\[\Phi_{\Delta t} = \varphi^{1}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t / 2} \circ \varphi^{1}_{\Delta t / 2}\]

For vector fields with two components, this is not the most efficient implementation. For such cases StrangA or StrangB should be used instead.

References:

Gilbert Strang.
On the construction and comparison of difference schemes.
SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
doi: 10.1137/0705041.

Gurij Ivanovich Marchuk.
Some applications of splitting-up methods to the solution of mathematical physics problems.
Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
doi: 10.21136/AM.1968.103142.

Strang Splitting A for a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$.

The splitting reads

\[\Phi_{\Delta t} = \varphi^{1}_{\Delta t / 2} \circ \varphi^{2}_{\Delta t} \circ \varphi^{1}_{\Delta t / 2}\]

References:

Gilbert Strang.
On the construction and comparison of difference schemes.
SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
doi: 10.1137/0705041.

Gurij Ivanovich Marchuk.
Some applications of splitting-up methods to the solution of mathematical physics problems.
Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
doi: 10.21136/AM.1968.103142.

Strang Splitting B for a vector field $\dot{x} = f_1 (t,x) + f_2 (t,x)$

The splitting reads

\[\Phi_{\Delta t} = \varphi^{2}_{\Delta t / 2} \circ \varphi^{1}_{\Delta t} \circ \varphi^{2}_{\Delta t / 2}\]

References:

Gilbert Strang.
On the construction and comparison of difference schemes.
SIAM Journal on Numerical Analysis, Volume 5, Pages 506-517, 1968.
doi: 10.1137/0705041.

Gurij Ivanovich Marchuk.
Some applications of splitting-up methods to the solution of mathematical physics problems.
Aplikace Matematiky, Volume 13, Pages 103-132, 1968.
doi: 10.21136/AM.1968.103142.

Suzuki's 4th order "fractal" composition method

The composition reads

\[\Phi_{\Delta t} = \varphi_{\gamma_5 \Delta t} \circ \varphi_{\gamma_4 \Delta t} \circ \varphi_{\gamma_3 \Delta t} \circ \varphi_{\gamma_2 \Delta t} \circ \varphi_{\gamma_1 \Delta t}\]

with

\[\gamma_1 = \gamma_2 = \gamma_4 = \gamma_5 = \frac{1}{4 - 4^{1/(p+1)}} , \qquad \gamma_3 = - \frac{4^{1/(p+1)}}{4 - 4^{1/(p+1)}} .\]

Reference:

Masuo Suzuki
Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations.
Physics Letters A, Volume 146, Pages 319-323, 1990.
doi: 10.1016/0375-9601(90)90962-N

4th order "Triple Jump" composition method.

The composition reads

\[\Phi_{\Delta t} = \varphi_{\gamma_3 \Delta t} \circ \varphi_{\gamma_2 \Delta t} \circ \varphi_{\gamma_1 \Delta t}\]

with

\[\gamma_1 = \gamma_3 = \frac{1}{2 - 2^{1/(p+1)}} , \qquad \gamma_2 = - \frac{2^{1/(p+1)}}{2 - 2^{1/(p+1)}} .\]

References:

Michael Creutz and Andreas Gocksch.
Higher-order hybrid Monte Carlo algorithms.
Physical Review Letters, Volume 63, Pages 9-12, 1989.
doi: 10.1103/PhysRevLett.63.9.

Etienne Forest.
Canonical integrators as tracking codes (or how to integrate perturbation theory with tracking).
AIP Conference Proceedings, Volume 184, Pages 1106-1136, 1989.
doi: 10.1063/1.38062.

Masuo Suzuki
Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations.
Physics Letters A, Volume 146, Pages 319-323, 1990.
doi: 10.1016/0375-9601(90)90962-N

Haruo Yoshida.
Construction of higher order symplectic integrators.
Physics Letters A, Volume 150, Pages 262-268, 1990.
doi: 10.1016/0375-9601(90)90092-3

Yoshida's 6th order symmetric composition method.

The method is the symmetric composition of $m = 7$ symmetric (2nd order) steps,

\[\Phi_{\Delta t} = S_{w_1 \Delta t} \circ S_{w_2 \Delta t} \circ S_{w_3 \Delta t} \circ S_{w_4 \Delta t} \circ S_{w_3 \Delta t} \circ S_{w_2 \Delta t} \circ S_{w_1 \Delta t} ,\]

with $w_1 = 0.78451361047755726381949763$, $w_2 = 0.23557321335935813368479318$, $w_3 = -1.17767998417887100694641568$, and $w_4 = 1 - 2 (w_1 + w_2 + w_3)$ (Yoshida's "solution A").

References:

Haruo Yoshida.
Construction of higher order symplectic integrators.
Physics Letters A, Volume 150, Pages 262-268, 1990.
doi: 10.1016/0375-9601(90)90092-3

Robert I. McLachlan.
On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
doi: 10.1137/0916010.

Yoshida's 8th order symmetric composition method.

The method is the symmetric composition of $m = 15$ symmetric (2nd order) steps with weights $w_1, \dotsc, w_8, w_7, \dotsc, w_1$ (Yoshida's "solution D"),

\[w_1 = 0.74167036435061295345 , \quad w_2 = -0.40910082580003159400 , \quad w_3 = 0.19075471029623837995 , \quad w_4 = -0.57386247111608226464 ,\]

\[w_5 = 0.29906418130365592384 , \quad w_6 = 0.33462491824529818378 , \quad w_7 = 0.31529309239676659663 , \quad w_8 = 1 - 2 \sum_{i=1}^{7} w_i .\]

References:

Haruo Yoshida.
Construction of higher order symplectic integrators.
Physics Letters A, Volume 150, Pages 262-268, 1990.
doi: 10.1016/0375-9601(90)90092-3

Robert I. McLachlan.
On the Numerical Integration of Ordinary Differential Equations by Symmetric Composition Methods
SIAM Journal on Scientific Computing, Volume 16, Pages 151-168, 1995.
doi: 10.1137/0916010.

Splitting integrator for the solution of initial value problems

\[\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,\]

whose vector field $v$ is given as a sum of vector fields

\[v (t) = v_1 (t) + ... + v_r (t) .\]

Splitting has two constructors:

Splitting(f::Vector{Int}, c::Vector)
Splitting(method::AbstractSplittingMethod, problem::SODEProblem)

The vectors f and c define the actual splitting method: f is a vector of indices of the flows in the split equation (that is the exact solution which can be obtained by solutions(problem)) to be solved and c is a vector of the same size as f that contains the coefficients for each splitting step, i.e., the resulting integrator has the form

\[\varphi_{\tau} = \phi_{c[s] \tau}^{v_{f[s]}} \circ \dotsc \circ \phi_{c[2] \tau}^{v_{f[2]}} \circ \phi_{c[1] \tau}^{v_{f[1]}} .\]

In the second constructor, these vectors are constructed from the splitting method and the problem.

Splitting integrator cache.

Composition integrator for the solution of initial value problems

\[\dot{q} (t) = v(t, q(t)) , \qquad q(t_{0}) = q_{0} ,\]

whose vector field $v$ is given as a sum of vector fields

\[v (t) = v_1 (t) + ... + v_r (t) .\]

Composition has one main constructor:

Composition(methods, splitting)

The first argument is tuple of methods that are used to solve each substep, and the second argument is a splitting method. A second convenience constructor uses the ExactSolution for all steps:

Composition(splitting) = Composition(ExactSolution(), splitting)

This constructs a composition method that is equivalent to a plain Splitting. In order to include exact solutions in the composition, the ExactSolution method implements the general integrator interface.