Base

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.

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

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.

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

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::NamedTuple, 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.

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

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

• AssertionError: 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]))

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

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

enforce_periodicity!(solstep::SolutionStep, ::OffsetVector{<:AbstractVariable})

No-op method for non-periodic variables.

This method is called for variable types that do not support periodicity constraints and performs no operation.

Arguments

• solstep::SolutionStep: The solution step (unused)
• ::OffsetVector{<:AbstractVariable}: The variable vector (unused)

enforce_periodicity!(solstep::SolutionStep, s::OffsetVector{<:Union{<:StateVariable,<:StateVariableWithError}})

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 containing the variable
• s::OffsetVector{<:Union{<:StateVariable,<:StateVariableWithError}}: The state variable to enforce periodicity on

Details

For each component i of the state variable:

• If isperiodic(s[0], i) is true, check if s[0][i] is within range(s[0], 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 s[j][i] for consistency

enforce_periodicity!(solstep::SolutionStep)

Enforce periodic boundary conditions on all solution variables.

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.

Arguments

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

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

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

history(solstep::SolutionStep)

Return the history field of the solution step, which provides convenient access to both solution and vectorfield data at the current and previous time steps.

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.

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.

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.