Base

The window, in iterations, after which a nonlinear solve that is not descending towards f_abstol is given up on; the f_stall_window of default_options.

Without it, a solve whose residual sits on a floor above f_abstol while its iterate keeps moving normally is invisible: no convergence branch fires, max_stalls cannot see it (the step is not small), and the solve spends max_iterations — 1000 by default — on every time step.

50 is the conservative point: it gives up only on a residual that is essentially flat, which no Newton() solve on these residuals is, and it stays above SimpleSolvers.F_STALL_REPORT_MINIMUM, below which SimpleSolvers itself holds the no-progress proportion to be no evidence of anything. See SimpleSolvers.F_STALL_WINDOW for the linear rate a given window corresponds to at a given f_stall_factor.

The one case it is not conservative for is a solve that converges linearly at a rate close to that threshold — a Picard() fixed-point iteration on a stiff problem near its convergence limit, say. Every method in this package uses Newton(), but default_options is reached with whatever solvermethod the caller hands to GeometricIntegrator, so such a caller has to raise this.

Crank-Nicolson Method, also known as the trapezoidal rule.

For an ordinary differential equation $\dot{q} = v(t,q)$, that is an ODEProblem, the method reads

\[q_{n+1} = q_{n} + \frac{h}{2} \, \big[ v (t_{n}, q_{n}) + v (t_{n+1}, q_{n+1}) \big]\]

The nonlinear solver solves for the vector field at the new time step, $V = v(t_{n+1}, q_{n+1})$, while $\bar{V} = v(t_{n}, q_{n})$ is computed once per time step. The method is symmetric and second order, but not symplectic (it is conjugate to a symplectic method).

For a partitioned differential equation, that is a PODEProblem or HODEProblem,

\[\dot{q} = v (t, q, p) , \qquad \dot{p} = f (t, q, p) ,\]

the trapezoidal rule is applied to both components,

\[\begin{aligned} q_{n+1} &= q_{n} + \frac{h}{2} \, ( \bar{V} + V ) , & p_{n+1} &= p_{n} + \frac{h}{2} \, ( \bar{F} + F ) , \end{aligned}\]

where $\bar{V} = v(t_{n}, q_{n}, p_{n})$ and $\bar{F} = f(t_{n}, q_{n}, p_{n})$ are computed once per time step, as they are given by the solution at the beginning of the time step, while $V = v(t_{n+1}, q_{n+1}, p_{n+1})$ and $F = f(t_{n+1}, q_{n+1}, p_{n+1})$ are solved for. The solver solution vector holds both of them, so the nonlinear system is twice the size of the one for an ordinary differential equation.

For an implicit differential equation, that is an IODEProblem or LODEProblem,

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

the trapezoidal rule is applied to the position and to the momentum alike, which amounts to the Lobatto IIIA method with two stages,

\[\begin{aligned} q_{n+1} &= q_{n} + \frac{h}{2} \, ( \bar{V} + V ) , & p_{n+1} &= p_{n} + \frac{h}{2} \, \big[ f (t_{n}, q_{n}, \bar{V}) + f (t_{n+1}, q_{n+1}, V) \big] . \end{aligned}\]

In contrast to the explicit case, the velocity at the beginning of the time step is not given by a function evaluation but implicitly by $\vartheta (t_{n}, q_{n}, \bar{V}) = p_{n}$. Both velocities are therefore solved for simultaneously, so that the nonlinear system reads

\[\begin{aligned} 0 &= \vartheta (t_{n}, q_{n}, \bar{V}) - p_{n} , \\ 0 &= \vartheta (t_{n+1}, q_{n+1}, V) - p_{n} - \frac{h}{2} \, \big[ f (t_{n}, q_{n}, \bar{V}) + f (t_{n+1}, q_{n+1}, V) \big] , \end{aligned}\]

and is twice the size of the one for an ordinary differential equation. Whenever $\vartheta$ is regular, so that the implicit equation is equivalent to a partitioned ordinary differential equation, this is the same map as the Crank-Nicolson method applied to that equation.

Crank-Nicolson integrator cache.

Fields

• x: nonlinear solver solution vector, holding the vector field $V = v(t_{n+1}, q_{n+1})$
• q: solution at the end of the time step
• v: vector field at the end of the time step
• v̄: vector field at the beginning of the time step, $v(t_{n}, q_{n})$, which is constant throughout the nonlinear solve

Crank-Nicolson integrator cache for implicit differential equations.

In contrast to CrankNicolsonCache, nothing is constant during the solve: the velocity at the beginning of the time step is part of the nonlinear solver solution vector, so every field is computed from it and has to be accessed through cache(int, ST).

Fields

• x: nonlinear solver solution vector, holding $\bar{V}$ and $V$
• q: solution at the end of the time step
• v, θ, f: velocity, momentum map and force at the end of the time step
• v̄, θ̄, f̄: velocity, momentum map and force at the beginning of the time step

Crank-Nicolson integrator cache for partitioned differential equations.

As in CrankNicolsonCache, the vector fields at the beginning of the time step are constant during the solve, so they are read from the cache at working precision.

Fields

• x: nonlinear solver solution vector, holding $V$ and $F$
• q, p: solution at the end of the time step
• v, f: vector fields at the end of the time step
• v̄, f̄: vector fields at the beginning of the time step, $v(t_{n}, q_{n}, p_{n})$ and $f(t_{n}, q_{n}, p_{n})$, which are constant throughout the nonlinear solve

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

Explicit Euler Method.

\[q_{n+1} = q_{n} + h \, v (t_{n}, q_{n})\]

Explicit Euler integrator cache.

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

GeometricMethod is the abstract supertype for all integration methods implemented in GeometricIntegrators.

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

See NormalizedHermiteExtrapolation for a version of this extrapolation that is normalised to the interval $[0,1]$.

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}\]

Implicit Euler Method.

Implicit Euler integrator cache.

Implicit Midpoint Method.

For an ordinary differential equation $\dot{q} = v(t,q)$, that is an ODEProblem, the method reads

\[q_{n+1} = q_{n} + h \, v \bigg( t_{n} + \frac{h}{2} , \frac{q_{n} + q_{n+1}}{2} \bigg)\]

The nonlinear solver solves for the stage vector field $V = v(t_{n} + h/2, q_{n} + h V / 2)$, so that the update reads $q_{n+1} = q_{n} + h V$.

For a partitioned differential equation, that is a PODEProblem or HODEProblem,

\[\dot{q} = v (t, q, p) , \qquad \dot{p} = f (t, q, p) ,\]

the same quadrature is applied to both components. With the stage time $\tilde{t} = t_{n} + h/2$ and the midpoints $Q = q_{n} + h V / 2$ and $P = p_{n} + h F / 2$, the nonlinear solver solves

\[V = v (\tilde{t}, Q, P) , \qquad F = f (\tilde{t}, Q, P) ,\]

for the stage vector fields $V$ and $F$, and the updates read

\[q_{n+1} = q_{n} + h \, V , \qquad p_{n+1} = p_{n} + h \, F .\]

Both components are coupled through $v$ and $f$, so the solver solution vector holds both stage vector fields and the nonlinear system is twice the size of the one for an ordinary differential equation. Even for a separable Hamiltonian the method does not decouple into two explicit substeps the way the symplectic Euler methods do: it is the Gauss method with a single stage, applied to the partitioned system.

For an implicit differential equation, that is an IODEProblem or LODEProblem,

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

the same quadrature is applied to the momentum map and to the force, which amounts to the Gauss method with a single stage. With the stage time $\tilde{t} = t_{n} + h/2$ and the midpoint $Q = q_{n} + h V / 2$, the nonlinear solver solves

\[\vartheta (\tilde{t}, Q, V) = p_{n} + \frac{h}{2} \, f (\tilde{t}, Q, V)\]

for the stage velocity $V$, and the updates read

\[\begin{aligned} q_{n+1} &= q_{n} + h \, V , & p_{n+1} &= p_{n} + h \, f (\tilde{t}, Q, V) . \end{aligned}\]

Whenever $\vartheta$ is regular, so that the implicit equation is equivalent to a partitioned ordinary differential equation, this is the same map as the implicit midpoint method applied to that equation. The solver solution vector holds $V$ in both cases, so the nonlinear system has the same size as for an ordinary differential equation.

Implicit midpoint integrator cache.

Fields

• x: nonlinear solver solution vector, holding the stage vector field $V$
• q: midpoint of the time step, $Q = q_{n} + h V / 2$
• v: stage vector field $V$

Implicit midpoint integrator cache for implicit differential equations.

Fields

• x: nonlinear solver solution vector, holding the stage velocity $V$
• q: midpoint of the time step, $Q = q_{n} + h V / 2$
• v: stage velocity $V$
• θ: momentum map $\vartheta$ evaluated at the stage
• f: force $f$ evaluated at the stage

Implicit midpoint integrator cache for partitioned differential equations.

Fields

• x: nonlinear solver solution vector, holding the stage vector fields $V$ and $F$
• q: midpoint of the time step, $Q = q_{n} + h V / 2$
• p: midpoint of the time step, $P = p_{n} + h F / 2$
• v: stage vector field $V$
• f: stage vector field $F$

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)$, 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$)

Normalized Hermite's Interpolating Polynomials

Implements the same two point Hermite inter-/extrapolation as HermiteExtrapolation, but normalised to the interval $[0,1]$, so that no sample times need to be passed. Instead of the time $t$ to extrapolate to, the normalised time $c$ is passed, which corresponds to

\[t = t_1 + c \, \Delta t , \qquad \Delta t = t_1 - t_0 ,\]

where $t_0$ and $t_1$ are the times of the first and second sample. Hence $c = -1$ reproduces the first and $c = 0$ the second sample.

As the interpolation is normalised, all derivative values are with respect to the normalised time $c$, that is they are scaled by $\Delta t$ compared to the vector field values of HermiteExtrapolation.

Call with one of the following methods

extrapolate!(x₀, ẋ₀, x₁, ẋ₁, c, x, NormalizedHermiteExtrapolation())
extrapolate!(x₀, ẋ₀, x₁, ẋ₁, c, x, ẋ, NormalizedHermiteExtrapolation())

where

• x₀: first solution value $x_0 = x(t_0)$
• ẋ₀: first derivative value $ẋ_0 = \Delta t \, v(t_0, x(t_0))$
• x₁: second solution value $x_1 = x(t_1)$
• ẋ₁: second derivative value $ẋ_1 = \Delta t \, v(t_1, x(t_1))$
• c: normalised time $c$ to extrapolate, corresponding to $t = t_1 + c \, \Delta t$
• x: extrapolated solution value $x(t)$
• ẋ: extrapolated derivative value $\Delta t \, ẋ(t)$

Basis functions

Substituting $x \to 1 + c$ into the basis functions $a_i(x)$ and $b_i(x)$ derived for HermiteExtrapolation gives

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

so that

\[g(c) = x_0 \, a_0(c) + x_1 \, a_1(c) + ẋ_0 \, b_0(c) + ẋ_1 \, b_1(c) ,\]

with derivatives

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

The basis functions satisfy

\[\begin{aligned} g(-1) &= x_0 , & g(0) &= x_1 , & g'(-1) &= ẋ_0 , & g'(0) &= ẋ_1 . \end{aligned}\]

A ProjectionMethod is an algorithm that is applied together with a geometric integrator to enforce constraints which are not automatically satisfied by the integrator. Examples include conservation of invariants like energy or the Dirac constraint in IODEs.

Holds the solution of a geometric equation at a single time step.

It stores all the information that is passed from one time step to the next. This includes the current solution, the vectorfield of the equation, and solution data from previous time steps.

Type Parameters

• equationType: type of the geometric equation
• solutionType: type of the solution tuple
• vectorfieldType: type of the vectorfield tuple
• historyType: type of the history tuple
• internalType: type of the internal variables tuple
• paramsType: type of the parameters
• nHistory: number of previous time steps to store

Fields

• solution: a NamedTuple of OffsetVectors holding the solution of the current and previous nHistory time steps. The indices of the OffsetVector are 0...nHistory. solution[k][0] is the current solution for variable k, solution[k][1] the solution at the previous time step, and so on.
• vectorfield: a NamedTuple of OffsetVectors holding the vector field of the current and previous nHistory time steps.
• history: a NamedTuple that provides convenient access to solution and vectorfield of the current and previous time steps.
• internal: a NamedTuple for integrator-specific internal variables.
• parameters: the parameters of the equation.

Constructors

SolutionStep{equType}(ics::State, parameters::OptionalParameters; nhistory=1, internal=NamedTuple())
SolutionStep(problem::GeometricProblem; kwargs...)

The constructor SolutionStep{equType}(...) automatically constructs the appropriate solution step object from the given initial conditions ics and parameters. The internal field of the solution step is for integrator-specific internal state. The solutionstep(integrator, ...) function is a convenient wrapper to construct a SolutionStep with the correct internal variables for a given integrator.

Symplectic Euler-A Method for separable Hamiltonians.

The momentum is updated explicitly and the position is updated with the new momentum,

\[\begin{aligned} p_{n+1} &= p_{n} + h \, f (t_{n}, q_{n}) , \\ q_{n+1} &= q_{n} + h \, v (t_{n+1}, p_{n+1}) . \end{aligned}\]

This is the adjoint of SymplecticEulerB. See SymplecticEulerMethod for the separability assumption.

Symplectic Euler-B Method for separable Hamiltonians.

The position is updated explicitly and the momentum is updated with the new position,

\[\begin{aligned} q_{n+1} &= q_{n} + h \, v (t_{n}, p_{n}) , \\ p_{n+1} &= p_{n} + h \, f (t_{n+1}, q_{n+1}) . \end{aligned}\]

This is the adjoint of SymplecticEulerA. See SymplecticEulerMethod for the separability assumption.

Symplectic Euler integrator cache.

Abstract supertype of the symplectic Euler methods SymplecticEulerA and SymplecticEulerB.

Both methods are implemented for a separable Hamiltonian, that is

\[H (t, q, p) = T (t, p) + V (t, q) ,\]

so that the vector fields satisfy $v = v(t,p)$ and $f = f(t,q)$. Under this assumption the otherwise implicit partitioned scheme decouples into two explicit substeps and no nonlinear solver is required. Separability cannot be checked at runtime, so applying these methods to a non-separable Hamiltonian silently computes something that is neither symplectic nor consistent with the symplectic Euler method.

Copy the initial conditions of a EquationProblem to the current state of a solution step.

copy!(solstep::SolutionStep, sol::NamedTuple)

Copy the values from a NamedTuple sol to the current time step of the solution step.

The keys of sol must be a subset of the keys of the solution step. Only the current time step (index 0) of the solution step is modified.

Arguments

• solstep: the solution step to copy into
• sol: the named tuple containing the solution values to copy

keys(solstep::SolutionStep)

Return the keys of the state variables in the solution step.

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

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

Solve for all time steps n:

integrate!(solution, integrator)

Parts of one integration step that are common to most if not all typical integrators

parameters(solstep::SolutionStep)

Return the parameters field of the solution step, which contains the parameters of the geometric equation.

reset!(solstep::SolutionStep, Δt)

Reset the solution step for the next time step by shifting the solution history backward and advancing the time by Δt.

This function moves the current solution to the previous time step position, the previous solution to the one before that, and so on. The time variable is incremented by Δt. This prepares the solution step for computing the next time step.

Arguments

• solstep: the solution step to reset
• Δt: the time step size to advance

solution(solstep::SolutionStep, i::Int)

Return a NamedTuple with the solution at time step i, where i=0 is the current time step, i=1 is the previous time step, etc.

solution(solstep::SolutionStep)

Return the solution field of the solution step, which contains the solution vectors for all variables at the current and previous time steps.

solutionkeys(solstep::SolutionStep)

Return the keys of the solution variables in the solution step.

state(solstep::SolutionStep, i::Int)

Return a State with the state at time step i, where i=0 is the current time step, i=1 is the previous time step, etc.

state(solstep::SolutionStep)

Return the state field of the solution step, which contains the state vectors for all variables at the current and previous time steps.

update!(solstep::SolutionStep, Δ::NamedTuple)

Update the current solution in the solution step by adding increments.

This function applies increments to the current solution variables (at index 0) in the solution step. The increments are added using the appropriate add! method for each variable type, which handles different variable types correctly (e.g., compensated summation for StateWithError variables).

Arguments

• solstep::SolutionStep: The solution step to update
• Δ::NamedTuple: Named tuple containing increments for each variable to update. Keys must be a subset of the solution step's variable keys.

Returns

• solstep: The updated solution step (for method chaining)

Throws

• ArgumentError: If any key in Δ is not present in the solution step

Examples

# Create a solution step
solstep = SolutionStep{MyEquation}(initial_conditions, parameters)

# Update with increments
update!(solstep, (q = [0.1, 0.2], p = [0.05, 0.1]))

vectorfield(solstep::SolutionStep, i::Int)

Return a NamedTuple with the vectorfield at time step i, where i=0 is the current time step, i=1 is the previous time step, etc.

vectorfield(solstep::SolutionStep)

Return the vectorfield field of the solution step, which contains the vectorfield evaluations for all variables at the current and previous time steps.

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

check_solver_status(status, int)

Act on the outcome of the nonlinear solve of one time step. Every integrate_step! that solves calls this with the SimpleSolvers.NonlinearSolverStatus that solve_with_status! returned, and returns whatever it returns.

It is silent by default, and that is the point of it rather than an omission. SimpleSolvers reports a solve that did not converge itself, from the end of its own solve!, and since 0.12 the status is the programmatic counterpart of that report rather than a replacement for it — a caller that wants to act on a rejected line search reads SimpleSolvers.dominant_linesearch_outcome instead of scraping the log. Warning here as well would say the same thing twice per time step, and SimpleSolvers' own report already backs off (occurrences 1, 2, 4, 8, …) precisely because a time-stepping loop is the case that floods.

What it buys is that the status is never silently dropped: it arrives somewhere named, on every one of the call sites, and there is one place — this one — to change the policy for all of them. Overriding it is how a downstream family of methods becomes strict about non-convergence, in the same way that overriding default_options is how it changes solver options in one place:

function GeometricIntegratorsBase.check_solver_status(status, int::GeometricIntegrator{<:MyMethod})
    SimpleSolvers.isconverged(status) || throw(NonlinearSolverException("step did not converge"))
    status
end

NonlinearSolverException is the type to throw, because integrate! already catches exactly that one, warns with the time step it failed at, and returns the trajectory computed so far instead of discarding it — see src/integrate.jl. Anything else propagates and loses the solution.

Both arguments of the fallback are untyped on purpose. status is, so that this package does not have to enumerate SimpleSolvers' status types to stay reachable — every solver method produces a NonlinearSolverStatus today, but the hook has no reason to care. int is, so that every integrator reaches the default; an override narrows one or both, which is what makes it a per-method policy rather than a global one.

Compute the stages of the Crank-Nicolson method for a partitioned differential equation from the nonlinear solver solution x.

Requires v̄ and f̄ in the cache at working precision to hold $v(t_{n}, q_{n}, p_{n})$ and $f(t_{n}, q_{n}, p_{n})$, which integrate_step! computes at the beginning of every time step. Calling this function before that would silently use stale or zero values.

Compute the stages of the Crank-Nicolson method from the nonlinear solver solution x.

Requires v̄ in the cache at working precision to hold $v(t_{n}, q_{n})$, which integrate_step! computes at the beginning of every time step. Calling this function before that would silently use a stale or zero v̄.

Returns a NamedTuple with the solution of the current time step.

default_linesearch([T,] method=nothing)

The line search a method uses by default, built at the working type T where one is given.

Backtracking() is always Backtracking{Float64}, so a Float32 integration makes the Linesearch constructor do a change_precision the caller could have avoided by building the method at the working type in the first place. Both forms exist: the untyped one is the downstream hook (GeometricIntegrators calls it for DIRK) and its signature must not change.

Neither form is called inside this package. They are the hook downstream methods reach for, which is why the tests assert both signatures and the choice below.

T is constrained so that the typed method cannot quietly intercept a call meant for the untyped one: a caller passing a method type rather than an instance — default_linesearch(ImplicitEuler) — falls through to the untyped hook and gets a Backtracking{Float64}, as it did before the typed method existed, rather than an error from inside Backtracking. Real rather than AbstractFloat, because Backtracking is buildable at any real type that T(::Float64) accepts, which includes the ForwardDiff.Dual a caller differentiating through an integration would arrive with.

expand = true is deliberately not set. The expansion phase of Backtracking is opt-in upstream because it costs a merit evaluation — a full residual evaluation here — on a direction whose scale is wrong, and gains nothing on one whose scale is right: a Newton step already sits at its model minimum. Every method in this package solves with Newton().

default_options(method, problem)

The options a method is solved with unless the caller says otherwise. They reach the nonlinear solver as the keyword arguments of initsolver.

GeometricIntegrator merges them with the options it is handed — merge(default_options(method, problem), options) in src/integrator.jl — so a caller who sets one of them keeps the rest. Replacing the whole bundle instead, as this package did up to 0.4.x, meant that passing any option at all silently dropped every default not restated alongside it.

The bundle is queried on the initmethod-specialised method, so solversize below sees the concrete integrator rather than the method the caller named: LobattoIIIA(2) on a 2-dof ODE arrives here as a DIRK with solversize == 2, RadauIIA(2) as an IRK with solversize == 4.

• min_iterations = 1 — SimpleSolvers tests its stopping criteria before the first step, so without this a solve whose initial guess already satisfies them takes no step at all, and the extrapolated initial guesses used throughout this package are regularly that good.
• f_abstol scales with the size of the solver system: max(8, solversize(method, problem)) residual components at eps(datatype(problem)) apiece, because the norm of a residual sitting on its round-off floor grows with the number of components it sums. f_reltol is no substitute — it is anchored at the initial residual, so a good initial guess makes the relative gate tighter rather than looser. The floor is 8 and not 2 because the smallest solver systems are genuinely tiny: at 2eps ≈ 4.4e-16 the DIRK case above sits below its own round-off floor and stagnates, whereas 8eps ≈ 1.78e-15 clears every floor measured downstream and is what this package defaulted to before the size factor existed.
• f_stall_window — see DEFAULT_F_STALL_WINDOW.

Overriding this is how a downstream family of methods changes solver options in one place rather than at every call site. The two overrides that existed for f_abstol — one for the implicit Runge-Kutta families, one for SPARK — were both removed once the size scaling landed, since the sized value lands above the residual floor of each of them.

enforce_periodicity!(solstep::SolutionStep)

Enforce periodic boundary conditions on a state variable.

This function checks each component of the state variable for periodic boundary conditions and adjusts values that fall outside the specified range by adding or subtracting the range size. The adjustment is applied to both the current solution and all historical values to maintain consistency for initial guesses in iterative solvers.

Arguments

• solstep::SolutionStep: The solution step to enforce periodicity on

Details

This function iterates through all variables in the solution step and calls enforce_periodicity! on each one. Variables that support periodicity will have their boundary conditions enforced, while others will be unaffected.

For each component i of the state variable:

• If isperiodic(current(solstep)[s], i) is true, check if current(solstep)[s][i] is within range(current(solstep)[s], i)
• If below the range, add the range size until within bounds
• If above the range, subtract the range size until within bounds
• Apply the same adjustment to all historical values state(solstep, j)[s][i] for consistency

Performs one integration step of an integrator.

internal(solstep::SolutionStep)

Return the internal field of the solution step, which contains integrator-specific internal variables.

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.

nhistory(solstep::SolutionStep)

Return the number of previous time steps stored in the solution step.

Returns a NamedTuple with the solution of the previous time step.

AbstractVariable{T,N} is a wrapper around a AbstractArray{T,N} that provides context for the nature of the variable.

AbstractScalarVariable{T} is a wrapper around a zero-dimensional AbstractArray{T,0} that provides context for the nature of the variable, e.g., time.

AbstractStateVariable{T,N,AT} is a wrapper around a AT <: AbstractArray{T,N} that provides context for the nature of the variable, e.g., a state or a vector field.

StateVariable{T,N,AT,RT,PT} is a wrapper around a AT <: AbstractArray{T,N} that holds one of the variables consituting the state of a dynamical system.

The value field holds the actual data, the range field holds the range of valid values, and the periodic field holds a bitmask indicating which dimensions are periodic.

StateVector{DT,VT} is a vector of StateVariables, where DT is the datatype of the state and VT is the type of the vector.

Holds the solution of a geometric equation at a single time step.

It stores all the information that is required to uniquely determine the state of a systen, in particular all state variables and their corresponding vector fields.

HistoryState(st::State)

Constructs a state whose symbols are decorated by a bar to indicate a previous value of the state.