Ensemble simulations with "paired" trajectories

Specific Domains Modelling & Simulations
1

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 reductions should look like in order to achieve this kind of result, so I thought I would ask here.

2

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

https://juliadynamics.github.io/DynamicalSystems.jl/dev/advanced/#DynamicalSystemsBase.parallel_integrator

(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
3

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?

4

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 Redirecting

5

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
6

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