Time-dependent events in ODE

question

#1

I recently started with Julia and wanted to implement one of my usual problems - implement time-depended events.

For now I have:

# Packages
using Plots
using DifferentialEquations

# Parameters
k21 = 0.14*24
k12 = 0.06*24
ke = 1.14*24
α = 0.5
β = 0.05
η = 0.477
μ = 0.218
k1 = 0.5
V1 = 6

# Time
maxtime = 10
tspan = (0.0, maxtime)

# Dose
stim = 100

# Initial conditions
x0 = [0 0 2e11 8e11]

# Model equations
function system(dy, y, p, t)
  dy[1] = k21*y[2] - (k12 + ke)*y[1]
  dy[2] = k12*y[1] - k21*y[2]
  dy[3] = (α - μ - η)*y[3] + β*y[4] - k1/V1*y[1]*y[3]
  dy[4] = μ*y[3] - β*y[4]
end

# Events
eventtimes = [2, 5]
function condition(y, t, integrator)
    t - eventtimes
end
function affect!(integrator)
    x0[1] = stim
end
cb = ContinuousCallback(condition, affect!)

# Solve
prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb)

# Plotting
plot(sol, layout = (2, 2))

But the output that is give is not correct. More specifically, the events are not taken into account and the initial condition doesn’t seems to be 0 for y1 but stim .

Any help would be greatly appreciated.


#2

Welcome to Julia!

I think you have to do to conditions

const eventtimes = [2, 5]
function condition1(y, t, integrator)
    t - eventtimes[1]
end
function condition1(y, t, integrator)
    t - eventtimes[2]
end

Two notes:

  • it is kinder to your helpers if you make your example as minimal as possible, so I don’t need to waste time an brain power to figure out that, for instance, all the parameters you define are irrelevant to the problem at hand.
  • check out the performance tips in the julia-docs. For instance you should declare your parameters const

#3

Thanks. But, I don’t think that this solves the problem… because I also tried it with just one eventtime a couldn’t get a proper output.


#4

You need to split into two callbacks (as your code errors on this), see http://docs.juliadiffeq.org/latest/features/callback_functions.html#CallbackSet-1 on how to combine them. Then the other error is

function affect!(integrator)
    integrator.u[1] = stim
end

as stated in the docs: http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-1:-Bouncing-Ball-1

So all in all, slightly modified:

using Plots
using DifferentialEquations

# Parameters
const k21 = 0.14*24
const k12 = 0.06*24
const ke = 1.14*24
const α = 0.5
const β = 0.05
const η = 0.477
const μ = 0.218
const k1 = 0.5
const V1 = 6

# Time
maxtime = 0.1
tspan = (0.0, maxtime)

# Dose
const stim = 100

# Initial conditions
x0 = [0 0 2e11 8e11]

# Model equations
function system(dy, y, p, t)
  dy[1] = k21*y[2] - (k12 + ke)*y[1]
  dy[2] = k12*y[1] - k21*y[2]
  dy[3] = (α - μ - η)*y[3] + β*y[4] - k1/V1*y[1]*y[3]
  dy[4] = μ*y[3] - β*y[4]
end

# Events
function condition(y, t, integrator)
    t - 0.05
end
function affect!(integrator)
    integrator.u[1] = stim
end
cb = ContinuousCallback(condition, affect!)

# Solve
prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb)

# Plotting
plot(sol, layout = (2, 2))

(But please, next time provide an example which is more minimal and runs faster.)


#5

Thanks! And I keep in my to post something faster and simpler next time.

But, regarding multiple eventtimes: having to write a condition for each event seem highly unpractical. I usually have the requirement to have 10 or even 100 events.


#6

The eventfunction needs to be zero at the events, so you could construct one which is zero at several points. Something like this should work:

const times = rand(100)
function condition(y, t, integrator)
    out = 1.0
    for tt in times
        out *= (t-tt)
    end
    out
end

#7

This is something on our mind: https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/331