Time-dependent events in ODE

Specific Domains Modelling & Simulations
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.)

1 Like
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.

1 Like
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: Vector of crossing functions for condition/affect! · Issue #331 · SciML/DifferentialEquations.jl · GitHub

8

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.

9

(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.

10

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
35%20PM1190×842 151 KB

11

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.

12

Does

integrator.u[1]

mean solution of the first ODE?

13

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

14

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.

15

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

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

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

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

also works just fine.

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