Continuous analog of SavingCallback from DifferentialEquations.jl

Specific Domains Modelling & Simulations
1

So, the problem that I wanted to solve is the following. For a system that I study I wanted to record time and phase state when F(x, t) = 0, where F(x, t) is some non-linear function. If F(x, t) is linear, then there is such functional in DynamicalSystems.jl, but in my case I need something like F(x, t) = \sin (x_k - x_k^{0}). My system is 2\pi periodic in some of the variables and this effectively replaces a family of surfaces of section when trajectory wanders on torus and in \mathbb{R}^n. As far as I’ve understood, SavingCallback from library of callbacks does not allow me to do this since DiscreteCallback operates on different type of “event” (pardon my MATLAB/Python terminology :)). I tried to mimic its implementation to use in a ContinuousCallback, and come up with the following code:

struct ContSavedValues{tType, savevalType}
    t::Vector{tType}
    saveval::Vector{savevalType}
end

function ContSavedValues(::Type{tType}, ::Type{savevalType}) where {tType,savevalType}
    ContSavedValues{tType, savevalType}(Vector{tType}(), Vector{savevalType}())
end

function Base.show(io::IO, saved_values::ContSavedValues)
    tType = eltype(saved_values.t)
    savevalType = eltype(saved_values.saveval)
    print(io, "SavedValues{tType=", tType, ", savevalType=", savevalType, "}",
                "\nt:\n", saved_values.t, "\nsaveval:\n", saved_values.saveval)
end

mutable struct ContSavingAffect{tType, savevalType}
    saved_values::ContSavedValues{tType, savevalType}
end

function (affect!::ContSavingAffect)(integrator)
    push!(affect!.saved_values.t, integrator.t)
    push!(affect!.saved_values.saveval, Tuple(Float64(v) for v in integrator.u))
end

and use it like

observer_array = ContSavedValues(Float64, Tuple{Float64, Float64})
affectFunc! = ContSavingAffect(observer_array)
cb = ContinuousCallback(myCondition, affectFunc!, nothing)
sol = solve(myProblem, callback=cb)

The part that I’m worried the most is this:

function (affect!::ContSavingAffect)(integrator)
    push!(affect!.saved_values.t, integrator.t)
    push!(affect!.saved_values.saveval, Tuple(Float64(v) for v in integrator.u))
end

I’m assuming that when affectFunc! will finally be called, the state of integrator corresponds to intersection with surface of interest. It runs okay on few tests that I’ve made, but is there any caveat that I might have overlooked?

2

This is fine. The only downside of course is that, as a ContinuousCallback, it will pull back to the time of the event. This will decrease the effective step size and, if the step size is large enough that multiple crossings would have happened in the time span, it will pull all the way back to the first and recompute, step, etc. So in theory there are cases where this can have decreased performance over a well-written DiscreteCallback that checks the local interpolation and saves that every crossing in a timespan, though in practice this may not make much of a difference.

1 Like
3

Thank you for your response! Is my understanding correct that this problem could be solved just in terms of already implemented DiscreteCallback and CallbackSet, without writing my own version?

4

Not really. You’d basically have to write an entire callback which matches the DiffEq callback code for how it finds each crossing in an interval. It would be a tedious callback to write and only better in the specific case where you would have gone far beyond crossings. It’s not clear it’s worth the work.

5

Aha, thank you for clarification! So, if I understood you correctly, I would have to repeat the code for affect!::SavingAffect in /src/saving.jl? Or it would be an entirely different and more tedious version of that callback?

6

no, you would have to essentially do a custom version of the ContinuousCallback code inside of the DiscreteCallback that ignores the pullback portion.

7

Thank you for your answer! I don’t know if it is suitable to ask another closely related question here, but I’ll try. Is event localization actively affects trajectory computation in DifferentialEquations.jl? We’ve encountered that computing trajectory with ContinuousCallback (even with the one that seems to not affect integrator in any way, like my version of SavingCallback) gives a significantly different trajectory in comparison with the case of no callbacks. Is it a bug or it works as intended?

8

That’s exactly what I’m pointing out. It’s working as intended. The ContinuousCallback version will be computed slightly more accurately because its dt will keep getting pulled back more than the adaptivity wants, in order to hit the crossing points exactly. That will make it more accurate, since dt = min(chosen_dt,event_dt), so it only decreases step sizes. The issue is that if you didn’t want more accuracy, you get it anyways, which is what I described the DiscreteCallback version for but is probably overkill in most situations.

1 Like