# How to force variable to zero upon meeting condition in Julia DifferentialEquations

**URL:** <https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655>\
**Category:** Numerics\
**Tags:** question\
**Created:** [December 11, 2020, 9:05am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655 "2020-12-11T09:05:11Z")\
**Posts on this page:** 20\
**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:** [December 11, 2020, 9:05am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/1 "2020-12-11T09:05:12Z")

</div>

Hello,  
I have this model of bacteria and phage interactions:

```julia
using DifferentialEquations
function payneJansenODE!(du, u, p, t)
    μ , φ, η, β, ω = p
    H = h = 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] = (η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end
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ϕ = 2.5 # time of inoculum, tφ
vϕ = 1e9 # amount of inoculum, vφ
tmax = 20.0 # duration
u0 = [s0, i0, v0]
parms = [mu, phi, eta, beta, omega]
tspan = (0.0, tmax)
condition(u, t, integrator) = t==tϕ # time of inoculum
affect!(integrator) = integrator.u[3] += vϕ # amount of inoculum
cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[tϕ])

```

The model works but from the plot (the horizontal line is y=1):

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/a/e/aed3ac79616ec1d66b365af5ea407b5d9d0c5914.png) it is apparent that the bacterium went extinct since there is no such thing as 1/10 of a bacterium.  
It is possible to force u[1] to zero (which will bring also u[2] = 0) when a condition is reached? (in this case, when u[1] \<1).  
Thank you

---

<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:** [December 11, 2020, 9:26am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/2 "2020-12-11T09:26:41Z")

</div>

> [@Luigi\_Marongiu](#):
>
> It is possible to force u[1] to zero (which will bring also u[2] = 0) when a condition is reached? (in this case, when u[1] \<1).

Yes, that would be a ContinuousCallback where `condition(u,t,integrator) = u[1]-1` and `affect!(integrator) = integrator.u[1] = 0`.

---

<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 11, 2020, 9:32am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/3 "2020-12-11T09:32:23Z")

</div>

Thank you, that is I need to create `cb1 = ContinuousCallback(condition(u,t,integrator) = u[1]-1)`, `modification = CallbackSet(cb, cb1)` and then `soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])`?

---

<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:** [December 11, 2020, 9:43am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/4 "2020-12-11T09:43:40Z")

</div>

Yes.

---

<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 11, 2020, 9:54am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/5 "2020-12-11T09:54:25Z")

</div>

> [@ChrisRackauckas](#):
>
> condition(u,t,integrator) = u[1]-1

gave an error but I ran:

```julia
julia> # inoculum
       condition(u, t, integrator) = t==tϕ # time of inoculum
condition (generic function with 1 method)
julia> affect!(integrator) = integrator.u[3] += vϕ # amount of inoculum
affect! (generic function with 1 method)
julia> cb1 = DiscreteCallback(condition, affect!)
DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}(condition, affect!, DiffEqBase.INITIALIZE_DEFAULT, Bool[1, 1])
julia> # extintion
       affect!(integrator) = integrator.u[1] = 0
affect! (generic function with 1 method)
julia> cb2 = ContinuousCallback(condition(u,t,integrator) = u[1]-1)
ERROR: syntax: invalid keyword argument name "condition(u, t, integrator)"
Stacktrace:
 [1] top-level scope at none:0
julia> cb2 = ContinuousCallback(condition, affect!)
ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}(condition, affect!, affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0)
julia> # combine
       modification = CallbackSet(cb1, cb2)
CallbackSet{Tuple{ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}},Tuple{DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}}}((ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing,Int64}(condition, affect!, affect!, DiffEqBase.INITIALIZE_DEFAULT, nothing, true, 10, Bool[1, 1], 1, 2.220446049250313e-15, 0),), (DiscreteCallback{typeof(condition),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}(condition, affect!, DiffEqBase.INITIALIZE_DEFAULT, Bool[1, 1]),))
julia> # run
       prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 20.0)
u0: [1000.0, 0.0, 0.0]
julia> soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])
retcode: Success

```

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/6/4/64d451d7c9a46ac53afb47c442a65b1f96aaf991.png)  
However, there are no u[2] and u[3]:

```julia
julia> getindex.(soln.u, 2)
12-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> getindex.(soln.u, 3)
12-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

```

The condition u[1]=0 must have kicked in too soon…

---

<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:** [December 11, 2020, 10:06am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/6 "2020-12-11T10:06:38Z")

</div>

This has nothing to do with DiffEq… don’t define all of your functions with the same name 🤦‍♂️. `condition1(...) = ...`, `condition2(...) = ...`

---

<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 11, 2020, 10:43am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/7 "2020-12-11T10:43:18Z")

</div>

> [@ChrisRackauckas](#):
>
> `condition(u,t,integrator) = u[1]-1`

ops… I tried with this, but the condition u[1]=0 did not work:

```julia
# inoculum
condition1(u, t, integrator) = t==tϕ # time of inoculum
affect1!(integrator) = integrator.u[3] += vϕ # amount of inoculum
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u,t,integrator) = u[1]-1 
affect2!(integrator) = integrator.u[1] = 0
cb2 = ContinuousCallback(condition2, affect2!)
# combine
modification = CallbackSet(cb1, cb2) 
# run
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])

```

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/9/6/96f9e2c61cdbd59e0c52ae25b8e59c39b820b9ce.png)

---

<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:** [December 11, 2020, 10:49am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/8 "2020-12-11T10:49:54Z")

</div>

Does `u[1]` ever equal 1? You didn’t clearly plot `u[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:** [December 11, 2020, 11:05am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/9 "2020-12-11T11:05:42Z")

</div>

I tried with `condition2(u,t,integrator) = u[1]<1 `  
and got:

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/5/f/5fe400ab3b656b29d88fa4b25ab1d8f03814dc61.png)  
where `susceptible = u[1]`, `infected=u[2]` and `phage=u[3].

---

<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:** [December 11, 2020, 11:12am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/10 "2020-12-11T11:12:25Z")

</div>

> [@Luigi\_Marongiu](#):
>
> condition2(u,t,integrator) = u[1]\<1

Yes, `u[1]<1` isn’t `u[1]-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:** [December 11, 2020, 11:57am UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/11 "2020-12-11T11:57:15Z")

</div>

But even if `u[1]<1` is correct, the condition does not happen. The blue and red lines should go down and then disappear (since the axis are log). The model instead continues till the end and in fact the bacterium resurrects…

---

<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:** [December 11, 2020, 12:02pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/12 "2020-12-11T12:02:53Z")

</div>

> [@Luigi\_Marongiu](#):
>
> But even if `u[1]<1` is correct, the condition does not happen.

`u[1]<1` is not correct and would not cause the condition to happen.

---

<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 11, 2020, 1:52pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/13 "2020-12-11T13:52:08Z")

</div>

So what would be the condition that is triggered when u[1] gets decimal?

---

<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:** [December 11, 2020, 2:49pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/14 "2020-12-11T14:49:14Z")

</div>

You have to pick a value to trigger it at. 0.9999999?

---

<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 11, 2020, 4:16pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/15 "2020-12-11T16:16:31Z")

</div>

Sorry, I don’t understand. So it has to be `u[1]<0.9999999`? What would be the difference with `u[1]<1` or `[u]= 0.0001`?

---

<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:** [December 11, 2020, 5:06pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/16 "2020-12-11T17:06:40Z")

</div>

That’s why I said `u[1] - 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:** [December 11, 2020, 5:26pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/17 "2020-12-11T17:26:00Z")

</div>

OK, so why it does not work?

```julia
condition1(u, t, integrator) = t==tϕ # time of inoculum
affect1!(integrator) = integrator.u[3] += vϕ # amount of inoculum
cb1 = DiscreteCallback(condition1, affect1!)
condition2(u,t,integrator) = u[1]-1 # extintion
affect2!(integrator) = integrator.u[1] == 0
cb2 = ContinuousCallback(condition2, affect2!)
modification = CallbackSet(cb1, cb2)
prob = ODEProblem(payneJansenODE!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[tϕ])

```

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/a/7/a7e269b4653171ff61cf926b0595af189f7a0cd6.png)

---

<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:** [December 11, 2020, 5:28pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/18 "2020-12-11T17:28:27Z")

</div>

> [@Luigi\_Marongiu](#):
>
> `affect2!(integrator) = integrator.u[1] == 0`

Typo. You made a boolean.

```julia
affect2!(integrator) = integrator.u[1] = 0

```

---

<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 11, 2020, 5:48pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/19 "2020-12-11T17:48:00Z")

</div>

🥴 gosh! Sorry about that. Now looks much better, in fact. However, u[2] and u[3] don’t go to zero so I also forced u[2] to zero but the phage does not disappear. But, if the syntax is OK, then is a problem with the model:

 ![ModelTherapy](https://global.discourse-cdn.com/julialang/original/3X/5/5/55104c08778f5744c9ed1f18015dbb24117f14bd.png)

---

<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:** [December 11, 2020, 5:54pm UTC](https://discourse.julialang.org/t/how-to-force-variable-to-zero-upon-meeting-condition-in-julia-differentialequations/51655/20 "2020-12-11T17:54:52Z")

</div>

Yeah no worries. I was confused why you were looking for other things and did see that you had tried it with an error haha. My inner voice was “what the heck is this guy asking for?”

But yeah, getting “exact” extinctions in an ODE-based model is always hard. ODEs generally don’t go extinct: it’s a downside to that modeling choice. If you want that kind of behavior, continuous-time Markov chains are where people go. Catalyst.jl is made for that.

[https://catalyst.sciml.ai/dev/](https://catalyst.sciml.ai/dev/)

Though then the results are stochastic and require more time to interpret.