Times of events with DifferentialEquations.jl

Specific Domains Numerics
1

I’m using DifferentialEquations.jl for the first time in a while (it’s changed a bit :slight_smile:) and I’m trying to work out how to get the times at which events happen with a ContinuousCallback function when integrating an ODE. Ultimately I’m trying to get a single period of a periodic orbit by defining a Poincare surface that it passes through once (twice) per period and extracting the time points at which that happens.

Is there any way to easily extract the time points that correspond to events found with a ContinuousCallback function? I’m could do it with a closure:

tt = Float64[]
cb = ContinuousCallback(poincare, (integrator) -> push!(tt, integrator.t))
sol = solve(prob, Vern6(), callback=cb, dtmax=0.02)

but I was wondering if there is an built in way (I can’t see one in the docs).

Thanks!

1 Like
2

Well, if you set save_everystep=false,save_start=false,save_end=false, then the only points which are saved correspond to the callbacks. If you use save_positions=(false,true) in the callback then only one time point per callback will be saved, and then sol.t will be what you want. That said, this setup would mean you’d save nothing else about the solution, which may or may not be what you’re looking for. If you want the solution to be undisturbed, then a closure is the way to go.

There seems to be quite a few people looking to do Poincure surfaces, so maybe adding a good callback for it to the callback library would be a good idea. A ContinuousCallback works (and is an easy way to get it working) but it restricts the stepping to only one crossing per step, which isn’t actually necessary here.

P.S. Tune into @Datseris 's DynamicalSystems.jl talk today!

1 Like
3

Thanks for the mention!

The function poincaresos of DynamicalSystems.jl computes the poincare section of abitrary continuous systems. The callback that does it is here: https://github.com/JuliaDynamics/ChaosTools.jl/blob/master/src/orbitdiagram.jl#L172

The function is “nice” in the sense that it can accept an arbitrary hyperplane as an argument and does everything else internally. I think you can either cite this function in your docs or just copy the code I’ve written. If you do the latter, please let me know so I can use your code, to avoid code duplication!

The current poincaresos function does not return the time of the crossings (it is available, just not returned). @dawbarton if you intend to use this function you can open up an feature request and the DynamicalSystems.jl repo and I’ll get to it!

4

It uses a ContinuousCallback too though, right? It would be best to try and come up with the DiscreteCallback implementation which doesn’t do unnecessary pullbacks and directly calls Roots.jl

5

Yeah it does.

What I mean is that you could of course use my code if it is of any use, but if you write anything better regardless, then please let me know!

6

Thanks both. I’d missed that DynamicalSystems.jl announcement - I’ll take a look.

Yes, I wanted the original solution intact rather than removing all output as I then do further calculations on it. I’ll keep using the closure for now. (For reference, in Matlab when you specify an event function it returns the original solution as sol.{x,y} and the solution at the events as a separate field sol.{xe,ye}.)

If I get chance I’ll take a look at using a DiscreteCallback - though I didn’t quite understand what you meant by the restriction to one crossing per time step for the ContinuousCallback.