# Time-dependent events in ODE

**URL:** <https://discourse.julialang.org/t/time-dependent-events-in-ode/14951>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [September 14, 2018, 11:32am UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951 "2018-09-14T11:32:17Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![hellopascal](https://avatars.discourse-cdn.com/v4/letter/h/a8b319/32.png) [@hellopascal](https://discourse.julialang.org/u/hellopascal)\
**Post date:** [September 14, 2018, 11:32am UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/1 "2018-09-14T11:32:17Z")

</div>

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

For now I have:

```julia
# 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.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [September 14, 2018, 11:40am UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/2 "2018-09-14T11:40:36Z")

</div>

Welcome to Julia!

I think you have to do to conditions

```julia
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`

---

<div class="post-metadata">

**Author:** ![hellopascal](https://avatars.discourse-cdn.com/v4/letter/h/a8b319/32.png) [@hellopascal](https://discourse.julialang.org/u/hellopascal)\
**Post date:** [September 14, 2018, 2:04pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/3 "2018-09-14T14:04:50Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [September 14, 2018, 2:21pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/4 "2018-09-14T14:21:59Z")

</div>

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](http://docs.juliadiffeq.org/latest/features/callback_functions.html#CallbackSet-1) on how to combine them. Then the other error is

```julia
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](http://docs.juliadiffeq.org/latest/features/callback_functions.html#Example-1:-Bouncing-Ball-1)

So all in all, slightly modified:

```julia
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.)

---

<div class="post-metadata">

**Author:** ![hellopascal](https://avatars.discourse-cdn.com/v4/letter/h/a8b319/32.png) [@hellopascal](https://discourse.julialang.org/u/hellopascal)\
**Post date:** [September 14, 2018, 2:33pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/5 "2018-09-14T14:33:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [September 14, 2018, 4:46pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/6 "2018-09-14T16:46:17Z")

</div>

> [@hellopascal](#):
>
> 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.

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:

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 14, 2018, 7:18pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/7 "2018-09-14T19:18:40Z")

</div>

> [@hellopascal](#):
>
> 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.

This is something on our mind: [Vector of crossing functions for condition/affect! · Issue #331 · SciML/DifferentialEquations.jl · GitHub](https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/331)

---

<div class="post-metadata">

**Author:** ![math](https://avatars.discourse-cdn.com/v4/letter/m/ea666f/32.png) [@math](https://discourse.julialang.org/u/math)\
**Post date:** [April 19, 2019, 5:46pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/8 "2019-04-19T17:46:28Z")

</div>

If I want the event to occur at two times, say, t=0.05 and t=5, what should I do? I read the links you posted but couldn’t figure out.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2019, 7:39pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/9 "2019-04-19T19:39:00Z")

</div>

> [@math](#):
>
> If I want the event to occur at two times, say, t=0.05 and t=5, what should I do? I read the links you posted but couldn’t figure out.

`(t-0.05)*(t-5.0)` is a good rootfinding condition if you don’t step too far. Or you can use flags to turn one on after the other for safety.

---

<div class="post-metadata">

**Author:** ![math](https://avatars.discourse-cdn.com/v4/letter/m/ea666f/32.png) [@math](https://discourse.julialang.org/u/math)\
**Post date:** [April 19, 2019, 7:50pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/10 "2019-04-19T19:50:48Z")

</div>

So, I tried to use two callbacks. Below is my code. Earlier it was atleast running and giving me plots but the plots didn’t make sense. It seemed like it was not doing what it was supposed to do. Now, the code stopped working and it says

> **ERROR:** UndefVarError: Pkg not defined
> 
> Stacktrace:
> 
> [1] top-level scope at **none:0**

Here is my code:

* * *

# Packages

using Plots  
using DifferentialEquations

# Parameters

i = 1.25 # 0.089;4/45/4=1.25  
f = .25  
u = .4  
w = 1.14  
a = 1.19  
d = .25  
s = .02  
b = 0.55  
e = .15  
r = 5  
m = 0  
q = .00001  
h =1  
v = .03

# Time

maxtime = 30  
tspan = (0.0, maxtime)

# Dose

stim = 5000

# Initial conditions

x0 = [10000 1 20.0 200.0 200.0]

# Model equations

function system(dy, y, p, t)

V = y[1]  
X = y[2]  
Y = y[3]  
P = y[4]  
Q = y[5]

V, X, Y, P, Q, = y  
T = X+Y  
E = (V)./(v+T+V)  
A = (V)./(s+V)  
B = (V)./(r+T+V)  
C = X./(X+Y+q)  
D = Y./(X+Y+q)

dy[1] = i - f.\*V - u.\*X.\*A- w.\*Y.\*B  
dy[2] = a.\*X.\*A - a._X._(1-A) + m.\*P - m.\*P.\*C + h.\*P.\*E  
dy[3] = b.\*Y.\*B -e.\*Y - b._Y._(1-B) - m.\*P.\*D  
dy[4] = a._X._(1-A) - m.\*P + m.\*P.\*C - d.\*P - h.\*P.\*E  
dy[5] = b._Y._(1-B) + m.\*P.\*D - e.\*Q  
end

# Events

function condition1(y, t, integrator)  
t - 6  
end  
function condition2(y, t, integrator)  
t-20  
end

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

cb1 = ContinuousCallback(condition1, affect!)  
cb2 = ContinuousCallback(condition2, affect!)  
cb = CallbackSet(cb1,cb2)

# Solve

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

# Plotting

plot(sol, layout = (2, 3))  
\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_\_–

 ![35%20PM](https://global.discourse-cdn.com/julialang/original/3X/3/9/399174f41f313501aae1c4dc3bca937c8d843eb2.jpeg)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 19, 2019, 7:59pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/11 "2019-04-19T19:59:50Z")

</div>

It looks like it’s doing what you told it to: it sets `integrator.u[4] = 500` at `t=6` and `t=20`. Did you mean for it to add 500? If so, then `integrator.u[4] += 500`.

But the easier way to do this is just a `DiscreteCallback` if you have a lot of doses.

---

<div class="post-metadata">

**Author:** ![math](https://avatars.discourse-cdn.com/v4/letter/m/ea666f/32.png) [@math](https://discourse.julialang.org/u/math)\
**Post date:** [April 19, 2019, 8:08pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/12 "2019-04-19T20:08:00Z")

</div>

Does

> integrator.u[1]

mean solution of the first ODE?

---

<div class="post-metadata">

**Author:** ![math](https://avatars.discourse-cdn.com/v4/letter/m/ea666f/32.png) [@math](https://discourse.julialang.org/u/math)\
**Post date:** [April 19, 2019, 8:10pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/13 "2019-04-19T20:10:02Z")

</div>

Also, I edited my code in the comments. Would you please look at t and tell why u1(t) is going to zero when I am actually defining it to be 5000 at t=6 and t=20? Thanks

---

<div class="post-metadata">

**Author:** ![math](https://avatars.discourse-cdn.com/v4/letter/m/ea666f/32.png) [@math](https://discourse.julialang.org/u/math)\
**Post date:** [April 19, 2019, 9:14pm UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/14 "2019-04-19T21:14:00Z")

</div>

Also, how can I plot the variables in separate windows? I tried sol (1,: ) or sol(:,1) similar to what’s been used in Matlab. But it didn’t work. All the examples that I came across use plot(sol). It won’t work for me because of the different scales of variables.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 20, 2019, 2:52am UTC](https://discourse.julialang.org/t/time-dependent-events-in-ode/14951/15 "2019-04-20T02:52:20Z")

</div>

It was fine, you just needed to reduce the tolerance.

```julia
using Plots
using DifferentialEquations

i = 1.25 # 0.089;4/45/4=1.25
f = .25
u = .4
w = 1.14
a = 1.19
d = .25
s = .02
b = 0.55
e = .15
r = 5
m = 0
q = .00001
h =1
v = .03

maxtime = 30
tspan = (0.0, maxtime)

stim = 5000

x0 = [10000 1 20.0 200.0 200.0]

function system(dy, y, p, t)

  V = y[1]
  X = y[2]
  Y = y[3]
  P = y[4]
  Q = y[5]

  V, X, Y, P, Q, = y
  T = X+Y
  E = (V)./(v+T+V)
  A = (V)./(s+V)
  B = (V)./(r+T+V)
  C = X./(X+Y+q)
  D = Y./(X+Y+q)

  dy[1] = i - f*V - u*X*A- w*Y*B
  dy[2] = a*X*A - a*X*(1-A) + m*P - m*P*C + h*P*E
  dy[3] = b*Y*B -e*Y - b*Y*(1-B) - m*P*D
  dy[4] = a*X*(1-A) - m*P + m*P*C - d*P - h*P*E
  dy[5] = b*Y*(1-B) + m*P*D - e*Q
end

function condition1(y, t, integrator)
  t - 6
end
function condition2(y, t, integrator)
  t-20
end

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

cb1 = ContinuousCallback(condition1, affect!)
cb2 = ContinuousCallback(condition2, affect!)
cb = CallbackSet(cb1,cb2)

prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb, reltol=1e-6, isoutofdomain=(u,p,t)->any(x->x<0,u))

plot(sol,layout=(2,3))

savefig("result.png")

```

![result](https://global.discourse-cdn.com/julialang/original/3X/3/0/30d61912a2447d114c759d192e370399bfee2940.png)

It was just hard to see what it was doing because of the divergence when values hit negative. You may want to do something like `isoutofdomain=(u,p,t)->any(x->x<0,u)`, so

```julia
sol = solve(prob, Rodas4(), callback = cb, isoutofdomain=(u,p,t)->any(x->x<0,u))

```

also works just fine.

> [@math](#):
>
> Also, how can I plot the variables in separate windows? I tried sol (1,: ) or sol(:,1) similar to what’s been used in Matlab. But it didn’t work. All the examples that I came across use plot(sol). It won’t work for me because of the different scales of variables.

Are you looking for `plot(sol,vars=1)`? The docs are here: [http://docs.juliadiffeq.org/latest/basics/plot.html](http://docs.juliadiffeq.org/latest/basics/plot.html)