Equations

ODE: Ordinary Differential Equation

Ordinary differential equations define an initial value problem of the form

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

with vector field $v$.

Parameters

• vType <: Callable: type of v
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the vector field
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

ODE(v, invariants, parameters, periodicity)
ODE(v; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

Function Definitions

The function v providing the vector field must have the interface

function v(v, t, q, params)
    v[1] = ...
    v[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the vector which holds the result of evaluating the vector field $v$ on t and q, and params is a NamedTuple of additional parameters on which the vector field may depend.

PODE: Partitioned Ordinary Differential Equation

A partitioned ordinary differential equation is an initial value problem of the form

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

with vector fields $v$ and $f$.

Parameters

• vType <: Callable: type of v
• fType <: Callable: type of f
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the vector field $v$
• f: function computing the vector field $f$
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

PODE(v, f, invariants, parameters, periodicity)
PODE(v, f; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

Function Definitions

The functions v and f must have the interface

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f(f, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

where t is the current time, q and p are the current solution vectors, v and f are the vectors which hold the result of evaluating the vector fields $v$ and $f$ on t, q and p, and params is a NamedTuple of additional parameters.

HODE: Hamiltonian Ordinary Differential Equation

A canonical Hamiltonian system of equations is special case of a partitioned ordinary differential equation,

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

with vector fields $v$ and $f$, given by

\[\begin{aligned} v &= \frac{\partial H}{\partial p} , & f &= - \frac{\partial H}{\partial q} . \end{aligned}\]

Parameters

• vType <: Callable: type of v
• fType <: Callable: type of f
• hamType <: Callable: Hamiltonian type
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the vector field $v$
• f: function computing the vector field $f$
• hamiltonian: function computing the Hamiltonian $H$ (usually the total energy of the system)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

HODE(v, f, hamiltonian, invariants, parameters, periodicity)
HODE(v, f, hamiltonian; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

Function Definitions

The functions v, f and hamiltonian must have the interface

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f(f, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function hamiltonian(t, q, p, params)
    return ...
end

where t is the current time, q and p are the current solution vectors, v and f are the vectors which hold the result of evaluating the vector fields on t, q and p, and params is a NamedTuple of additional parameters.

IODE: Implicit Ordinary Differential Equation

An implicit ordinary differential equations is an initial value problem of the form

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

with force field $f$, the momentum defined by $p$. This is a special case of a differential algebraic equation with dynamical variables $(q,p)$ and algebraic variable $v$, that is determined such that the constraint $p(t) = ϑ(t, q(t), v(t))$ is satisfied.

Many integrators perform a projection step in order to enforce this constraint. To this end, the system is extended to

\[\begin{aligned} \dot{q} (t) &= v(t) + λ(t) , \\ \dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), v(t), λ(t)) , \\ p(t) &= ϑ(t, q(t), v(t)) , \end{aligned}\]

where the vector field defining the projection step is usually given as

\[\begin{aligned} g(t, q(t), v(t), λ(t)) &= λ(t) \cdot \nabla ϑ(t, q(t), v(t)) . \end{aligned}\]

Parameters

• ϑType <: Callable: type of ϑ
• fType <: Callable: type of f
• gType <: Callable: type of g
• v̄Type <: Callable: type of v̄
• f̄Type <: Callable: type of f̄
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• ϑ: function determining the momentum
• f: function computing the vector field
• g: function determining the projection, given by $\nabla \vartheta (t,q,v) \cdot \lambda$
• v̄: function computing an initial guess for the velocity field $v$ (optional)
• f̄: function computing an initial guess for the force field $f$ (optional)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

IODE(ϑ, f, g, v̄, f̄, invariants, parameters, periodicity)
IODE(ϑ, f, g; v̄ = _iode_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

where

_iode_default_v̄(v, t, q, params) = nothing

The functions ϑ, f and g compute the momentum and the vector fields, respectively.

Function Definitions

The functions ϑ and f must have the interface

function ϑ(p, t, q, v)
    p[1] = ...
    p[2] = ...
    ...
end

and

function f(f, t, q, v)
    f[1] = ...
    f[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the current velocity and f and p are the vectors which hold the result of evaluating the functions $f$ and $ϑ$ on t, q and v. In addition, the functions g, v̄ and f̄ are specified by

function g(g, t, q, v, λ)
    g[1] = ...
    g[2] = ...
    ...
end

function v̄(v, t, q, p)
    v[1] = ...
    v[2] = ...
    ...
end

function f̄(f, t, q, v)
    f[1] = ...
    f[2] = ...
    ...
end

The function g is used in projection methods that enforce $p = ϑ(q)$. The functions v̄ and f̄ are used for initial guesses in nonlinear implicit solvers.

LODE: Lagrangian Ordinary Differential Equation

A Lagrangian system of equations is a special case of an implicit ordinary differential equations, that is an implicit initial value problem of the form

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

with momentum $p$ and force field $f$, given by

\[\begin{aligned} p &= \frac{\partial L}{\partial v} , & f &= \frac{\partial L}{\partial q} . \end{aligned}\]

This is a special case of an implicit ordinary differential equation, that is defined by a Lagrangian, as well as a special case of a differential algebraic equation with dynamical variables $(q,p)$ and algebraic variable $v$, that is determined such that the constraint $p(t) = ϑ(t, q(t), v(t))$ is satisfied.

Many integrators perform a projection step in order to enforce this constraint. To this end, the system is extended to

\[\begin{aligned} \dot{q} (t) &= v(t) + \lambda(t) , \\ \dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), v(t), \lambda(t)) , \\ p(t) &= ϑ(t, q(t), v(t)) , \end{aligned}\]

where the vector field defining the projection step is usually given as

\[\begin{aligned} g(t, q(t), v(t), λ(t)) &= λ(t) \cdot \nabla ϑ(t, q(t), v(t)) . \end{aligned}\]

Parameters

• ϑType <: Function: type of ϑ
• fType <: Function: type of f
• gType <: Function: type of g
• ωType <: Function: type of ω
• v̄Type <: Function: type of v̄
• f̄Type <: Function: type of f̄
• lagType <: Function: Lagrangian type
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• ϑ: function determining the momentum
• f: function computing the vector field
• g: function determining the projection, given by $\nabla \vartheta (q) \cdot \lambda$
• ω: function computing the symplectic matrix
• v̄: function computing an initial guess for the velocity field $v$ (optional)
• f̄: function computing an initial guess for the force field $f$ (optional)
• lagrangian: function computing the Lagrangian $L$
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

LODE(ϑ, f, g, ω, l, v̄, f̄, invariants, parameters, periodicity)
LODE(ϑ, f, g, ω, l; v̄ = _lode_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

where

_lode_default_v̄(v, t, q, params) = nothing

Function Definitions

The functions ϑ and f must have the interface

function ϑ(p, t, q, v)
    p[1] = ...
    p[2] = ...
    ...
end

and

function f(f, t, q, v)
    f[1] = ...
    f[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the current velocity and f and p are the vectors which hold the result of evaluating the functions $f$ and $ϑ$ on t, q and v. The funtions g, v̄ and f̄ are specified by

function g(g, t, q, λ)
    g[1] = ...
    g[2] = ...
    ...
end

and

function v̄(v, t, q, p)
    v[1] = ...
    v[2] = ...
    ...
end

function f̄(f, t, q, v)
    f[1] = ...
    f[2] = ...
    ...
end

The function g is used in projection methods that enforce $p = ϑ(q)$. The functions v̄ and f̄ are used for initial guesses in nonlinear implicit solvers. Finally, the functions ω and l, computing the symplectic matrix and the Lagrangian, have the following signature

function ω(f, t, q, v, params)
    ω[1,1] = ...
    ω[1,2] = ...
    ...
end

function l(t, q, v, params)
    return ...
end

SODE: Split Ordinary Differential Equation

Defines an initial value problem

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

with vector field $v$, initial condition $q_{0}$ and the solution $q$ taking values in $\mathbb{R}^{d}$. Here, the vector field $v$ is given as a sum of vector fields

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

The dynamical variables $q$ with initial condition $q_{0}$ take values in $\mathbb{R}^{d}$.

Parameters

• vType <: Union{Tuple,Nothing}: type of v
• qType <: Union{Tuple,Nothing}: type of q
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: tuple of functions computing the vector fields for each substep
• q: tuple of functions computing the solutions for each substep
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

SODE(v, invariants, parameters, periodicity)
SODE(v; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity())
SODE(v, q, invariants, parameters, periodicity)
SODE(v, q; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity())

Function Definitions

The functions v_i providing the vector field must have the interface

function v_i(v, t, q, params)
    v[1] = ...
    v[2] = ...
    ...
end

and the functions q_i providing the solutions must have the interface

function q_i(q₁, t₁, q₀, t₀, params)
    q₁[1] = q₀[1] + ...
    q₁[2] = q₀[2] + ...
    ...
end

where t₀ is the current time, q₀ is the current solution vector, q₁ is the new solution vector at time t₁, holding the result of computing one substep

The fact that the function v returns the solution and not just the vector field for each substep increases the flexibility for the use of splitting methods, e.g., it allows to use another integrator for solving substeps. with the vector field $v_i$.

DAE: Differential Algebraic Equation

Defines a differential algebraic initial value problem

\[\begin{aligned} \dot{q} (t) &= v(t, q(t)) + u(t, q(t), \lambda(t)) , \\ 0 &= \phi (t, q(t)) , \end{aligned}\]

with vector field $v$, projection $u$, algebraic constraint $\phi=0$.

Some integrators also enforce the secondary constraint $\psi$, that is the time derivative of the algebraic constraint $\phi$. In this case, the system of equations is modified as follows

\[\begin{aligned} \dot{q} (t) &= v(t, q(t)) + u(t, q(t), \lambda(t)) + \bar{u} (t, q(t), \dot{q} (t), \dot{p} (t), \mu(t)) , \\ 0 &= \phi (t, q(t)) , \\ 0 &= \psi (t, q(t), \dot{q} (t)) . \end{aligned}\]

Parameters

• vType <: Callable: type of v
• uType <: Callable: type of u
• ϕType <: Callable: type of ϕ
• ūType <: OptionalCallable: type of ū
• ψType <: OptionalCallable: type of ψ
• v̄Type <: Callable: type of v̄
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the vector field v(v, t, q, params)
• u: function computing the projection u(u, t, q, λ, params)
• ϕ: algebraic constraint ϕ(ϕ, t, q, params)
• ū: function computing the secondary projection field ū(ū, t, q, λ, params) (optional)
• ψ: secondary constraint ψ(ψ, t, q, v, params) (optional)
• v̄: function computing an initial guess for the velocity field $v$ (defaults to v)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

DAE(v, u, ϕ, ū, ψ, v̄, invariants, parameters, periodicity)
DAE(v, u, ϕ, ū, ψ; kwargs...)
DAE(v, u, ϕ; kwargs...)

The functions v and u compute the vector field and the projection, respectively, ϕ provides the algebraic constraint. The functions ψ and ū are optional and provide the secondary constraint, that is the time derivative of the algebraic constraint, and the corresponding projection.

Function Definitions

The functions are defined by

The functions v, u and ϕ must have the interface

function v(v, t, q, params)
    v[1] = ...
    v[2] = ...
    ...
end

function u(u, t, q, λ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ϕ(ϕ, t, q, params)
    ϕ[1] = ...
end

where t is the current time, q and λ are the current solution vectors, and v, u and ϕ are the vectors which hold the result of evaluating the vector field $v$, the projection $u$ and the algebraic constraint $\phi$ on t, q and λ.

Some integrators also enforce the secondary constraint $\psi$ and require the following additional functions

function ū(u, t, q, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ψ(ψ, t, q, v, params)
    ψ[1] = ...
end

The DAE is created by

equ = DAE(v, u, ϕ)

or

equ = DAE(v, u, ϕ, ū, ψ)

PDAE: Partitioned Differential Algebraic Equation

A partitioned differential algebraic equation has the form

\[\begin{aligned} \dot{q} (t) &= v(t, q(t), p(t)) + u(t, q(t), p(t), \lambda(t)) , \\ \dot{p} (t) &= f(t, q(t), p(t)) + g(t, q(t), p(t), \lambda(t)) , \\ 0 &= \phi (t, q(t), p(t)) , \end{aligned}\]

with vector fields $v$ and $f$, projection $u$ and $g$, algebraic constraint $\phi=0$.

Some integrators also enforce the secondary constraint $\psi$, that is the time derivative of the algebraic constraint $\phi$. In this case, the system of equations is modified as follows

\[\begin{aligned} \dot{q} (t) &= v(t, q(t), p(t)) + u(t, q(t), p(t), \lambda(t)) + \bar{u} (t, q(t), p(t), \mu(t)) , \\ \dot{p} (t) &= f(t, q(t), p(t)) + g(t, q(t), p(t), \lambda(t)) + \bar{g} (t, q(t), p(t), \mu(t)) , \\ 0 &= \phi (t, q(t), p(t)) , \\ 0 &= \psi (t, q(t), p(t), \dot{q} (t), \dot{p} (t)) . \end{aligned}\]

Parameters

• vType <: Callable: type of v
• fType <: Callable: type of f
• uType <: Callable: type of u
• gType <: Callable: type of g
• ϕType <: Callable: type of ϕ
• ūType <: Callable: type of ū
• ḡType <: Callable: type of ḡ
• ψType <: Callable: type of ψ
• v̄Type <: Callable: type of v̄
• f̄Type <: Callable: type of f̄
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the vector field $v$
• f: function computing the vector field $f$
• u: function computing the projection for $q$
• g: function computing the projection for $p$
• ϕ: algebraic constraints
• ū: function computing the secondary projection field $\bar{u}$ (optional)
• ḡ: function computing the secondary projection field $\bar{g}$ (optional)
• ψ: secondary constraints (optional)
• v̄: function computing an initial guess for the velocity field $v$ (optional, defaults to v)
• f̄: function computing an initial guess for the force field $f$ (optional, defaults to f)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

PDAE(v, f, u, g, ϕ, ū, ḡ, ψ, v̄, f̄, invariants, parameters, periodicity)
PDAE(v, f, u, g, ϕ, ū, ḡ, ψ; kwargs...)
PDAE(v, f, u, g, ϕ; kwargs...)

The functions v and f compute the vector field, u and g compute the projections, and ϕ provides the algebraic constraint. The functions ψ, ū and ḡ are optional and provide the secondary constraint and the corresponding projection.

Function Definitions

The functions are defined by

The functions v, f, u, g and ϕ must have the interface

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f(g, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function u(u, t, q, p, λ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function g(g, t, q, p, λ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ϕ(ϕ, t, q, p, params)
    ϕ[1] = ...
end

where t is the current time, q, p and λ are the current solution vectors, v, f, u and g are the vectors which hold the result of evaluating the vector fields $v$ and $f$, the projections $u$ and $g$, and ϕ holds the algebraic constraint $\phi$, evaluated on t, q, p and λ.

Some integrators also enforce the secondary constraint $\psi$ and require the following additional functions

function ū(u, t, q, p, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ḡ(g, t, q, p, μ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ψ(ψ, t, q, p, v, f, params)
    ψ[1] = ...
end

The PDAE is created by

equ = PDAE(v, f, u, g, ϕ)

or

equ = PDAE(v, f, u, g, ϕ, ū, ḡ, ψ)

HDAE: Hamiltonian Differential Algebraic Equation

A Hamiltonian differential algebraic is an initial value problem, that is a canonical Hamiltonian system of equations subject to Dirac constraints,

\[\begin{aligned} \dot{q} (t) &= v(t, q(t), p(t)) + u(t, q(t), p(t), \lambda(t)) + \bar{u} (t, q(t), p(t), \mu(t)) , \\ \dot{p} (t) &= f(t, q(t), p(t)) + g(t, q(t), p(t), \lambda(t)) + \bar{g} (t, q(t), p(t), \mu(t)) , \\ 0 &= \phi (t, q(t), p(t)) , \\ 0 &= \psi (t, q(t), p(t), \dot{q}(t), \dot{p}(t)) , \end{aligned}\]

with vector fields $v$, $u$, $\bar{u}$ and $f$, $g$, $\bar{g}$, primary constraint $\phi(q,p)=0$ and secondary constraint $\psi(q,p,\dot{q},\dot{p})=0$.

Parameters

• vType <: Callable: type of v
• fType <: Callable: type of f
• uType <: Callable: type of u
• gType <: Callable: type of g
• ϕType <: Callable: type of ϕ
• ūType <: Callable: type of ū
• ḡType <: Callable: type of ḡ
• ψType <: Callable: type of ψ
• v̄Type <: Callable: type of v̄
• f̄Type <: Callable: type of f̄
• hamType <: Callable: Hamiltonian type
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the Hamiltonian vector field $v$
• f: function computing the Hamiltonian vector field $f$
• u: function computing the projection for $q$
• g: function computing the primary projection field $g$
• ϕ: primary constraints
• ū: function computing the secondary projection field $\bar{u}$ (optional)
• ḡ: function computing the secondary projection field $\bar{g}$ (optional)
• ψ: secondary constraints (optional)
• v̄: function computing an initial guess for the velocity field $v$ (optional, defaults to v)
• f̄: function computing an initial guess for the force field $f$ (optional, defaults to f)
• hamiltonian: function computing the Hamiltonian $H$ (usually the total energy of the system)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

HDAE(v, f, u, g, ϕ, ū, ḡ, ψ, h, v̄, f̄, invariants, parameters, periodicity)
HDAE(v, f, u, g, ϕ, ū, ḡ, ψ, h; kwargs...)
HDAE(v, f, u, g, ϕ, h; kwargs...)

The functions v and f compute the vector field, u and g compute the projections, ϕ provides the algebraic constraint and h the Hamiltonian. The functions ψ, ū and ḡ are optional and provide the secondary constraint, that is the time derivative of the algebraic constraint, and the corresponding projection.

Function Definitions

The functions are defined by

The functions v, f, u, g, ϕ and h must have the interface

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f(g, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function u(u, t, q, p, λ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function g(g, t, q, p, λ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ϕ(ϕ, t, q, p, params)
    ϕ[1] = ...
end

function h(t, q, p, params)
    ...
end

where t is the current time, q, p, λ and μ are the current solution vectors, v, f, u and g are the vectors which hold the result of evaluating the vector fields $v$ and $f$, the projections on the primary constraint $u$ and $g$, ϕ holds the algebraic constraint $\phi$, and h returns the Hamiltonian of the system, all evaluated on t, q, p and λ.

Some integrators also enforce the secondary constraint $\psi$ and require the following additional functions

function ū(u, t, q, p, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ḡ(g, t, q, p, μ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ψ(ψ, t, q, p, v, f, params)
    ψ[1] = ...
end

The HDAE is created by

equ = HDAE(v, f, u, g, ϕ, h)

or

equ = HDAE(v, f, u, g, ϕ, ū, ḡ, ψ, h)

IDAE: Implicit Differential Algebraic Equation

An implicit differential algebraic initial value problem takes the form

\[\begin{aligned} \dot{q} (t) &= v(t) + u(t, q(t), v(t), p(t), \lambda(t)) , \\ \dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), v(t), p(t), \lambda(t)) , \\ p(t) &= ϑ(t, q(t), v(t)) , \\ 0 &= \phi (t, q(t), v(t), p(t)) , \end{aligned}\]

with force field $f$, the momentum defined by $ϑ$, projections $u$ and $g$, algebraic constraint $\phi(t,q,v,p)=0$.

Some integrators also enforce the secondary constraint $\psi$, that is the time derivative of the algebraic constraint $\phi$. In this case, the system of equations is modified as follows

\[\begin{aligned} \dot{q} (t) &= v(t) + u(t, q(t), v(t), p(t), \lambda(t)) + \bar{u} (t, q(t), v(t), p(t), \mu(t)) , \\ \dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), v(t), p(t), \lambda(t)) + \bar{g} (t, q(t), v(t), p(t), \mu(t)) , \\ p(t) &= ϑ(t, q(t), v(t)) , && \\ 0 &= \phi (t, q(t), v(t), p(t)) , \\ 0 &= \psi (t, q(t), v(t), p(t), \dot{q} (t), \dot{p} (t)) . \end{aligned}\]

Parameters

• ϑType <: Callable: type of ϑ
• fType <: Callable: type of f
• uType <: Callable: type of u
• gType <: Callable: type of g
• ϕType <: Callable: type of ϕ
• ūType <: Callable: type of ū
• ḡType <: Callable: type of ḡ
• ψType <: Callable: type of ψ
• v̄Type <: Callable: type of v̄
• f̄Type <: Callable: type of f̄
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• ϑ: function determining the momentum
• f: function computing the vector field $f$
• u: function computing the projection for $q$
• g: function computing the projection for $p$
• ϕ: algebraic constraints
• ū: function computing the secondary projection field $\bar{u}$ (optional)
• ḡ: function computing the secondary projection field $\bar{g}$ (optional)
• ψ: secondary constraints (optional)
• v̄: function computing an initial guess for the velocity field $v$ (optional)
• f̄: function computing an initial guess for the force field $f$ (optional)
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

IDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, v̄, f̄, invariants, parameters, periodicity)
IDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ; v̄ = _idae_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())
IDAE(ϑ, f, u, g, ϕ; v̄ = _idae_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

where

_idae_default_v̄(v, t, q, params) = nothing

The function ϑ computes the momentum, f computes the force field, u and g compute the projections, and ϕ provides the algebraic constraint. The functions ψ, ū and ḡ are optional and provide the secondary constraint, that is the time derivative of the algebraic constraint, and the corresponding projection.

Function Definitions

The functions ϑ and f must have the interface

function ϑ(p, t, q, v, params)
    p[1] = ...
    p[2] = ...
    ...
end

function f(f, t, q, v, params)
    f[1] = ...
    f[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the current velocity and f and p are the vectors which hold the result of evaluating the functions $f$ and $ϑ$ on t, q and v. The funtions g, v̄ and f̄ are specified by

function u(u, t, q, v, p, λ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function g(g, t, q, v, p, λ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function v̄(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f̄(f, t, q, v, params)
    f[1] = ...
    f[2] = ...
    ...
end

Some integrators also enforce the secondary constraint $\psi$ and require the following additional functions

function ū(u, t, q, v, p, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ḡ(g, t, q, v, p, μ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ψ(ψ, t, q, v, p, q̇, ṗ, params)
    ψ[1] = ...
end

LDAE: Lagrangian Differential Algebraic Equation

A special case of an implicit initial value problem is a Lagrangian differential algebraic equation of the form

\[\begin{aligned} \dot{q} (t) &= v(t) + u(t, q(t), v(t), p(t), \lambda(t)) + \bar{u} (t, q(t), v(t), p(t), \mu(t)) , \\ \dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), v(t), p(t), \lambda(t)) + \bar{g} (t, q(t), v(t), p(t), \mu(t)) , \\ p(t) &= ϑ(t, q(t), v(t)) , \\ 0 &= \phi (t, q(t), v(t), p(t)) , \\ 0 &= \psi (t, q(t), v(t), p(t), \dot{q}(t), \dot{p}(t)) , \end{aligned}\]

with momentum $p$ and force field $f$, given by

\[\begin{aligned} p &= \frac{\partial L}{\partial v} (q,v) , & f &= \frac{\partial L}{\partial q} (q,v) , \end{aligned}\]

projection fields $u$, $\bar{u}$ and $g$, $\bar{g}$. This is a special case of a differential algebraic equation with dynamical variables $(q,p)$ and algebraic variables $v$, $\lambda$ and $\mu$.

Parameters

• ϑType <: Callable: type of ϑ
• fType <: Callable: type of f
• uType <: Callable: type of u
• gType <: Callable: type of g
• ϕType <: Callable: type of ϕ
• ūType <: Callable: type of ū
• ḡType <: Callable: type of ḡ
• ψType <: Callable: type of ψ
• ωType <: Callable: type of ω
• v̄Type <: Callable: type of v̄
• f̄Type <: Callable: type of f̄
• lagType <: Callable: Lagrangian type
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• f: function computing the vector field
• u: function computing the projection for $q$, for a degenerate system given by $\lambda$
• g: function computing the projection for $p$, for a degenerate system given by $\nabla \vartheta (q) \cdot \lambda$
• ϕ: primary constraints, for a degenerate system given by $p - \vartheta (t,q)$
• ū: function computing the secondary projection field $\bar{u}$, for a degenerate system given by $\lambda$ (optional)
• ḡ: function computing the secondary projection field $\bar{g}$, for a degenerate system given by $\lambda \cdot \nabla \vartheta (t,q)$ (optional)
• ψ: secondary constraints, for a degenerate system given by $\dot{p} - \dot{q} \cdot \nabla \vartheta (t,q)$ (optional)
• ω: function computing the symplectic matrix
• v̄: function computing an initial guess for the velocity field $v$ (optional)
• f̄: function computing an initial guess for the force field $f$ (optional)
• lagrangian: function computing the Lagrangian $L$
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, lagrangian, invariants, parameters, periodicity)
LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, lagrangian; v̄ = _ldae_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())
LDAE(ϑ, f, u, g, ϕ, ω, lagrangian; v̄ = _ldae_default_v̄, f̄ = f, invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

where

_ldae_default_v̄(v, t, q, params) = nothing

The function ϑ computes the momentum, f computes the force field, u and g compute the projections, and ϕ provides the algebraic constraint. The functions ψ, ū and ḡ are optional and provide the secondary constraint, that is the time derivative of the algebraic constraint, and the corresponding projection.

Function Definitions

The functions ϑ and f must have the interface

function ϑ(p, t, q, v, params)
    p[1] = ...
    p[2] = ...
    ...
end

function f(f, t, q, v, params)
    f[1] = ...
    f[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the current velocity and f and p are the vectors which hold the result of evaluating the functions $f$ and $ϑ$ on t, q and v. The funtions g, v̄ and f̄ are specified by

function u(u, t, q, v, p, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function g(g, t, q, v, p, μ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function v̄(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f̄(f, t, q, v, params)
    f[1] = ...
    f[2] = ...
    ...
end

and the functions ω and l, computing the symplectic matrix and the Lagrangian, have the following signature

function ω(f, t, q, v, params)
    ω[1,1] = ...
    ω[1,2] = ...
    ...
end

function l(t, q, v, params)
    return ...
end

Some integrators also enforce the secondary constraint $\psi$ and require the following additional functions

function ū(u, t, q, v, p, μ, params)
    u[1] = ...
    u[2] = ...
    ...
end

function ḡ(g, t, q, v, p, μ, params)
    g[1] = ...
    g[2] = ...
    ...
end

function ψ(ψ, t, q, v, p, q̇, ṗ, params)
    ψ[1] = ...
end

SDE: Stratonovich Stochastic Differential Equation

Defines a stochastic differential initial value problem

\[\begin{aligned} dq (t) &= v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} , \end{aligned}\]

with drift vector field $v$, diffusion matrix $B$, initial conditions $q_{0}$, the dynamical variable $q$ taking values in $\mathbb{R}^{d}$, and the m-dimensional Wiener process W

Parameters

• vType <: Callable: type of v
• BType <: Callable: type of B
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the deterministic vector field
• B: function computing the d x m diffusion matrix
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

SDE(v, B, invariants, parameters, periodicity)
SDE(v, B; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

The functions v and B, providing the drift vector field and diffusion matrix. The function v must have the interface

function v(v, t, q, params)
    v[1] = ...
    v[2] = ...
    ...
end

where t is the current time, q is the current solution vector, v is the vector which holds the result of evaluating the vector field $v$ on t and q, and params are additional parameters. The function B should have a method with interface

function B(B, t, q, params)
    B[1,1] = ...
    ...
end

PSDE: Stratonovich Partitioned Stochastic Differential Equation

A partitioned stochastic differential equations is an initial value problem of the form

\[\begin{aligned} dq (t) &= v(t, q(t), p(t)) \, dt + B(t, q(t), p(t)) \circ dW , & q(t_{0}) &= q_{0} , \\ dp (t) &= f(t, q(t), p(t)) \, dt + G(t, q(t), p(t)) \circ dW , & p(t_{0}) &= p_{0} \end{aligned}\]

with the drift vector fields $v$ and $f$, diffusion matrices $B$ and $G$, initial conditions $q_{0}$ and $p_{0}$, the dynamical variables $(q,p)$ taking values in $\mathbb{R}^{d} \times \mathbb{R}^{d}$, and the m-dimensional Wiener process W

Parameters

• vType <: Callable: type of v
• fType <: Callable: type of f
• BType <: Callable: type of B
• GType <: Callable: type of G
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v: function computing the drift vector field for the position variable $q$
• f: function computing the drift vector field for the momentum variable $p$
• B: function computing the d x m diffusion matrix for the position variable $q$
• G: function computing the d x m diffusion matrix for the momentum variable $p$
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

PSDE(v, f, B, G, invariants, parameters, periodicity)
PSDE(v, f, B, G; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

The functions v, f, B and G, providing the drift vector fields and diffusion matrices, each take five arguments, v(v, t, q, p, params), f(f, t, q, p, params), B(B, t, q, p, params) and G(G, t, q, p, params), where t is the current time, (q, p) is the current solution, and v, f, B and G are the variables which hold the result of evaluating the vector fields $v$, $f$ and the matrices $B$, $G$ on t and (q,p), and params are optional parameters.

The corresponding methods should have the following signatures:

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f(f, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function B(B, t, q, p, params)
    B[1,1] = ...
    ...
end

function G(G, t, q, p, params)
    G[1,1] = ...
    ...
end

SPSDE: Stratonovich Split Partitioned Stochastic Differential Equation

Defines a partitioned stochastic differential initial value problem

\[\begin{aligned} dq (t) &= v(t, q(t), p(t)) \, dt + B(t, q(t), p(t)) \circ dW , & q(t_{0}) &= q_{0} , \\ dp (t) &= [ f_1(t, q(t), p(t)) + f_2(t, q(t), p(t)) ] \, dt + [ G_1(t, q(t), p(t)) + G_2(t, q(t), p(t)) ] \circ dW , & p(t_{0}) &= p_{0} , \end{aligned}\]

with the drift vector fields $v$ and $f_i$, diffusion matrices $B$ and $G_i$, initial conditions $q_{0}$ and $p_{0}$, the dynamical variables $(q,p)$ taking values in $\mathbb{R}^{d} \times \mathbb{R}^{d}$, and the m-dimensional Wiener process W

Parameters

• vType <: Function: type of v
• f1Type <: Function: type of f1
• f2Type <: Function: type of f2
• BType <: Function: type of B
• G1Type <: Function: type of G1
• G2Type <: Function: type of G2
• invType <: OptionalInvariants: invariants type
• parType <: OptionalParameters: parameters type
• perType <: OptionalPeriodicity: periodicity type

Fields

• v : function computing the drift vector field for the position variable $q$
• f1: function computing the drift vector field for the momentum variable $p$
• f2: function computing the drift vector field for the momentum variable $p$
• B : function computing the d x m diffusion matrix for the position variable $q$
• G1: function computing the d x m diffusion matrix for the momentum variable $p$
• G2: function computing the d x m diffusion matrix for the momentum variable $p$
• invariants: functions for the computation of invariants, either a NamedTuple containing the equation's invariants or NullInvariants
• parameters: type constraints for parameters, either a NamedTuple containing the equation's parameters or NullParameters
• periodicity: determines the periodicity of the state vector q for cutting periodic solutions, either a AbstractArray or NullPeriodicity

Constructors

SPSDE(v, f1, f2, B, G1, G2, invariants, parameters, periodicity)
SPSDE(v, f1, f2, B, G1, G2; invariants = NullInvariants(), parameters = NullParameters(), periodicity = NullPeriodicity())

The functions v, f1, f2, B, G1 and G2, providing the drift vector fields and diffusion matrices, all take five arguments, (out, t, q, p, params).

function v(v, t, q, p, params)
    v[1] = ...
    v[2] = ...
    ...
end

function f1(f, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function f2(f, t, q, p, params)
    f[1] = ...
    f[2] = ...
    ...
end

function B(B, t, q, p, params)
    B[1,1] = ...
    ...
end

function G1(G, t, q, p, params)
    G[1,1] = ...
    ...
end

function G2(G, t, q, p, params)
    G[1,1] = ...
    ...
end

where t is the current time, (q,p) is the current solution vector, and v, f, B and G are the variables which hold the result of evaluating the vector fields $v$, $f$ and the matrices $B_i$, $G_i$ on (t,q,p).