Equation and Problem Types

In GeometricIntegrators.jl we support three basic types of equations:

• ordinary differential equations (ODEs),
• differential algebraic equations (DAEs),
• stochastic differential equations (SDEs).

For each type, there are several subtypes:

• standard equations (ODEProblem, DAEProblem, SDEProblem),
• implicit equations (IODEProblem, IDAEProblem),
• partitioned equations (PODEProblem, PDAEProblem, PSDEProblem),
• Hamiltonian equations (HODEProblem, HDAEProblem),
• Lagrangian equations (LODEProblem, LDAEProblem),
• split equations (SODEProblem, SPSDEProblem).

Each equation holds a number of functions determining the vector field, constraints, and possibly additional information like parameters, periodicity, invariants and the Hamiltonian or Lagrangian.

To each equation type there exists a corresponding problem type, which holds the equation, initial conditions, parameters, a time span to integrate over, as well as a time step (which is typically fixed in GeometricIntegrators). In addition, these problem types provide some convenience constructors to generate the equation and the problem at once.

All of the equation and problem types are defined in the GeometricEquations.jl package. The GeometricProblems.jl package implements a number of example problems for testing and benchmarking.

Keyword Arguments

All equation and problem types take the following keyword arguments:

• invariants = NullInvariants()
• parameters = NullParameters()
• periodicity = NullPeriodicity()

If not set to their corresponding Null types, the user needs to pass a NamedTuple whose values are

• functions for invariants,
• arbitrary data structures for parameters,
• the same data structure as the solution for periodicity.

The latter should be zero everywhere, except for those components, that are periodic, i.e., whose value are supposed to stay within a range (0, max). Support for ranges starting with other values than zero is currently missing but can be added if demand arises.

Ordinary Differential Equations (ODEs)

Ordinary differential equations define an initial value problem of the form

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

with vector field $v$.

Example: Harmonic Oscillator

As an example, let us consider the harmonic oscillator. The dynamical equations are given by

\[\dot{q} (t) = \begin{pmatrix} 0 & 1 \\ -k & 0 \\ \end{pmatrix} q(t) , \qquad q \in \mathbb{R}^{2} .\]

In order to create an ODEProblem for the harmonic oscillator, we need to write the following code:

function v(v, t, x, params)
    v[1] = x[2]
    v[2] = - params.k * x[1]
end

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [0.5, 0.0]

prob = ODEProblem(v, tspan, tstep, q₀; parameters = (k = 0.5,))

Geometric Equation Problem for Ordinary Differential Equation (ODE)

 with vector field
   v = v

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0])

 Parameters: 
   (k = 0.5,)

The energy of the harmonic oscillator is preserved, so we can add it as an invariant,

energy(t, q, params) = q[2]^2 / 2 + params.k * q[1]^2 / 2

prob = ODEProblem(v, tspan, tstep, q₀; parameters = (k = 0.5,), invariants = (h=energy,))

Geometric Equation Problem for Ordinary Differential Equation (ODE)

 with vector field
   v = v

 Invariants: 
   (h = Main.energy,)

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0])

 Parameters: 
   (k = 0.5,)

Partitioned Ordinary Differential Equations (PODEs)

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

Example: Harmonic Oscillator

As an example, let us consider the harmonic oscillator. The dynamical equations are given by

\[\begin{aligned} \dot{q} (t) &= p (t) \\ \dot{p} (t) &= - k \, q(t) . \end{aligned}\]

In order to create a PODEProblem for the harmonic oscillator, we need to write the following code:

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

function f(f, t, q, p, params)
    f[1] = - params.k * q[1]
end

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [0.5]
p₀ = [0.0]

prob = PODEProblem(v, f, tspan, tstep, q₀, p₀; parameters = (k = 0.5,))

Geometric Equation Problem for Partitioned Ordinary Differential Equation (PODE)

 with vector fields
   v = v
   f = f

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0])

 Parameters: 
   (k = 0.5,)

The energy of the harmonic oscillator is preserved, so we can add it as an invariant,

energy(t, q, p, params) = p[1]^2 / 2 + params.k * q[1]^2 / 2

prob = PODEProblem(v, f, tspan, tstep, q₀, p₀; parameters = (k = 0.5,), invariants = (h=energy,))

Geometric Equation Problem for Partitioned Ordinary Differential Equation (PODE)

 with vector fields
   v = v
   f = f

 Invariants: 
   (h = Main.energy,)

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0])

 Parameters: 
   (k = 0.5,)

Hamiltonian Ordinary Differential Equations (HODEs)

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

Example: Harmonic Oscillator

As an example, let us consider the harmonic oscillator. The dynamical equations are given by

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

which can also be written as

\[\begin{pmatrix} \dot{q} (t) \\ \dot{p} (t) \\ \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \\ \end{pmatrix} \nabla H( q(t) , p(t) ) , \qquad H(q,p) = \frac{p^2}{2} + k \, \frac{q^2}{2} ,\]

where $H$ is the Hamiltonian, i.e., the total energy of the system.

In order to create a HODEProblem for the harmonic oscillator, we need to write the following code:

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

function f(f, t, q, p, params)
    f[1] = - params.k * q[1]
end

h(t, q, p, params) = p[1]^2 / 2 + params.k * q[1]^2 / 2

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [0.5]
p₀ = [0.0]

prob = HODEProblem(v, f, h, tspan, tstep, q₀, p₀; parameters = (k = 0.5,))

Geometric Equation Problem for Hamiltonian Ordinary Differential Equation (HODE)

 with vector fields
   v = v
   f = f

 Hamiltonian: H = h

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0])

 Parameters: 
   (k = 0.5,)

Implicit Ordinary Differential Equations (IODEs)

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

Example: Harmonic Oscillator

As an example, let us consider the harmonic oscillator. In implicit form, its equations are given as follows,

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

Here, v acts as a Lagrange multiplier that enforces the "constraint" $p(t) = v(t)$.

In order to create an IODEProblem for the harmonic oscillator, we thus need to write the following code:

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

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

function g(f, t, q, v, params)
    p[1] = - params.k * q[1]
end

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [0.5]
p₀ = [0.0]

prob = IODEProblem(p, f, tspan, tstep, q₀, p₀; parameters = (k = 0.5,))

Geometric Equation Problem for Implicit Ordinary Differential Equation (IODE)

 with vector fields
   v = p
   f = f
   g = _iode_default_g

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0], v = [0.0])

 Parameters: 
   (k = 0.5,)

The energy of the harmonic oscillator is preserved, so we can add it as an invariant,

energy(t, q, v, params) = v[1]^2 / 2 + params.k * q[1]^2 / 2

prob = IODEProblem(p, f, tspan, tstep, q₀, p₀; parameters = (k = 0.5,), invariants = (h=energy,))

Geometric Equation Problem for Implicit Ordinary Differential Equation (IODE)

 with vector fields
   v = p
   f = f
   g = _iode_default_g

 Invariants: 
   (h = Main.energy,)

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0], v = [0.0])

 Parameters: 
   (k = 0.5,)

Lagrangian Ordinary Differential Equations (LODEs)

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

Example: Harmonic Oscillator

As an example, let us consider the harmonic oscillator. Its Lagrangian is given by

\[L(q, \dot{q}) = \frac{\dot{q}^2}{2} - k \, \frac{q^2}{2} ,\]

so that the Euler-Lagrange equations

\[\frac{d}{dt} \frac{\partial L}{\partial \dot{q}} = \frac{\partial L}{\partial q} ,\]

become

\[\ddot{q} (t) = - k q(t) .\]

Most integrators for Lagrangian systems do not solve this second order system (semi-spray form), but instead use a reformulation as an implicit ordinary differential equation. This formulation can most easily be obtained from a Hamilton-Pontryagin principle

\[\delta \int \limits_{t_0}^{t_1} \big[ L(q, v) + \left< p , \dot{q} - v \right> \big] = 0 ,\]

as follows,

\[\begin{aligned} \dot{q} (t) &= v(t) , \\ p(t) &= \frac{\partial L}{\partial v} (q(t),v(t)) = v(t) , \\ \dot{p} (t) &= \frac{\partial L}{\partial q} (q(t),v(t)) = - k q(t) . \end{aligned}\]

Here, v acts as a Lagrange multiplier that enforces the "constraint" $p(t) = \partial L / \partial v$.

In order to create an LODEProblem for the harmonic oscillator, we thus need to write the following code:

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

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

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

function l(t, q, v, params)
    v[1]^2 / 2 - params.k * q[1]^2 / 2
end

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [0.5]
p₀ = [0.0]

prob = LODEProblem(p, f, ω, l, tspan, tstep, q₀, p₀; parameters = (k = 0.5,))

Geometric Equation Problem for Lagrangian Ordinary Differential Equation (LODE)

 with vector fields
   ϑ = p
   f = f
   g = _lode_default_g

 Lagrangian: L = l

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0], v = [0.0])

 Parameters: 
   (k = 0.5,)

The energy of the harmonic oscillator is preserved, so we can add it as an invariant,

energy(t, q, v, params) = v[1]^2 / 2 + params.k * q[1]^2 / 2

prob = LODEProblem(p, f, ω, l, tspan, tstep, q₀, p₀; parameters = (k = 0.5,), invariants = (h=energy,))

Geometric Equation Problem for Lagrangian Ordinary Differential Equation (LODE)

 with vector fields
   ϑ = p
   f = f
   g = _lode_default_g

 Lagrangian: L = l

 Invariants: 
   (h = Main.energy,)

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0], v = [0.0])

 Parameters: 
   (k = 0.5,)

Differential Algebraic Equation (DAE)

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

Example: Harmonic Oscillator

As an example we consider the harmonic oscillator, with an additional constraint that enforces energy conservation. While the system itself is energy conserving, most integrators do not respect this property. A possible way of remedying this flaw is to explicitly add energy conservation as an algebraic constraint.

The dynamical equations are given by

\[\dot{q} (t) = \begin{pmatrix} 0 & 1 \\ -k & 0 \\ \end{pmatrix} q(t) + \nabla \phi (q(t)) \lambda , \qquad \phi (q) = \frac{q_2^2}{2} + k \, \frac{q_1^2}{2} \qquad q \in \mathbb{R}^{2} , \qquad \lambda \in \mathbb{R}^{1} .\]

In order to create an DAEProblem for the harmonic oscillator including the projection on the constant energy manifold, we need to write the following code:

hamiltonian(t, q, params) = q[2]^2 / 2 + params.k * q[1]^2 / 2

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

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

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

tspan = (0.0, 1.0)
tstep = 0.1
q₀ = [1., 1.]
λ₀ = [0.]
params = (k=0.5,)

prob = DAEProblem(v, u, ϕ, tspan, tstep, q₀, λ₀; parameters = params)

Geometric Equation Problem for Differential Algebraic Equation (DAE)

 with vector fields
   v = v
   u = u
   ū = nothing

 and constraints
   ϕ = ϕ
   ψ = nothing

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [1.0, 1.0], λ = [0.0], μ = [0.0])

 Parameters: 
   (k = 0.5,)

Partitioned Differential Algebraic Equation (PDAE)

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

Example: Pendulum

As an example we consider the pendulum in cartesian coordinates, with a constraint that enforces the solution to respect the length of the pendulum.

The dynamical equations are given by

\[\begin{aligned} \dot{q} (t) &= p (t) , & q &\in \mathbb{R}^{2}, \\ \dot{p} (t) &= \begin{pmatrix} - \lambda q_1 (t) \\ 1 - \lambda q_2 (t) \\ \end{pmatrix} , & p &\in \mathbb{R}^{2} , \\ \phi (q) &= q^2 - l^2 , & \lambda &\in \mathbb{R}^{1} . \end{aligned}\]

In order to create an PDAEProblem for the pendulum including the projection on the constant length constraint, we need to write the following code:

function v(v, t, q, p, params)
    v .= p
end

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

function u(u, t, q, p, λ, params)
    u .= 0
end

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

function ϕ(ϕ, t, q, p, params)
    ϕ[1] = q[1]^2 + q[2]^2 - params.l^2
end

tspan = (0.0, 1.0)
tstep = 0.1
params = (l=0.5,)

q₀ = [params.l, 0.0]
p₀ = [0.0, 0.1]
λ₀ = [0.0]

prob = PDAEProblem(v, f, u, g, ϕ, tspan, tstep, q₀, p₀, λ₀; parameters = params)

Geometric Equation Problem for Partitioned Differential Algebraic Equation (PDAE)

 with vector fields
   v = v
   u = u
   ū = nothing
   f = f
   g = g
   ḡ = nothing

 and constraints
   ϕ = ϕ
   ψ = nothing

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0], p = [0.0, 0.1], λ = [0.0], μ = [0.0])

 Parameters: 
   (l = 0.5,)

Hamiltonian Differential Algebraic Equation (HDAE)

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

Example: Pendulum

As an example we consider the pendulum in cartesian coordinates, with a constraint that enforces the solution to respect the length of the pendulum.

The Hamiltonian is given by

\[H(q,p) = \frac{1}{2} p^2 + l - q_2 ,\]

and dynamical equations read

\[\begin{aligned} \dot{q} (t) &= p (t) , & q &\in \mathbb{R}^{2}, \\ \dot{p} (t) &= \begin{pmatrix} - \lambda q_1 (t) \\ 1 - \lambda q_2 (t) \\ \end{pmatrix} , & p &\in \mathbb{R}^{2} , \\ \phi (q) &= q^2 - l^2 , & \lambda &\in \mathbb{R}^{1} . \end{aligned}\]

In order to create an HDAEProblem for the pendulum including the projection on the constant length constraint, we need to write the following code:

function v(v, t, q, p, params)
    v .= p
end

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

function u(u, t, q, p, λ, params)
    u .= 0
end

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

function ϕ(ϕ, t, q, p, params)
    ϕ[1] = q[1]^2 + q[2]^2 - params.l^2
end

function h(t, q, p, params)
    return (p[1]^2 + p[2]^2)/2 + params.l - q[2]
end

tspan = (0.0, 1.0)
tstep = 0.1
params = (l=0.5,)

q₀ = [params.l, 0.0]
p₀ = [0.0, 0.1]
λ₀ = [0.0]

prob = HDAEProblem(v, f, u, g, ϕ, h, tspan, tstep, q₀, p₀, λ₀; parameters = params)

Geometric Equation Problem for Hamiltonian Differential Algebraic Equation (HDAE)

 with vector fields
   v = v
   f = f
   u = u
   g = g
   ū = nothing
   ḡ = nothing

 and constraints
   ϕ = ϕ
   ψ = nothing

 Hamiltonian: H = h

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0], p = [0.0, 0.1], λ = [0.0], μ = [0.0])

 Parameters: 
   (l = 0.5,)

Implicit Differential Algebraic Equation (IDAE)

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

Example: Pendulum

As an example we consider the pendulum in cartesian coordinates, with a constraint that enforces the solution to respect the length of the pendulum.

The implicit equations are given by

\[\begin{aligned} \dot{q} (t) &= v(t) , & q &\in \mathbb{R}^{2}, \\ \dot{p} (t) &= \begin{pmatrix} 0 \\ 1 \\ \end{pmatrix} - \lambda(t) q(t) , & v &\in \mathbb{R}^{2} , \\ p(t) &= v(t) , & p &\in \mathbb{R}^{2} , \\ \phi (q) &= q^2 - l^2 , & \lambda &\in \mathbb{R}^{1} . \end{aligned}\]

In order to create an IDAEProblem for the pendulum including the projection on the constant length constraint, we need to write the following code:

function p(p, t, q, v, params)
    p .= v
end

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

function u(u, t, q, v, p, λ, params)
    u .= 0
end

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

function ϕ(ϕ, t, q, v, p, params)
    ϕ[1] = q[1]^2 + q[2]^2 - params.l^2
end

tspan = (0.0, 1.0)
tstep = 0.1
params = (l=0.5,)

q₀ = [params.l, 0.0]
p₀ = [0.0, 0.1]
λ₀ = [0.0]

prob = IDAEProblem(p, f, u, g, ϕ, tspan, tstep, q₀, p₀, λ₀; parameters = params)

Geometric Equation Problem for Implicit Differential Algebraic Equation (IDAE)

 with vector fields
   ϑ = p
   f = f
   u = u
   g = g
   ū = nothing
   ḡ = nothing

 and constraints
   ϕ = ϕ
   ψ = nothing

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0], p = [0.0, 0.1], v = [0.0, 0.0], λ = [0.0], μ = [0.0])

 Parameters: 
   (l = 0.5,)

Lagrangian Differential Algebraic Equation (LDAE)

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

Example: Pendulum

As an example we consider the pendulum in cartesian coordinates, with a constraint that enforces the solution to respect the length of the pendulum.

The Hamilton-Pontryagin principle for this system is given by

\[\delta \int \limits_{t_0}^{t_1} \big[ L(q, v) + \left< p , \dot{q} - v \right> + \lambda \, \phi (q) \big] = 0 ,\]

with Lagrangian

\[L(q,v) = \frac{1}{2} v^2 - (l - q_2) ,\]

and constraint

\[\phi(q) = q^2 - l^2 .\]

The resulting implicit equations read

\[\begin{aligned} \dot{q} (t) &= v(t) , & q &\in \mathbb{R}^{2}, \\ \dot{p} (t) &= \begin{pmatrix} 0 \\ 1 \\ \end{pmatrix} - \lambda(t) q(t) , & v &\in \mathbb{R}^{2} , \\ p(t) &= v(t) , & p &\in \mathbb{R}^{2} , \\ \phi (q) &= q^2 - l^2 , & \lambda &\in \mathbb{R}^{1} . \end{aligned}\]

In order to create an LDAEProblem for the pendulum including the projection on the constant length constraint, we need to write the following code:

function p(p, t, q, v, params)
    p .= v
end

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

function u(u, t, q, v, p, λ, params)
    u .= 0
end

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

function ϕ(ϕ, t, q, v, p, params)
    ϕ[1] = q[1]^2 + q[2]^2 - params.l^2
end

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

function l(t, q, v, params)
    return (v[1]^2 + v[2]^2)/2 - (params.l - q[2])
end

tspan = (0.0, 1.0)
tstep = 0.1
params = (l=0.5,)

q₀ = [params.l, 0.0]
p₀ = [0.0, 0.1]
λ₀ = [0.0]

prob = LDAEProblem(p, f, u, g, ϕ, ω, l, tspan, tstep, q₀, p₀, λ₀; parameters = params)

Geometric Equation Problem for Lagrangian Differential Algebraic Equation (LDAE)

 with vector fields
   ϑ = p
   f = f
   u = u
   g = g
   ū = nothing
   ḡ = nothing

 and constraints
   ϕ = ϕ
   ψ = nothing

 Lagrangian: L = l

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.1 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0], p = [0.0, 0.1], v = [0.0, 0.0], λ = [0.0], μ = [0.0])

 Parameters: 
   (l = 0.5,)

Stochastic Differential Equations (SDEs)

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

Example: Kubo Oscillator

function v(t, q, v, params)
    v[1] = + params.λ * q[2]
    v[2] = - params.λ * q[1]
end

function B(t, q, B, params)
    for j in axes(B, 2)
        B[1,j] = + params.ν * q[2]
        B[2,j] = - params.ν * q[1]
    end
end

tspan = (0.0, 1.0)
tstep = 0.01
q₀ = [0.5, 0.0]

prob = SDEProblem(v, B, WienerProcess(), tspan, tstep, q₀; parameters = (λ=2.0, μ=1.0))

Geometric Equation Problem for Stratonovich Stochastic Differential Equation (SDE)

 with vector field
   v = v

 and diffusion matrix
   B = B

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.01 

 Initial conditions: 
   (t = 0.0, q = [0.5, 0.0])

 Parameters: 
   (λ = 2.0, μ = 1.0)

Partitioned Stochastic Differential Equation (PSDE)

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

Example: Kubo Oscillator

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

function f(f, t, q, p, params)
    f[1] = - params.λ * q[1]
end

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

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

tspan = (0.0, 1.0)
tstep = 0.01
q₀ = [0.5]
p₀ = [0.0]

prob = PSDEProblem(v, f, B, G, WienerProcess(), tspan, tstep, q₀, p₀; parameters = (λ=2.0, μ=1.0))

Geometric Equation Problem for Stratonovich Partitioned Stochastic Differential Equation (PSDE)

 with vector fields
   v = v
   f = f

 and diffusion matrices
   B = B
   G = G

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.01 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0])

 Parameters: 
   (λ = 2.0, μ = 1.0)

Split Partitioned Stochastic Differential Equation (SPSDE)

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

Example: Kubo Oscillator

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

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

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

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

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

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

tspan = (0.0, 1.0)
tstep = 0.01
q₀ = [0.5]
p₀ = [0.0]

prob = SPSDEProblem(v, f1, f2, B, G1, G2, WienerProcess(), tspan, tstep, q₀, p₀; parameters = (λ=2.0, μ=1.0))

Geometric Equation Problem for Stratonovich Split Partitioned Stochastic Differential Equation (PSDE)

 with vector fields
   v  = v
   f1 = f1
   f2 = f2

 and diffusion matrices
   B  = B
   G1 = G1
   G2 = G2

 Invariants: 
   NullInvariants()

 Timespan: (0.0, 1.0) 
 Timestep: 0.01 

 Initial conditions: 
   (t = 0.0, q = [0.5], p = [0.0])

 Parameters: 
   (λ = 2.0, μ = 1.0)