6

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.

Rolfe Power
  • 138
  • 1
  • 5

1 Answers1

7

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))
Korsbo
  • 724
  • 5
  • 5
  • 1
    Great, thank you! Just on a side note, for this to run on my machine I had to add the parameter input to f, i.e. change the function definition to `f(du, u, p, t)`. But other than that it solves my problem perfectly, much appreciated. – Rolfe Power Apr 09 '19 at 17:28