Compute Invariants With DifferentialEquations.jl

Before reading this, it's worth familiarising yourself with the Ensemble Simulations section of the DifferentialEquations.jl docs.

To compute invariants with DifferentialEquations.jl, we'll need to create a PIEnsembleProblem. This works as follows

PIEnsembleProblem(init, prob, pinv;
    prob_func = (prob,i,repeat)->prob,
    output_func = (sol,i) -> (sol,false),
    reduction = (u,data,I)->(append!(u,data),false),
    u_init = [], safetycopy = false
)

init is the phasespace curve or surface parameterisation, which determines the initial conditions; prob is a problem from DifferentialEquations.jl, which specifies the time span and differential equation; pinv is an AbstractPoincareInvariant; prob_func allows the user to remake the problem for each trajectory in the ensemble; all further arguments are exactly like the EnsembleProblem type from DifferentialEquations.jl. You need not add a prob_func of your own to have the correct initial conditions. These will be entered automaticaly according ot the given setup objet and parameterisation. However, by default, the type of the initial condition is Vector{T}. If you want to use an ArrayPartion or a StaticArray or other type to represent a point in phase space, you should specify a prob_func to do the conversion.

All in all, this can look like so:

using OrdinaryDiffEq
using RecursiveArrayTools: ArrayPartition
using PoincareInvariants

dt = 0.1
prob = SecondOrderODEProblem((p, θ, params, t) -> [-sin(θ[1])], 0.0, 0.0, (0.0, 2.0))
pf(prob, i, repeat) = remake(prob; u0 = ArrayPartition((prob.u0[1:1], prob.u0[2:2])))

pi1 = CanonicalFirstPI{Float64, 2}(1_000)
ens_prob = PIEnsembleProblem(ϕ -> (sinpi(2ϕ), 3 * cospi(2ϕ)), prob, pi1; prob_func=pf)

Once, we have our PIEnsembleProblem, we can call solve on it, just like the EnsembleProblem it wraps. The general pattern is

solve(prob::PIEnsembleProblem, alg, ensemblealg; kwargs...)

Most arguments are simply passed to the solve function called on the wrapped EnsembleProblem. The trajectories keyword argument is already set as getpointnum(pinv) and the adaptive keyword is set to false by default. Otherwise the behaviour is identical. For our example, we would have

sol = solve(ens_prob, SymplecticEuler(), EnsembleSerial(); dt=dt)

Finally, we can call compute! on the result, which returns a Vector{T} of invariant values for every saved timestep. You must make sure all trajectories save at the same time steps.

compute!(pi1, sol[, p])

The times and points are taken from the solution. An optional parameter p may be given, which is passed to the differential form.