Unexpected behaviour from Continuous Callback in DifferentialEquations.jl simulating Lotka-Volterra Dynamics

[colebrookson]

Hi there,

So I’m using a ContinuousCallback method to simulate a simple Lotka-Volterra consumer-resource system of equations:

\frac{dC}{dt} = (\alpha - \beta R)C

\frac{dR}{dt} = (-\gamma + \delta C)R

My goal is to simply create a callback such that when the resource R goes below a certain value (in my example, 60), the callback function takes the integrator to zero.

Using the bouncing ball example as a basis, I have formatted the code as following:

function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

function zero_condition!(u,t,integrator) # Event when event_f(u,t) <= 60
  u[1] <= 60
end

function zero_affect!(integrator)
  integrator.u[1] = 0
end

zero_event = ContinuousCallback(zero_condition!,zero_affect!)

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,50.0)
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = zero_event)
plot(sol)

The code runs, and produces the following solution:

print(sol)
u: [[100.0, 100.0], [80.80845686345101, 114.00221332347269], [59.438166455464156, 126.7071869446583], [38.37959894566831, 134.27941573269763], [18.843731217381873, 130.52474821796022], [11.656653966686788, 121.35951343513864], [6.436587909054811, 105.673103443787], [3.874110828007275, 89.5763973309722], [2.2403782888041226, 69.74988819115794], [1.4937160732283663, 52.82058943966623], [1.1299759092760473, 38.650593036020666], [0.9615473179089438, 26.694047651872364], [0.9337419507717553, 17.161970095532308], [1.054390144771701, 9.912042531836985], [1.347009098584734, 5.531988127100285], [1.8952712402506138, 2.955913455350605], [2.8618179586434507, 1.5337792854669305], [4.579730793457142, 0.7807746687626642], [7.765207862902192, 0.3974028644034909], [14.716051089510543, 0.20591489694526696], [26.934456181374337, 0.14530645535833772], [51.10710701790373, 0.17148599156666613], [97.11195013965043, 0.5789527461182298], [0.0, 0.5789527461182298], [0.0, 0.2140524219697758], [0.0, 0.07916734689187103], [0.0, 0.030792621980557648], [0.0, 0.012118725928445456], [0.0, 0.004872909916774048], [0.0, 0.0019705425705026595], [0.0, 0.0013744561046710527]]

Where I am confused is that the callback is doing something, I can see it take the value of u[1] to zero, but it’s not at the value 60.

I’m sure this is a fundamental misunderstanding on my part of how to handle the ContinuousCallback, but I can’t find an answer to this in the docs, and would appreciate any help!

[tomaklutfu]

@colebrookson

function zero_condition!(u,t,integrator) # Event when event_f(u,t) <= 60
  u[1] <= 60
end

You need a zero crossing function in continuous callback. It is triggered by sign change(-1 to 1, up-crossing or 1 to -1, down-crossing) of the return value. Try below function

function zero_condition!(u,t,integrator) # Event when event_f(u,t) <= 60
  60-u[1] #Detection for down-crossing change
end

Note: Sorry, I made a mistake in early version. Updated the function as suggested in the next comment.

[hzgzh]

try

function zero_condition!(u,t,integrator) # Event when event_f(u,t) <= 60
  u[1]  - 60
end

zero_event = ContinuousCallback(zero_condition!,zero_affect!, nothing)

affect! is upcross, affect_neg! is downcross, so need set affect_neg! to nothing

[colebrookson]

Follow-up question:

If I wanted to check if both u[1] and u[2] go below 60 what would be the best way to do that, specifically that would scale if there were a large number of equations (i.e. u had 50 elements)?

I can do it this explicit way, wherein I simply write it out twice:

function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

function u1_zero_condition!(u,t,integrator) # Event when event_f(u,t) == 0
  60-u[1]
end
function u2_zero_condition!(u,t,integrator) # Event when event_f(u,t) == 0
  60-u[2]
end

function u1_zero_affect!(integrator)
  integrator.u[1] = 0
end
function u2_zero_affect!(integrator)
  integrator.u[2] = 0
end

u1_zero_event = ContinuousCallback(u1_zero_condition!,u1_zero_affect!)
u2_zero_event = ContinuousCallback(u2_zero_condition!,u2_zero_affect!)

cb_set = CallbackSet(u1_zero_event, u2_zero_event)

With the same run section:

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = cb_set)
plot(sol)
@show sol[end]
print(sol)
minimum(sol)

Which produces this:

So I know I’m getting exactly what I want. However, writing it out manually for a large system (i.e. 50) seems inefficient. Is there a standard way of doing this? Or alternatively, could I somehow have one condition! and one affect! function check all of u[i] where there are i=50 elements?

Any thoughts much appreciated!

[tomaklutfu]

You can use VectorContinuousCallback. The docs are here

[colebrookson]

Sorry to come back on this, using VectorContinuousCallback but couldn’t find any examples of having multiple forms of callbacks, wonder if you have insight on this. So I’m trying to perform the above task (taking the function to zero if it crosses some value), but ALSO, change the parameter values at some specific times. I could previously achieve that with DiscreteCallbacks but now I’m trying to get it to play nice with the vector of callbacks.

I tried to implement both through the out[i] component of the condition, but could only get it to do either taking the function to zero OR change the parameters, but not both.

Code here:

function dynamics!(du,u,p,t)
  α, β, γ, δ = p
  du[1] = (α - β * u[2]) * u[1]
  du[2] = (-γ + δ * u[1]) * u[2]
end

# set constant stop values
const stop1 = [1.0]
const stop2 = [3.0]

function condition(out,u,t,integrator)
  out[1] = t in stop1
  out[2] = t in stop2
  out[3] = 20-u[1]
  out[4] = 20-u[2]
end

function affect!(integrator, idx)
  if idx == 1
    integrator.p[4] += 4
  end
  if idx == 2
    integrator.p[4] -= 4
  end
  if idx == 3
    integrator.u[1] = 0
  end
  if idx == 4
    integrator.u[2] = 0
  end
end

cb = VectorContinuousCallback(condition,affect!,4,save_positions=(true,true))

u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
const tstop = [1.;3.]
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = cb, tstops=tstop)
plot(sol)

which gives me this:

However, if i simply remove the second two indices of out[i] like this,

function condition(out,u,t,integrator)
  out[1] = t in stop1
  out[2] = t in stop2
end

but change nothing else, I see the parameter change take effect, and it looks like this:

How can I get BOTH of these callbacks to work together in the callback set?

Again, any help much appreciated and thanks for all your help so far!

[tomaklutfu]

Your first two condition is boolean(true, false). They also should be sign changingin values, afaik. Something like

out[1] = 1.0 - t
out[3] = 3.0 - t

Instead, you can combine these 2 to a DiscreteCallback, while its condition checks and returns true in case of any discrete event and affect reacts appropriate action depending on time. In this case, you need to add tstops to integrator while making them into VectorContinuousCallback catches regardless of tstops.