Search code examples
juliadifferential-equations

Is there an idiomatic way to terminate integration after n callbacks in DifferentialEquations.jl


First of all, I am using the DifferentialEquations.jl library, which is fantastic! Anyway, my question is as follows:

Say for example, I have the following differential equation:

function f(du, u, t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -u[1] - 2 * u[1] * u[2]
    du[4] = -u[2] - u[1]^2 + u[2]^2
end

and I have a callback which is triggered every time the trajectory crosses the y axis:

function condition(u, t, integrator)
    u[2]
end

However, I need the integration to terminate after exactly 3 crossings. I am aware that the integration can be terminated by using the effect:

function affect!(integrator)
    terminate!(integrator)
end

but what is the proper way to allow for a counting of the number of callbacks until the termination criterion is met. Furthermore, is there a way to extend this methodology to n events with n different counts?

In my research I often need to look at Poincare maps and the first, second, third, etc. return to the map so I am in need of a framework that allows me to perform this counting termination. I am still new to Julia and so am trying to reinforce good idiomatic code early on. Any help is appreciated and please feel free to ask for clarification.


Solution

  • There is a userdata keyword argument to solve which can be useful for this. It allows you to pass objects to the integrator. These objects can be used in creative ways by the callback functions.

    If you pass userdata = Dict(:my_key=>:my_value) to solve, then you can access this from integrator.opts.userdata[:my_key].

    Here is a minimal example which controls how many times the callback is triggered before it actually terminates the simulation:

    function f(du, u, t)
        du[1] = sin(t)
    end
    function condition(u, t, integrator)
        u[1] 
    end
    
    function affect!(integrator)
        integrator.opts.userdata[:callback_count] +=1
        if integrator.opts.userdata[:callback_count] == integrator.opts.userdata[:max_count]
            terminate!(integrator)
        end
    end
    
    
    callback = ContinuousCallback(condition, affect!)
    
    u0 = [-1.]
    tspan = (0., 100.)
    
    prob = ODEProblem(f, u0, tspan)
    sol = solve(prob; callback=callback, userdata=Dict(:callback_count=>0, :max_count=>3))