SPARK

Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems.

Holds the coefficients of an additive Runge-Kutta method.

Holds the coefficients of an additive Runge-Kutta method.

Holds the multiplier Runge-Kutta coefficients.

Holds the coefficients of a projective Runge-Kutta method.

Holds the coefficients of a SPARK method.

Cache of a Specialised Partitioned Additive Runge-Kutta integrator.

Fields

• n: time step number
• t: time of current time step
• t̄: time of previous time step
• q: current solution of q
• q̄: previous solution of q
• p: current solution of p
• p̄: previous solution of p
• v: vector field of q
• v̄: vector field of q̄
• f: vector field of p
• f̄: vector field of p̄
• q̃: initial guess of q
• p̃: initial guess of p
• ṽ: initial guess of v
• f̃: initial guess of f
• s̃: holds shift due to periodicity of solution
• Q: internal stages of q
• P: internal stages of p
• V: internal stages of v
• F: internal stages of f
• Y: vector field of internal stages of q
• Z: vector field of internal stages of p

Partitioned Additive Runge-Kutta integrator for Hamiltonian systems subject to Dirac constraints EXPERIMENTAL.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r , \end{aligned}\]

with definitions

\[\begin{aligned} V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= - \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for Hamiltonian systems EXPERIMENTAL.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \end{aligned}\]

with definitions

\[\begin{aligned} V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= - \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for Hamiltonian systems EXPERIMENTAL.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \end{aligned}\]

with definitions

\[\begin{aligned} V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= - \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for Hamiltonian systems subject to Dirac constraints with projection on secondary constraint EXPERIMENTAL.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \end{aligned}\]

with definitions

\[\begin{aligned} V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= - \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) . \end{aligned}\]

Specialised Lobatto Runge-Kutta integrator for degenerate variational systems with projection on the secondary constraint.

SLRK targets degenerate Lagrangians of the form

\[L (q, v) = \vartheta (q) \cdot v - H (q) ,\]

whose fibre derivative $p = \vartheta(q)$ is not invertible for $v$, so that the equations of motion take the form of an index-two DAE with the primary constraint $\phi (q,p) = p - \vartheta (q) = 0$ and the secondary constraint $\psi = \dot{p} - \dot{q} \cdot \nabla \vartheta (q) = \tfrac{d}{dt} \phi = 0$.

Unlike the SPARK methods of Jay, which average the primary constraints over the internal stages, SLRK imposes the primary constraint at every stage and averages the secondary constraints instead. This is what makes the primary constraint hold to round-off along the whole trajectory.

In contrast to VSPARKsecondary there is only a single set of $s$ stages: the projective stages coincide with the internal stages, so $\sigma = s$ and the $\Lambda$ multipliers live on the same nodes as the velocities $V$.

One step $(q_n, p_n) \mapsto (q_{n+1}, p_{n+1})$ solves

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{s} a^2_{ij} \Lambda_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} G_{n,j} , & i &= 1, ..., s , \\ 0 &= \Phi_{n,i} - \frac{d_i}{\bar{b}^1_i} \, \mu , & i &= 1, ..., s , \\ 0 &= \sum \limits_{j=1}^{s} \omega_{ij} \, \Psi_{n,j} + \omega_{i,s+1} \, \phi (q_{n+1}, p_{n+1}) , & i &= 1, ..., s , \\ 0 &= \sum \limits_{i=1}^{s} d_{i} \, V_{n,i} , \end{aligned}\]

with the update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} \Lambda_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}^1_{i} F_{n,i} + h \sum \limits_{i=1}^{s} \bar{b}^2_{i} G_{n,i} , \end{aligned}\]

and the definitions

\[\begin{aligned} F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) = - \nabla H (Q_{n,i}) + \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , \\ G_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot \Lambda_{n,i} , \\ \Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , \\ \Psi_{n,i} &= \psi (Q_{n,i}, V_{n,i}, P_{n,i}, F_{n,i}) = F_{n,i} - V_{n,i} \cdot \nabla \vartheta (Q_{n,i}) . \end{aligned}\]

The four coefficient tableaus are stored as q $\to a^1, b^1$, q̃ $\to a^2, b^2$, p $\to \bar{a}^1, \bar{b}^1$ and p̃ $\to \bar{a}^2, \bar{b}^2$. All SLRKLobattoIII* constructors set q̃ = q and p̃ = p, so a method is fixed by one Lobatto pair $(a, \bar{a})$.

The $\omega$ constraint

$\omega$ is the $s \times (s+1)$ matrix returned by lobatto_ω_matrix(s). The system $\omega \, [\Psi_{n,1}, \ldots, \Psi_{n,s}, \phi(q_{n+1},p_{n+1})]^T = 0$ is equivalent to imposing the $s-1$ Lobatto-IIIA-averaged secondary constraints $\sum_j a^{\mathrm{IIIA}}_{ij} \Psi_{n,j} = 0$ for $i = 2, ..., s$ together with $\phi (q_{n+1}, p_{n+1}) = 0$.

The null vector

The Lobatto stage system is rank deficient by one in the $V$-direction. The multiplier $\mu$ relaxes the primary constraint along the null vector $d$ of the Lagrange-derivative matrix (get_lobatto_nullvector), and the deficiency is removed by the extra condition $\sum_i d_i V_{n,i} = 0$. $\mu$ enters the primary-constraint equation and only that one, as in VPARK/VSPARK/SPARK. (VPRK implements the same idea with a different residual layout: it has no primary-constraint row, so it perturbs the momentum-stage equation instead — see residual_correction!.)

Requirements on the problem

Properties

• The primary constraint $\phi(q_n, p_n) = 0$ is preserved to round-off at every step (measured: $\max |\phi| \approx 2 \times 10^{-15}$ over $10^3$ steps on Lotka–Volterra), which is the point of the method.
• The order is $2s-2$ for all six SLRKLobattoIII* families, confirmed empirically at $s = 2, 3$ in test/verification/spark_convergence_tests.jl.
• Symplecticity is approximate, not exact. The underlying manuscript claims exact preservation of the noncanonical symplectic form, but its proof leaves an uncontrolled term in the multiplier block: constraint $\sum_i d_i V_{n,i} = 0$ has no counterpart for $\Lambda$, so a residual $h \sum_j b^2_j (d_j / \bar{b}^1_j) \, \mathrm{d}\mu \wedge \mathrm{d}\Lambda_{n,j}$ survives. Measured on Lotka–Volterra, the first Poincaré invariant $\oint p \, dq$ drifts secularly at $O(h^{p+1})$ per step, whereas a genuinely symplectic variational integrator holds the same (canonical) invariant to round-off out to $3 \times 10^3$ steps. See the "Fifth pass" section of the verification report, docs/src/audit.md, for the numbers.
• The defect is gauge invariant: two Lagrangians whose one-forms differ by an exact form give the same trajectory and the same drift, to round-off.

Specialised Partitioned Additive Runge-Kutta integrator for index-two DAE systems.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \end{aligned}\]

with definitions

\[\begin{aligned} P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) . \end{aligned}\]

Variational partitioned additive Runge-Kutta integrator.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r , \end{aligned}\]

with definitions

\[\begin{aligned} P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for Variational systems.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \end{aligned}\]

with definitions

\[\begin{aligned} P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for Variational systems.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\ \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\ 0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 , \\ 0 &= \sum \limits_{i=1}^{r} \tilde{d}_i \, \Lambda_{n,i} , \end{aligned}\]

with definitions

\[\begin{aligned} P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ G_{n,i} &= - \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\ \tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} . \end{aligned}\]

Specialised Partitioned Additive Runge-Kutta integrator for degenerate variational systems with projection on secondary constraint.

This integrator solves the following system of equations for the internal stages,

\[\begin{aligned} Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{\sigma} a^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., s , \\ P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} F^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{a}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., s , \\ 0 &= \Phi_{n,i} , & i &= 1, ..., s , \end{aligned}\]

the projective stages

\[\begin{aligned} \tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \alpha^1_{ij} \tilde{V}_{n,j} + h \sum \limits_{j=1}^{\sigma} \alpha^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., \sigma , \\ \tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{\alpha}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{\alpha}^2_{ij} \tilde{F}^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{\alpha}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., \sigma , \\ 0 &= \tilde{\Phi}_{n,i} , & i &= 1, ..., \sigma , \\ 0 &= \sum \limits_{j=1}^{\sigma} \omega_{ij} \tilde{\Psi}_{n,i} , & i &= 1, ..., \sigma-1 , \end{aligned}\]

and update rule

\[\begin{aligned} q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{\sigma} b^2_{i} \tilde{\Lambda}_{n,i} , \\ p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b^1_{i} F^1_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} F^2_{n,i} + h \sum \limits_{i=1}^{\sigma} b^3_{i} \tilde{F}^3_{n,i} , \\ 0 &= \phi (q_{n+1}, p_{n+1}) , \end{aligned}\]

with definitions

\[\begin{aligned} F^1_{n,i} &= \nabla H (Q_{n,i}) , & i &= 1, ..., s , \\ F^2_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , & i &= 1, ..., s , \\ \tilde{F}^3_{n,i} &= \nabla \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) \cdot \tilde{\Lambda}_{n,i} , & i &= 1, ..., \sigma , \\ \Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , & i &= 1, ..., s , \\ \tilde{\Phi}_{n,i} &= \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = \tilde{P}_{n,i} - \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \\ \tilde{\Psi}_{n,i} &= \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) = \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \end{aligned}\]

so that

\[\begin{aligned} F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ \phi (Q_{n,i}, P_{n,i}) &= P_{n,i} - \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\ \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) &= \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \frac{\partial L}{\partial v} (\tilde{Q}_{n,i}, \tilde{V}_{n,i}) , & i &= 1, ..., \sigma . \end{aligned}\]

Holds all parameters of an Specialised Partitioned Additive Runge-Kutta method for variational systems subject to constraints.

Holds the tableau of a Hamiltonian Partitioned Additive Runge-Kutta methods.

Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method.

Holds the tableau of an Variational Partitioned Additive Runge-Kutta method.

Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems.

Print additive Runge-Kutta coefficients.

Print additive Runge-Kutta coefficients.

Print multiplier Runge-Kutta coefficients.

Print projective Runge-Kutta coefficients.

Print SPARK coefficients.

SPARK tableau for Gauss-Legendre Runge-Kutta method with s stages.

SPARK tableau for Gauss-Legendre/Gauss-Lobatto-IIIA-IIIB methods.

SPARK tableau for Gauss-Legendre/Gauss-Lobatto-IIIB-IIIA methods.

Tableau for Variational Gauss-Legendre method with s stages.

Tableau for Gauss-Lobatto IIIA-IIIB-IIIC method with s stages.

Tableau for Gauss-Lobatto IIIA-IIIB-IIID method with s stages.

SPARK tableau for Gauss-Lobatto methods.

SPARK tableau for Gauss-Lobatto IIIA-IIIB method with s stages.

SPARK tableau for Gauss-Lobatto IIIB-IIIA method with s stages.

SPARK Tableau for Variational Partitioned Runge-Kutta Methods.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Legendre HPARK method with s stages.

SPARK tableau for Gauss-Lobatto IIIA-IIIB HPARK method with s stages.

SPARK tableau for Gauss-Lobatto IIIB-IIIA method with s stages.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symmetric projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symmetric projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symplectic projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symplectic projection.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Legendre method with s stages and Lobatto-IIIA-IIIB projection.

Tableau for Gauss-Legendre method with s stages and Lobatto-IIIB-IIIA projection.

Tableau for Gauss-Legendre method with s stages and midpoint projection.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Tableau for Gauss-Legendre method with s stages and midpoint projection.

Tableau for Gauss-Legendre method with s stages and symplectic projection.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. Use the same tableaus for $\tilde{a}^{1}$ and $\tilde{a}^{3}$, so that $\tilde{s} = s$, as well as

\[\begin{aligned} \begin{array}{c|cc} & \tfrac{1}{2} b^{1} \\ & \vdots \\ & \tfrac{1}{2} b^{1} \\ \hline a^{2} & \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{3} \\ & \vdots \\ & \tfrac{1}{2} b^{3} \\ \hline a^{4} & \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{1} \\ c^{1} & \vdots \\ & \tfrac{1}{2} b^{1} \\ \hline \tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty)) \, b^{1} \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{3} \\ c^{3} & \vdots \\ & \tfrac{1}{2} b^{3} \\ \hline \tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty)) \, b^{3} \\ \end{array} \end{aligned}\]

Set $\omega = [0, ..., 0, 1]$ and

\[\delta_{ij} = \begin{cases} +1 & j = i , \\ -1 & j = \tilde{s} , \\ 0 & \text{else} , \end{cases}\]

so that $\Lambda_{1} = \Lambda_{2} = ... = \Lambda_{\tilde{s}}$.

This methods is constructed to satisfy the constraint on the projective stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$. Note, however, that it violates the symplecticity conditions $b^{1}_{i} b^{4}_{j} = b^{1}_{i} a^{4}_{ij} + b^{4}_{j} \tilde{a}^{1}_{ji}$ and $b^{2}_{i} b^{3}_{j} = b^{2}_{i} \tilde{a}^{3}_{ij} + b^{3}_{j} a^{2}_{ji}$.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. For the projection, choose the Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{4}, b^{4})$ and $(\tilde{a}^{2}, b^{2})$, respectively, and choose $\tilde{a}^{1}$ and $\tilde{a}^{3}$ such that the projective stages correspond to the initial condition and the solution, i.e.,

\[\begin{aligned} \begin{array}{c|cc} 0 & 0 \\ 1 & b^{1} \\ \hline \tilde{a}^{1} & \\ \end{array} && \begin{array}{c|cc} 0 & \tfrac{1}{2} & 0 \\ 1 & \tfrac{1}{2} & 0 \\ \hline \tilde{a}^{2} & \\ \end{array} && \begin{array}{c|cc} 0 & 0 \\ 1 & b^{3} \\ \hline \tilde{a}^{3} & \\ \end{array} && \begin{array}{c|cc} 0 & 0 & 0 \\ 1 & \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{4} & \\ \end{array} \end{aligned}\]

and compute $a^{2}$ and $a^{4}$ by the symplecticity conditions, that is $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$ and $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$. Finally choose $\omega = [0, 0, 1]$ and $\delta = [-1, R_{\infty}]$, implying that $\rho = 1$. By construction, this method satisfies all symplecticity conditions, but the constraint on the projection stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, is not satisfied exactly, but only approximately, although with bounded error.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIA-IIIB projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIB-IIIA projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and midpoint projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIA-IIIB projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and Lobatto-IIIB-IIIA projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and midpoint projection.

Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symmetric projection.

This methods is the same as TableauVSPARKLobattoIIIAIIIBProjection, except for using Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{2}, b^{2})$, and $(\tilde{a}^{4}, b^{4})$ respectively, instead of the other way around.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIA-IIIB projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIB-IIIA projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and midpoint projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIA-IIIB projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and Lobatto-IIIB-IIIA projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and midpoint projection.

Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symmetric projection.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. For the projection, choose the tableau with $\tilde{s} = 1$ and $\rho = 0$, such that $\tilde{Q}_{n,1} = \tfrac{1}{2} ( q_{n} + q_{n+1})$, $\tilde{P}_{n,1} = \tfrac{1}{2} ( p_{n} + p_{n+1})$, i.e.,

\[\begin{aligned} \begin{array}{c|cc} & \tfrac{1}{2} \\ & \vdots \\ & \tfrac{1}{2} \\ \hline a^{2} & \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} \\ & \vdots \\ & \tfrac{1}{2} \\ \hline a^{4} & \\ \end{array} && \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty))\\ \end{array} && \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty))\\ \end{array} \end{aligned}\]

The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, and $\omega = [0, 1]$.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$, and set

\[\begin{aligned} \begin{array}{c|cc} & \tfrac{1}{2} b^{1} \\ & \vdots \\ & \tfrac{1}{2} b^{1} \\ \hline a^{2} & \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{4} \\ & \vdots \\ & \tfrac{1}{2} b^{4} \\ \hline a^{4} & \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{1} \\ c^{1} & \vdots \\ & \tfrac{1}{2} b^{1} \\ \hline \tilde{a}^{2} & \tfrac{1}{2} (1 + R(\infty)) \, b^{1} \\ \end{array} && \begin{array}{c|cc} & \tfrac{1}{2} b^{3} \\ c^{3} & \vdots \\ & \tfrac{1}{2} b^{3} \\ \hline \tilde{a}^{4} & \tfrac{1}{2} (1 + R(\infty)) \, b^{3} \\ \end{array} \end{aligned}\]

Note that by this definition $\tilde{s} = s$. The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the (modified) symplecticity conditions, specifically $a^{4}_{ij} = b^{3}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{1}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, where $b^{2}$ has been replaced with $b^{1}$ and $b^{4}$ with $b^{3}$, respectively. Set $\omega = [0, ..., 0, 1]$ and

\[\delta_{ij} = \begin{cases} +1 & j = i , \\ -1 & j = \tilde{s} , \\ 0 & \text{else} , \end{cases}\]

so that $\Lambda_{1} = \Lambda_{2} = ... = \Lambda_{\tilde{s}}$.

Note that this method satisfies the symplecticity conditions $b^{1}_{i} b^{4}_{j} = b^{1}_{i} a^{4}_{ij} + b^{4}_{j} \tilde{a}^{1}_{ji}$ and $b^{2}_{i} b^{3}_{j} = b^{2}_{i} \tilde{a}^{3}_{ij} + b^{3}_{j} a^{2}_{ji}$ only if $R(\infty) = 1$ due to the definitions of $b^{2}$ and $b^{4}$. Moreover, it does usually not satisfy the constraint on the projective stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, exactly, but only approximately with bounded error, thus implying a residual in the symplecticity equation even if $R(\infty) = 1$.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. For the projection, choose the Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{4}, b^{4})$ and $(\tilde{a}^{2}, b^{2})$, respectively.

The coefficients $\tilde{a}^{1}$ and $\tilde{a}^{3}$ are determined by the relations

\[\begin{aligned} \sum \limits_{j=1}^{s} \tilde{a}^{1}_{ij} (c_{j}^{1})^{k-1} &= \frac{(c_{i}^{2})^k}{k} , \qquad & \sum \limits_{j=1}^{s} \tilde{a}^{3}_{ij} (c_{j}^{3})^{k-1} &= \frac{(c_{i}^{4})^k}{k} , \qquad & i &= 1 , \, ... , \, \tilde{s} , \qquad & k &= 1 , \, ... , \, s . \end{aligned}\]

The coefficients $a^{2}$ and $a^{4}$ by the symplecticity conditions, that is $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$ and $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$. Finally choose $\omega = [0, 0, 1]$ and $\delta = [-1, R_{\infty}]$, implying that $\rho = 1$. By construction, this method satisfies all symplecticity conditions, but the constraint on the projection stages, $\phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = 0$ for $i = 1, \, ..., \, \tilde{s}$, is not satisfied exactly, but only approximately, although with bounded error.

This methods is the same as TableauVSPARKModifiedLobattoIIIAIIIBProjection, except for using Lobatto-IIIA and IIIB tableaus with $\tilde{s} = 2$ stages for $(\tilde{a}^{2}, b^{2})$, and $(\tilde{a}^{4}, b^{4})$ respectively, instead of the other way around.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. For the projection, choose the tableau with $\tilde{s} = 1$ and $\rho = 0$, such that $\tilde{Q}_{n,1} = \tfrac{1}{2} ( q_{n} + q_{n+1})$, $\tilde{P}_{n,1} = \tfrac{1}{2} ( p_{n} + p_{n+1})$, i.e.,

\[\begin{aligned} \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} b^{1} \\ \hline \tilde{a}^{1} & \\ \end{array} && \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{2} & \tfrac{1}{2} ( 1 + R (\infty) ) \\ \end{array} && \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} b^{3} \\ \hline \tilde{a}^{3} & \\ \end{array} && \begin{array}{c|cc} \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{4} & \tfrac{1}{2} ( 1 + R (\infty) ) \\ \end{array} \end{aligned}\]

The coefficients $a^{2}$ and $a^{4}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$, and $\omega = [0, 1]$.

Consider a symplectic pair of tableaus $(a^{1}, b^{1}, c^{1})$ and $(a^{3}, b^{3}, c^{3})$, i.e., satsifying $b^{1}_{i} b^{3}_{j} = b^{1}_{i} a^{3}_{ij} + b^{3}_{j} a^{1}_{ji}$, with an arbitrary number of stages $s$. For the projection, choose the tableau with $\tilde{s} = 2$ and $\rho = 1$, such that $\tilde{Q}_{n,1} = q_{n}$, $\tilde{Q}_{n,2} = q_{n+1}$, $\tilde{P}_{n,1} = p_{n}$, $\tilde{P}_{n,2} = p_{n+1}$, i.e.,

\[\begin{aligned} \begin{array}{c|cc} 0 & 0 \\ 1 & b^{1} \\ \hline \tilde{a}^{1} & \\ \end{array} && \begin{array}{c|cc} & 0 & 0 \\ & \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{2} & \tfrac{1}{2} & \tfrac{1}{2} \\ \end{array} && \begin{array}{c|cc} 0 & 0 \\ 1 & b^{3} \\ \hline \tilde{a}^{3} & \\ \end{array} && \begin{array}{c|cc} & 0 & 0 \\ & \tfrac{1}{2} & \tfrac{1}{2} \\ \hline \tilde{a}^{4} & \tfrac{1}{2} & \tfrac{1}{2} \\ \end{array} \end{aligned}\]

The coefficients $a^{2}$ and $a^{4}$ are determined by the symplecticity conditions, specifically $a^{4}_{ij} = b^{4}_{j} ( b^{1}_{i} - \tilde{a}^{1}_{ji}) / b^{1}_{i}$ and $a^{2}_{ij} = b^{2}_{j} ( b^{3}_{i} - \tilde{a}^{3}_{ji} ) / b^{3}_{i}$. Further choose $\omega = [1, 1, 0]$ and $\delta = [-1, R_{\infty}]$, so that $\tilde{\Lambda}_{n,1} = R_{\infty} \tilde{\Lambda}_{n,2}$ and

\[\tilde{P}_{n,1} - \vartheta (\tilde{Q}_{n,1}) + R_{\infty} ( \tilde{P}_{n,2} - \vartheta (\tilde{Q}_{n,2}) ) = 0 .\]

Due to the particular choice of projective stages, this is equivalent to

\[p_{n} - \vartheta (q_{n}) + R_{\infty} ( p_{n+1} - \vartheta (q_{n+1}) ) = 0 ,\]

so that the constraint $\phi(q_{n+1}, p_{n+1}) = 0$ is satisfied if $\phi(q_{n}, p_{n}) = 0$. Note that the choice of $\tilde{a}^{2}$ and $\tilde{a}^{4}$ violates the symplecticity condition $b^{2}_{i} b^{4}_{j} = b^{2}_{i} \tilde{a}^{4}_{ij} + b^{4}_{j} \tilde{a}^{2}_{ji}$.

The projective Lobatto-GLRK coefficients are implicitly given by

\[\sum \limits_{j=1}^{s} a_{ij} c_{j}^{k-1} = \frac{\bar{c}_i^k}{k} \qquad i = 1 , \, ... , \, \sigma , \; k = 1 , \, ... , \, s ,\]

where $c$ are Gauß-Legendre nodes with $s$ stages and $\bar{c}$ are Gauß-Lobatto nodes with $\sigma$ stages.

nullvectorsize(problem, method)

Number of unknowns the null-vector multiplier μ contributes to the nonlinear system: one component per degree of freedom when the tableau carries a null vector d, zero otherwise.

Every solversize method in this submodule adds this term, so solversize is the full length of the solver's unknown vector and the cache constructor needs no correction of its own. Keep it that way — the two used to be computed in different files, and nothing checked that they agreed (see the SPARK solver size testset).