# Ensemble simulations with "paired" trajectories

I’m interested in doing something roughly like the following:

1. Choose a (possibly random) initial condition u0
2. From u0, take a (possibly random) “perturbation” of u0, say, x0 = u0 + s
3. Integrate both trajectories, and evaluate some condition at the endpoint of the trajectories (e.g whether they approached the same attractor.)
4. Repeat 1-3 until my endpoint condition is satisfied N times, or maximum iterations is reached.
5. Terminate

I have a vague idea about how I might accomplish this serially (or possibly even with something like `Threads.@threads`), e.g:

``````using OrdinaryDiffEq

function eom(dx,x,p,t)
foo()
end

tspan = (0.0,100.0)

N_conditions = 0
N_total = 0
N_maxiters = 100

while N_conditions < N_target && N_total < N_maxiters
u0 = rand()
x0 = u0 + rand()

prob = ODEProblem(eom,u0,tspan)
sol_u0 = solve(prob)

remake(prob,u0 = x0)
sol_x0 = solve(prob)

if check_condition(sol_u0,sol_x0)
N_conditions += 1
N_total += 1
else
N_total += 1
end
end
``````

But this seems inefficient. I am wondering if this is possible to achieve as some kind of `EnsembleProblem` setup, with a `reduction` that will end the simulation once the target condition is reached. It’s not clear to me what the `prob_func` or `reduction`s should look like in order to achieve this kind of result, so I thought I would ask here.

The parallel integrator of DynamicalSystems.jl offers a framework to evolve multiple initial conditions in parallel at exactly same times

(it uses DiffEq under the hood)

If you do not care about having both solutions at exactly the same time at all times (e.g. only care that end points match), then `EnsembleProblem` is the way to go.

2 Likes

Just out of curiosity, are you doing this to calculate return times to an attractor, e.g. Poincare return times, to get the distribution of return times?

Thanks for your replies! I have previously played around with DynamicalSystems.jl, but the way I’ve defined my system (and the information I typically want out of the solution) probably lends itself a bit better to a DifferentialEquations setup.

If you do not care about having both solutions at exactly the same time at all times (e.g. only care that end points match), then `EnsembleProblem` is the way to go.

The only part of the solution I need for this problem is…well, none of it: the quantities I want to look at here are `N_conditions` and `N_total` (i.e, what are they when the simulation stops?) Hence, I’ll only need the endpoints as far as checking the condition.

Just out of curiosity, are you doing this to calculate return times to an attractor, e.g. Poincare return times, to get the distribution of return times?

No, I’m interested in a fractal dimension related to final-state uncertainty, as in https://doi.org/10.1016/0375-9601(83)90945-3

Oh very nice paper, thanks for sharing, I just went through this and I guess it would be nice to have a method that computes the fraction `f` in ChaosTools.jl, if you are interested! c.f. https://github.com/JuliaDynamics/ChaosTools.jl/issues/121

The way I imagine it is with either by initializing several integrators in parallel (what EnsembleProblem does) or by doing `reinit!` for every new initial condition (what I typically do in DynamicalSystems.jl, because you typically parallelize at a higher level, e.g. for different parameters).

The user can provide a `terminate` function that checks the state of the integrator and terminates integration on request (e.g. `terminate(integ) = norm(integ.u[3:4]) < 1e-3` a typical condition that terminates when velocity becomes too small). Once terminated, you check the final state of the integrator (always accessible, even if you don’t save any states) and check whether it matches the requested fixed point or not.

Regarding optimization, the function can take a series of sizes `ε` as input, because for `ε2 > ε1` you don’t have to recompute what initial conditions starting in ball `ε1` do; you have already done that in the `ε1` step.

1 Like

That would be very nice indeed! I’m happy to try to contribute once I get something working for my own case.