[feature request] Linear event handling

[samuela]

The event handling support in DifferentialEquations.jl is very robust and so far so great! However I have a situation in which my dynamics function is very expensive to compute. Therefore I'd like to minimize the number of root finding steps to the bare minimum. But the good news is that all my events are linear in the state of the system. So with a interpolant between t1 and t2 it should be possible to analytically compute the event time; no root finding necessary! My situation is analogous to the example here: https://diffeq.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball-with-multiple-walls.

Would it be possible to accomplish something like this with DifferentialEquations.jl? Either in DifferentialEquations.jl itself or with some custom hackery?

[ChrisRackauckas]

I think the best thing to do would be to use a DiscreteCallback and basically do event location yourself inside of it. The function that you need to do the effect correctly is https://diffeq.sciml.ai/stable/basics/integrator/#DiffEqBase.change_t_via_interpolation! which is how the ContinuousCallback moves back in time to the event location.

[ChrisRackauckas]

The event handling support in DifferentialEquations.jl is very robust and so far so great!

Unfortunately it's not 😢 and we're working on it at least. There are some weird bracketing issues in the rootfinding library where it can give roots that are outside of the original brackets... and other weird behaviors... so we're releasing our own rootfinding library relatively soon (in the next month or so).

[samuela]

Unfortunately it's not cry and we're working on it at least. There are some weird bracketing issues in the rootfinding library where it can give roots that are outside of the original brackets... and other weird behaviors... so we're releasing our own rootfinding library relatively soon (in the next month or so).

Oh :grimacing: well at least it looked like it was working for me...

[ChrisRackauckas]

Yeah it works for most cases. People do weird stuff with floating point numbers though... library writing teaches you crazy stuff about things like 1e-6 - 1f-6...

[samuela]

I think the best thing to do would be to use a DiscreteCallback and basically do event location yourself inside of it. The function that you need to do the effect correctly is https://diffeq.sciml.ai/stable/basics/integrator/#DiffEqBase.change_t_via_interpolation! which is how the ContinuousCallback moves back in time to the event location.

So DiffEqBase.change_t_via_interpolation! will walk time backward in the internal state of the solver and set the position at the new point in time according to the interpolant? That def seems handy.

Is there a convenient way to get the interpolant between the "current" time step when the DiscreteCallback fired, and the previous time step right before the DiscreteCallback fired?

[ChrisRackauckas]

integrator(t)

[samuela]

Is there a way to get t_prev, so that I know how far I should be backtracking?

[samuela]

Also is there a way to get the parameters of the interpolant? That will be necessary to do this in closed-form.

[ChrisRackauckas]

The straddle interval is [integrator.t,integrator.t+integrator.dt].

Also is there a way to get the parameters of the interpolant? That will be necessary to do this in closed-form.

A lot of the methods have a specific interpolation method, each defined here: https://github.com/SciML/OrdinaryDiffEq.jl/blob/master/src/dense/interpolants.jl

[samuela]

The straddle interval is [integrator.t,integrator.t+integrator.dt].

I'm trying this with a simple example but it's not quite operating the way I'd expect:

cb = DiscreteCallback((u, t, integrator) -> u[1] <= 0, (integrator) -> begin
    @show integrator.t
    @show integrator.t + integrator.dt
end)
solve(ODEProblem((u, _, _) -> [u[2], 0], [0.75, -1], (0.0, 1.0)), callback = cb)

shows

integrator.t = 1.0
integrator.t + integrator.dt = 2.0

But what I'm really after is an interval that straddles the zero-crossing at t = 0.75. Is there a way to get that interval?

Edit: Ok, I think integrator.tprev is more or less what I want.

[ChrisRackauckas]

However I have a situation in which my dynamics function is very expensive to compute. Therefore I'd like to minimize the number of root finding steps to the bare minimum.

Rootfinding does not require evaluating the dynamics function.

tprev is what you want. And I think someone is contributing a callback for this case.

[samuela]

Rootfinding does not require evaluating the dynamics function.

Ah, ok I was mistaken! Happy to hear that. I may be able to get away with using VectorContinuousCallback for the time being then.