# How to debug ODE

**URL:** <https://discourse.julialang.org/t/how-to-debug-ode/50576>\
**Category:** Numerics\
**Tags:** package, debugging, calculus\
**Created:** [November 22, 2020, 10:39am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576 "2020-11-22T10:39:43Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 22, 2020, 10:39am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/1 "2020-11-22T10:39:43Z")

</div>

Hello,  
I have a set of ODEs to model a bacterial infection:

```julia
function kinInfect!(du, u, p, t)
    μ , φ, η, β, ω, χ, ε = p
    H = h = 0
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    μ = growth rate (replication coef), a
    φ = adsorption rate (transmission coef), b
    η = lysis rate, k
    β = burst size, L
    ω = outflow (decay rate phage), m
    H = host response against bacterium
    h = host response against phage
    =#
    x = χ(t, ε)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = (x * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end
Chi(x, y) = (x < y) ? 0.0 : 1.0
# parms
tmax = 20.0 # time span 0-tmax
mu = 0.5 # growth rate (replication coef), a
phi = 1e-7 # adsorption rate (transmission coef), b
eta = 5 # lysis rate, k
beta = 100 # burst size, L
omega = 5 # outflow (decay rate phage), m
tmax = 20.0 # time span 0-tmax
t_in = 2.5 # time of inoculum
v_in = 1e8 # amount of inoculum
s0 = 1000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 0.0 # initial phage population
tspan = (0.0, tmax)
parms = [mu, phi, eta, beta, omega, Chi, t_in]
u0 = [s0, i0, v0]
# delay infection
condition(u, t, integrator) = t==t_in # time of inoculum
affect!(integrator) = integrator.u[3] += 1e8 # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# implement model
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

```

the function Chi is there to neutralize the first term of du[3] until the inoculum time is reached. The first time I ran it I got:

 ![Model_2](https://global.discourse-cdn.com/julialang/original/3X/3/2/3274624fc988cbe6cd15208f41c906c1b08c93a2.png)  
But the second time I got:  
 ![Model_2](https://global.discourse-cdn.com/julialang/original/3X/1/e/1ef71691e673909351c6fa24ea4f8be9bfb56e60.png)  
I don’t remember having changed anything between the first and second time (apart from removing the vertical line).  
How can I debug an ODE function? Are there some techniques?  
Thank you

---

<div class="post-metadata">

**Author:** ![MatFi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/matfi/32/10002_2.png) [@MatFi](https://discourse.julialang.org/u/MatFi)\
**Post date:** [November 23, 2020, 1:16pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/2 "2020-11-23T13:16:50Z")

</div>

I guess that it’s the `plot(...)` command messing it up (probably an erroneous assignment `t_in=...` )

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [November 23, 2020, 5:10pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/3 "2020-11-23T17:10:34Z")

</div>

I don’t know the reason but one thing to care about is that you have to stop at 2.5sec(or any time with this kind of discontinuity) to correctly get the changes in `Chi(x, y) = (x < y) ? 0.0 : 1.0`. Otherwise, the solver allowed to pass that time by large margin if accepted by its error metrics.

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [November 23, 2020, 9:19pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/4 "2020-11-23T21:19:59Z")

</div>

that may be truly the thing. How would I do that? Tx

---

<div class="post-metadata">

**Author:** ![thisrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thisrod/32/10642_2.png) [@thisrod](https://discourse.julialang.org/u/thisrod)\
**Post date:** [November 23, 2020, 9:59pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/5 "2020-11-23T21:59:43Z")

</div>

The `t==t_in` looks dodgy. Unless you get lucky with rounding errors, that will never be true at a discrete time step, and the inoculation won’t happen. The plotted phage number never gets to 1e8, so it looks like the inoculation isn’t happening.

If I were doing this, I would solve for 2.5s, add 1e8 to the returned solution, then solve again for the remaining time, with that as the initial condition, and with the parameter chi set to 1 instead of 0.

---

<div class="post-metadata">

**Author:** ![tomaklutfu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tomaklutfu/32/2411_2.png) [@tomaklutfu](https://discourse.julialang.org/u/tomaklutfu)\
**Post date:** [November 24, 2020, 3:54am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/6 "2020-11-24T03:54:47Z")

</div>

Simply adding more to `tstops=[2.5,1000 #=And more of the discontinuity time=#]` or by continuous callbacks.

---

<div class="post-metadata">

**Author:** ![thisrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thisrod/32/10642_2.png) [@thisrod](https://discourse.julialang.org/u/thisrod)\
**Post date:** [November 24, 2020, 6:33am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/7 "2020-11-24T06:33:08Z")

</div>

Here’s the simplest way to do what I think you want to do.

```
tmax = 20.0 # time span 0-tmax
t_in = 2.5 # time of inoculum
v_in = 1e8 # amount of inoculum

μ = 0.5 # growth rate (replication coef), a
φ = 1e-7 # adsorption rate (transmission coef), b
η = 5 # lysis rate, k
β = 100 # burst size, L
ω = 5 # outflow (decay rate phage), m
H = 0 # host response against bacterium
h = 0 # host response against phage

s0 = 1000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 0.0 # initial phage population

function kinInfect!(du, u, _, t)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = ((t > t_in) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    @assert du[3] ≈ (η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end

prob = ODEProblem(kinInfect!, [s0, i0, v0], (0, t_in))
soln1 = solve(prob, AutoVern7(Rodas5()))
s1, i1, v1 = soln1[end]
v1 += v_in
prob = ODEProblem(kinInfect!, [s1, i1, v1], (t_in, tmax))
soln2 = solve(prob, AutoVern7(Rodas5()))

tt = [soln1.t; soln2.t]
uu = zeros(length(tt), 3)
for (j,u) in pairs([soln1.u; soln2.u])
    uu[j,:] = u
end
plot(tt, uu)

```

Note:

- Parameters defined at the top level are still available in the body of `kinInfect!`

- The boolean expression `t>t_in` can be promoted to a number, 0 if false, 1 if true. This makes a lot of simulation code simpler. Edit: But it can be left out of this code. If the initial state is in the subspace with no phages and no infections, `u[2] == u[3] == 0`, then it stays there.

- You can use the solution of one problem as the initial value for another.

- The plotting code should be `plot(soln1); plot!(soln2)`, but `plot!` has a bug and replaces the existing plot instead of adding to it.

 ![yours](https://global.discourse-cdn.com/julialang/original/3X/5/3/530fa494a926e9c75495c2e2cbe1e81ee0b7aad5.png)

 ![mine](https://global.discourse-cdn.com/julialang/original/3X/b/0/b09cb53114c91e1feb7fcb810260909c4d010e43.png)

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 8, 2020, 2:43pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/8 "2020-12-08T14:43:31Z")

</div>

I tried but did not work. I re-arranged as:

```julia
function kinInfect!(du, u, p, t)
    μ, φ, η, β, ω, ε = p
    H = h = 0
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    μ = growth rate (replication coef), a
    φ = adsorption rate (transmission coef), b
    η = lysis rate, k
    β = burst size, L
    ω = outflow (decay rate phage), m
    H = host response against bacterium
    h = host response against phage
    =#
    X(x, y) = (x < y) ? 0.0 : 1.0
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = (X(t,ε) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    println("t= ", t, "> ", X(t,ε), ": ", du[3])
end 
# parameters
mu = 0.5 # growth rate, a
phi = 1e-7 # adsorption rate, b
eta = 5 # lysis rate, k
beta = 100 # burst size, L
omega = 5 # outflow, m
s0 = 1000 # initial susceptible population, x0
i0 = 0.0 # initial infected population, y0
v0 = 0.0 # initial phage population
t_in = 2.5 # time of inoculum, tφ
v_in = 100 # amount of inoculum, vφ
tmax = 20.0 # duration
# instantiate
tspan = (0.0, tmax)
parms = [mu, phi, eta, beta, omega, t_in]
u0 = [s0, i0, v0]
condition(u, t, integrator) = t==t_in # time of inoculum
affect!(integrator) = integrator.u[3] += v_in # amount of inoculum
cb = DiscreteCallback(condition,affect!)
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

```

and I got:

```julia
julia>
t= 0.02> 0.0: 0.0
t= 0.0> 0.0: 0.0
t= 0.0011528391772627539> 0.0: 0.0
t= 0.02510627541594441> 0.0: 0.0
t= 0.03765941312391662> 0.0: 0.0
t= 0.10502364904863687> 0.0: 0.0
t= 0.1405332743433084> 0.0: 0.0
t= 0.20382196654005486> 0.0: 0.0
t= 0.21327524779360946> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23500388065939562> 0.0: 0.0
t= 0.32717504217939464> 0.0: 0.0
t= 0.37547864554281657> 0.0: 0.0
t= 0.6346915537961186> 0.0: 0.0
t= 0.7713301294032545> 0.0: 0.0
t= 1.0148606280227241> 0.0: 0.0
t= 1.051236198718852> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.124538601794153> 0.0: 0.0
t= 1.2650322206698672> 0.0: 0.0
t= 1.3386598925940385> 0.0: 0.0
t= 1.7337700218279322> 0.0: 0.0
t= 1.942043929289803> 0.0: 0.0
t= 2.3132498519821483> 0.0: 0.0
t= 2.368695996757698> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.4779124745446595> 0.0: 0.0
t= 2.639783896538048> 1.0: 0.0
t= 2.7246149091335132> 1.0: 0.0
t= 3.1798458226432613> 1.0: 0.0
t= 3.419811123806132> 1.0: 0.0
t= 3.847500516672675> 1.0: 0.0
t= 3.9113834628925255> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.037152472484574> 1.0: 0.0
t= 4.2222816752185315> 1.0: 0.0
t= 4.319301257399964> 1.0: 0.0
t= 4.8399400152491> 1.0: 0.0
t= 5.114383671461072> 1.0: 0.0
t= 5.603523726846436> 1.0: 0.0
t= 5.676585412203474> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.820092888234482> 1.0: 0.0
t= 6.024921293647171> 1.0: 0.0
t= 6.132264522152217> 1.0: 0.0
t= 6.708303337160416> 1.0: 0.0
t= 7.011949956102049> 1.0: 0.0
t= 7.55313835905778> 1.0: 0.0
t= 7.633974382156477> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.792595310607519> 1.0: 0.0
t= 8.015959631677706> 1.0: 0.0
t= 8.133016869457805> 1.0: 0.0
t= 8.761184230157331> 1.0: 0.0
t= 9.092309304086916> 1.0: 0.0
t= 9.682472469951104> 1.0: 0.0
t= 9.77062373676918> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.94339866563282> 1.0: 0.0
t= 10.182812927080112> 1.0: 0.0
t= 10.30828136356051> 1.0: 0.0
t= 10.981585962999713> 1.0: 0.0
t= 11.336504173041954> 1.0: 0.0
t= 11.969073769707178> 1.0: 0.0
t= 12.063559184117928> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.248606138382897> 1.0: 0.0
t= 12.502256710845652> 1.0: 0.0
t= 12.63518588785822> 1.0: 0.0
t= 13.348527257153213> 1.0: 0.0
t= 13.724549999534794> 1.0: 0.0
t= 14.394734137049063> 1.0: 0.0
t= 14.494837945860569> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.690752829766621> 1.0: 0.0
t= 14.95667016398437> 1.0: 0.0
t= 15.096027911328429> 1.0: 0.0
t= 15.84386708814722> 1.0: 0.0
t= 16.23807461659283> 1.0: 0.0
t= 16.940669440458176> 1.0: 0.0
t= 17.045614356315234> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.250894453033002> 1.0: 0.0
t= 17.52739867156927> 1.0: 0.0
t= 17.672304625668488> 1.0: 0.0
t= 18.449917291594815> 1.0: 0.0
t= 18.859819282055923> 1.0: 0.0
t= 19.590386295643068> 1.0: 0.0
t= 19.69950935087289> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.899629072326462> 1.0: 0.0
t= 19.91010889559556> 1.0: 0.0
t= 19.915600995490593> 1.0: 0.0
t= 19.94507339686609> 1.0: 0.0
t= 19.960609146882575> 1.0: 0.0
t= 19.988298464713438> 1.0: 0.0
t= 19.992434352185413> 1.0: 0.0
t= 20.0> 1.0: 0.0
t= 20.0> 1.0: 0.0
julia> 

```

Looks like the test worked but the condition did not…

---

<div class="post-metadata">

**Author:** ![thisrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thisrod/32/10642_2.png) [@thisrod](https://discourse.julialang.org/u/thisrod)\
**Post date:** [December 9, 2020, 1:05am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/9 "2020-12-09T01:05:28Z")

</div>

However you rearrange the code, the way you’re using `condition` and `affect!` simply won’t work. In my last post, I showed you a different way that does work.

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 9, 2020, 12:24pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/10 "2020-12-09T12:24:34Z")

</div>

Sorry, but I do not understand the code. For starter, why `affect` and `condition` don’t work, since they work in other instances. Then, if I understand the code, you are splitting the analysis in two: one before t=2.5 and one after. Why not a single one with a split in the middle (the role of having `condition`)? And what is the role of `@assert`? Thank you

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 9, 2020, 12:46pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/11 "2020-12-09T12:46:16Z")

</div>

> adding more to `tstops=[2.5,1000 #=And more of the discontinuity time=#]` or by continuous callbacks.

How would I do that?

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 9, 2020, 1:05pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/12 "2020-12-09T13:05:36Z")

</div>

> [@thisrod](#):
>
> However you rearrange the code, the way you’re using `condition` and `affect!` simply won’t work.

I wrote this:

```julia
using DifferentialEquations
using PyPlot

# function
function baseInfect!(du, u, p, t)
    μ, κ, φ, ω, η, β = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end

# parameters
mu = 0.47 # maximum growth rate susceptible (B. longum)
kappa = 2.2*10^7 # maximum population density
phi = 10.0^-9 # adsorption rate
omega = 0.05 # outflow
eta = 1.0 # lyse rate
beta = 50.0 # burst size
tmax = 4000.0 # time span 0-tmax
t_in = 1000.0 # injection time
s0 = 50000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 0.0 # initial phage population
v_in = 80.0 # amount of injection

#=
SCENARIO 1
    one host one predator
=#
# instantiate solver
u0 = [s0, i0, v0]
parms = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)
# delay infection
condition(u, t, integrator) = t==t_in # time of inoculum
affect!(integrator) = integrator.u[3] += v_in # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = ODEProblem(baseInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

```

 ![Model_1](https://global.discourse-cdn.com/julialang/original/3X/4/b/4bded72a5e97c60747dc15b6049d01b76651497c.png)  
In this case, `condition` and `affect` do work. Is the problem of time scale? In this example, the injection is at t=1000, in the original problem it was 2.5. In that case, would [tomaklutfu]'s suggestion be the easiest?

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 9, 2020, 1:48pm UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/13 "2020-12-09T13:48:32Z")

</div>

Interestingly, I tried with t=900 or 1100, and I get this:

 ![Model_1](https://global.discourse-cdn.com/julialang/original/3X/e/d/ed2bb1188618d153e84d26de78bde5303ab05ae8.png)

So I got lucky the first time choosing 1000 as t\_in…

---

<div class="post-metadata">

**Author:** ![thisrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/thisrod/32/10642_2.png) [@thisrod](https://discourse.julialang.org/u/thisrod)\
**Post date:** [December 10, 2020, 12:12am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/14 "2020-12-10T00:12:35Z")

</div>

That way could work too, when you set `tstops=[t_in]`.

---

<div class="post-metadata">

**Author:** ![Luigi\_Marongiu](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/luigi_marongiu/32/7909_2.png) [@Luigi\_Marongiu](https://discourse.julialang.org/u/Luigi_Marongiu)\
**Post date:** [December 10, 2020, 6:59am UTC](https://discourse.julialang.org/t/how-to-debug-ode/50576/15 "2020-12-10T06:59:22Z")

</div>

> [@thisrod](#):
>
> `tstops=[t_in]`

Looks like it is fine this way. Thanks