Event functions

[tomasaschan]

For some problems, it is interesting to run the solver until something happens - specified by (part of) the solution at that point. In those cases, it can be hard to know on beforehand exactly when this happens, so it's difficult to give an upper limit to the solution domain.

For example, think about a program modelling a bouncing ball: it's just solving F = ma for a free-falling object, except every time the ball hits the floor, we want to change sign of the y-component of its velocity. Without event functions, one has to solve for one parabola at a time, obtain the final components of the velocity, change the sign and use as initial values for the next - with them, one can let the solver do the work.

It could also be nice to be able to stop the simulation when some condition is met, rather than after a final time - for a system that you expect to stabilize but you don't know how long it's going to take, you can then stop when the rate of change in the solution is small enough.

Basically, the solvers should support some way of saying every time X happens, do Y where Y can be anything from changing some part of the problem to stopping the solver. As noted in #7 it could also be interesting to support only emitting solution points when some condition is met, which could also be specified by an event function.

What is a good way of implementing this in the current solvers? What should the user interface look like?

[stevengj]

Changing the problem as a function of the solution can already be accomplished in the ODE callback function, unless I am misunderstanding you.

To halt the solver, maybe the ODE callback could throw some special exception?

[tomasaschan]

Yeah - you're probably right. Changing the problem along the way can be accomplished without this.

Throwing an exception to halt the solver seems like a code smell to me - it's not really an error as such. I'd rather have the ODE callback return some special symbol, but that would be harmful for performance since it's return type wouldn't be fixed... It also wouldn't work if we want the event function to specify when to output solution points, but not halt the solver - that depends on the ODE callback returning as usual, but something else gets evaluated too.

Maybe we could have a keyword argument event that takes a function g(t, y) which returns either :halt, :out (to output solution, but probably with a better name) or :donothing (also with a better name), that gets evaluated once at the end of each iteration?

[stevengj]

Premature termination seems like an exception to me...

I'd prefer not to have to write two callback functions when one will do, especially since the termination criteria may depend on internal state of the dx/dt callback.

[tomasaschan]

"Premature" is probably a bad term for what I'm trying to get at - I'm rather looking for cases when it is not possible to state the end time in advance.

[acroy]

I would also say that exceptions provide a clean way to terminate the solver. However, I am not sure that triggering events in the dx/dt callback is a good idea. Firstly, because one only gets local information (x and dx/dt at the current time) and secondly, because one does not know in which context the function was called. For instance, one cannot know whether the time-step will actually be accepted or not and why the solver actually called the function. The locality might be a problem if one wants to check if a given function of x and/or dx/dt is passing through zero during the integration step (this is basically how events work in MATLAB).

[ivarne]

Maybe not directly on topic, but I have sometimes wanted to have some means of telling the solver to reduce the step size if some of the integrated properties is out of bounds. If i know that all values in the x vector will be positive it would be nice if I could throw an exception (ReduceStepSizeException?) if the algorithm overshoots and tries with a negative value.

[sglyon]

Any update on this? Does anyone have a working example of using callbacks to emulate the events functionality of Matlab?

[mauro3]

I don't think so. The new RK solvers have dense output (although only 3rd order polynomials). This could then be used to detect events.

[sglyon]

I needed this, so I implemented a routine that uses any integrator in ODE.jl. The routine is in this gist

I'm quite positive this is one of the most inefficient ways of doing it, but it works for my needs right now and I thought I'd share.

[acroy]

Since it just came up in #33: At some point I suggested a simple way to have a progress callback for Sundials, which might in principle also work for our purposes?

[geandersoncarvalho]

I am new in julia programming language.
I am trying to solve an ODE but I need to set a stop condition. spencerlyon2 I try to use your routine as a tentative to set the this, but when I ran your routine even with your example I get the following error messsage
ERROR: LoadError: Zero time span
in hinit at /home/carvalho/.julia/v0.4/ODE/src/ODE.jl:67
while loading /home/carvalho/Downloads/conditionstopJULIA/teste2.jl, in expression starting on line 125
where line 125 was
X, Y, xe, ye, ie = ode_events(yp, y0, xspan, [event], ode45)
I don't understand what can I have to do to fix this problem

[sglyon]

Hi @geandersoncarvalho I'm sorry but I wasn't ever able to fully finish the events code I sent a link to here so I can't make any guarantees that it will work.

[ChrisRackauckas]

Thank you for this thread. I used a lot of ideas from here. I just finished implementing callbacks in DifferentialEquations.jl and made sure the easy interface addressed at least all of the concerns here (it does a lot more, but at least everything here). For example, the easy interface signature is (in full, only event_f and apply_event! don't have defaults)

callback = @ode_callback begin
  @ode_event event_f apply_event! rootfind_event_loc interp_points terminate_on_event Δt_safety
end

1. The option terminate_on_event doesn't use exceptions because those have many complications (including causing performance issues if you're solving a lot of ODEs). Instead it makes the final timepoint equal to the current timepoint, so the loop will simply end itself and undergo any postamble
2. This form of termination allows one to solve on a timespan [0;Inf] for cases where the end is not known.
3. This callback is only called when the step is accepted, addressing @acroy's concern.
4. This type of event has Δt_safety which will be a multiplier to Δt which can be used to do things like @ivarne was talking about. But using the event to directly go to the time of the event is already automatic if rootfind_event_loc is true (it will rootfind to the event), so that would actually handle @ivarne 's case better.
5. A lot of the DifferentialEquations.jl algorithms have higher than 3rd order interpolations, with algorithm-specific interpolations up to 9th order which are used for the rootfinding and event detection. This is why a DSL with macros was created because there's a lot of "algorithm specific code".
6. This implementation is pretty well optimized, and since the DifferentialEquations.jl algorithms benchmark at around 10x faster than the ODE.jl algorithms, this should address your performance concerns @spencerlyon2. I say "pretty well" since there is one place where there's allocations, but it has a pretty trivial fix: in-place odeinterpolant DifferentialEquations.jl#66. Then the implementation will be allocation-free and very swift. On all of the problems I've checked the callbacks inline as well.
7. I didn't add the progress callback as @acroy mentions since I already do that somewhere else. I could move that in here, but it just seems unnecessary. Someone could add their own progressbar in the callback function though.
8. One other cool feature I'd like to mention is that you can change the size of the ODE along the way by using this. I think it's cool.

Once again, thank you for this thread of ideas. It was a helpful source of "user-concerns" throughout the design/implementation.

[rtburns-jpl]

Are there any plans to implement callbacks for ODE.jl? I'd love to use @ChrisRackauckas's new callback functionality with the higher-order ODE algorithms here.

[ChrisRackauckas]

No plans. Just use https://github.com/JuliaDiffEq/DifferentialEquations.jl . The higher order methods there are more efficient anyways.