GeometricIntegrator

Collects all data structures needed by an integrator:

• problem: EquationProblem to solve
• method: integration method
• cache: temprary data structures needed by method
• solver: linear or nonlinear solver needed by method
• iguess: initial guess for implicit methods
• projection: optional projection method

Constructors:

GeometricIntegrator(problem::EquationProblem, method::GeometricMethod; solver = default_solver(method), iguess = default_iguess(method), projection = default_projection(method))

Solve for time steps n with n₁ ≤ n ≤ n₂.

integrate!(solution, integrator, n₁, n₂)

Solve for all time steps n:

integrate!(solution, integrator)

internal_variables(::Integrator) = NamedTuple()

Returns a NamedTuple containing all internal variables of an integrator that shall be stored in an SolutionStep. If there is no method for a specific integrator implemented an empty NamedTuple() is returned.

Compute p(x) where p is the unique polynomial of degree length(xi), such that p(x[i]) = y[i]) for all i. Call with

aitken_neville!(x::AbstractVector, t::Real, ti::AbstractVector, xi::AbstractMatrix)

where

• x: evaluation value
• t: evaluation point
• ti: interpolation nodes
• xi: interpolation values

Euler extrapolation method with arbitrary order p.

Solves the ordinary differential equation

\[\begin{aligned} \dot{x} &= v(t, x) , & x(t_0) &= x_0 , \end{aligned}\]

for $x_1 = x(t_1)$, and is called with

extrapolate!(t₀, x₀, t₁, x₁, problem, EulerExtrapolation(s))

where

• t₀: initial time
• t₁: final time
• x₀: initial value $x_0 = x(t_0)$
• x₁: final value $x_1 = x(t_1)$
• problem: ODEProblem whose solution to extrapolate
• s: number of interpolations (order $p=s+1$)

Hermite's Interpolating Polynomials

Implements a two point Hermite inter-/extrapolation function which passes through the function and its first derivative for the interval $[0,1]$. The polynomial is determined by four constraint equations, matching the function and its derivative at the points $0$ and $1$.

Call with one of the following methods

extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, t, x, HermiteExtrapolation())
extrapolate!(t₀, x₀, ẋ₀, t₁, x₁, ẋ₁, t, x, ẋ, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, v, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, ẋ, v, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, problem, HermiteExtrapolation())
extrapolate!(t₀, x₀, t₁, x₁, t, x, ẋ, problem, HermiteExtrapolation())

where

• t₀: first sample time $t_0$
• x₀: first solution value $x_0 = x(t_0)$
• ẋ₀: first vector field value $ẋ_0 = v(t_0, x(t_0))$
• t₁: second sample time $t_1$
• x₁: second solution value $x_1 = x(t_1)$
• ẋ₁: second vector field value $ẋ_1 = v(t_1, x(t_1))$
• t: time $t$ to extrapolate
• x: extrapolated solution value $x(t)$
• ẋ: extrapolated vector field value $ẋ(t)$
• v: function to compute vector field with signature v(ẋ,t,x)
• problem: ODEProblem whose vector field to use

Derivation

The interpolation works as follows: Start by defining the 3rd degree polynomial and its derivative by

\[\begin{aligned} g(x) &= a_0 + a_1 x + a_2 x^2 + a_3 x^3 , \\ g'(x) &= a_1 + 2 a_2 x + 3 a_3 x^2 , \end{aligned}\]

and apply the constraints

\[\begin{aligned} g(0) &= f_0 & & \Rightarrow & a_0 &= f_0 , \\ g(1) &= f_1 & & \Rightarrow & a_0 + a_1 + a_2 + a_3 &= f_1 , \\ g'(0) &= f'_0 & & \Rightarrow & a_1 &= f'_0 , \\ g'(1) &= f'_1 & & \Rightarrow & a_1 + 2 a_2 + 3 a_3 &= f'_1 . \\ \end{aligned}\]

Solving for $a_0, a_1, a_2, a_3$ leads to

\[\begin{aligned} a_0 &= f_0 , & a_1 &= f'_0 , & a_2 &= - 3 f_0 + 3 f_1 - 2 f'_0 - f'_1 , & a_3 &= 2 f_0 - 2 f_1 + f'_0 + f'_1 , \end{aligned}\]

so that the polynomial $g(x)$ reads

\[g(x) = f_0 + f'_0 x + (- 3 f_0 + 3 f_1 - 2 f'_0 - f'_1) x^2 + (2 f_0 - 2 f_1 + f'_0 + f'_1) x^3 .\]

The function and derivative values can be factored out, so that $g(x)$ can be rewritten as

\[g(x) = f_0 (1 - 3 x^2 + 2 x^3) + f_1 (3 x^2 - 2 x^3) + f'_0 (x - 2 x^2 + x^3) + f'_1 (- x^2 + x^3) ,\]

or in generic form as

\[g(x) = f_0 a_0(x) + f_1 a_1(x) + f'_0 b_0(x) + f'_1 b_1(x) ,\]

with basis functions

\[\begin{aligned} a_0 (x) &= 1 - 3 x^2 + 2 x^3 , & b_0 (x) &= x - 2 x^2 + x^3 , \\ a_1 (x) &= 3 x^2 - 2 x^3 , & b_1 (x) &= - x^2 + x^3 . \end{aligned}\]

The derivative $g'(x)$ accordingly reads

\[g'(x) = f_0 a'_0(x) + f_1 a'_1(x) + f'_0 b'_0(x) + f'_1 b'_1(x) ,\]

with

\[\begin{aligned} a'_0 (x) &= - 6 x + 6 x^2 , & b'_0 (x) &= 1 - 4 x + 3 x^2 , \\ a'_1 (x) &= 6 x - 6 x^2 , & b'_1 (x) &= - 2 x + 3 x^2 . \end{aligned}\]

The basis functions $a_0$and $a_1$ are associated with the function values at $x_0$ and $x_1$, respectively, while the basis functions $b_0$ and $b_1$ are associated with the derivative values at $x_0$ and $x_1$. The basis functions satisfy the following relations,

\[\begin{aligned} a_i (x_j) &= \delta_{ij} , & b_i (x_j) &= 0 , & a'_i (x_j) &= 0 , & b'_i (x_j) &= \delta_{ij} , & i,j &= 0, 1 , \end{aligned}\]

where $\delta_{ij}$ denotes the Kronecker-delta, so that

\[\begin{aligned} g(0) &= f_0 , & g(1) &= f_1 , & g'(0) &= f'_0 , & g'(1) &= f'_1 . \end{aligned}\]

Midpoint extrapolation method with arbitrary order p.

For an ODEProblem, this solves the ordinary differential equation

\[\begin{aligned} \dot{x} &= v(t, x) , & x(t_0) &= x_0 , \end{aligned}\]

for $x_1 = x(t_1)$, and is called with

extrapolate!(t₀, x₀, t₁, x₁, ::ODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• x₀: initial value $x_0 = x(t_0)$
• t₁: final time
• x₁: final value $x_1 = x(t_1)$
• s: number of interpolations (order $p=2s+2$)

For a PODEProblem or HODEProblem, this solves the partitioned ordinary differential equation

\[\begin{aligned} \dot{q} &= v(t, q, p) , & q(t_0) &= q_0 , \\ \dot{p} &= f(t, q, p) , & p(t_0) &= p_0 , \end{aligned}\]

for $q_1 = q(t_1)$ and $p_1 = p(t_1)$m and is called with

extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::PODEProblem, MidpointExtrapolation(s))
extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::HODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• q₀: initial position $q_0 = q(t_0)$
• p₀: initial momentum $p_0 = p(t_0)$
• t₁: final time
• q₁: final position $q_1 = q(t_1)$
• p₁: final momentum $p_1 = p(t_1)$
• s: number of interpolations (order $p=2s+2$)

Similarly, for a IODEProblem or LODEProblem, this solves the explicit dynamical equation

\[\begin{aligned} \dot{q} &= v(t, q) , & q(t_0) &= q_0 , \\ \dot{p} &= f(t, q, v) , & p(t_0) &= p_0 , \end{aligned}\]

corresponding to the implicit problem, for $q_1 = q(t_1)$ and $p_1 = p(t_1)$, and is called with

extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::IODEProblem, MidpointExtrapolation(s))
extrapolate!(t₀, q₀, p₀, t₁, q₁, p₁, ::LODEProblem, MidpointExtrapolation(s))

where

• t₀: initial time
• q₀: initial position $q_0 = q(t_0)$
• p₀: initial momentum $p_0 = p(t_0)$
• t₁: final time
• q₁: final position $q_1 = q(t_1)$
• p₁: final momentum $p_1 = p(t_1)$
• s: number of interpolations (order $p=2s+2$)

Explicit Euler Method.

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.

Implicit Euler integrator cache.

Implicit Euler Method.

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.

Implicit Euler integrator cache.

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

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} + b_{j} a_{ji} &= b_{i} 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}$.

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 noncanonical symplectic equations 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}\]

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

DVRK(tableau::Tableau)
DVRK(method::RKMethod)

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

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

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

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.