# How to bring differential equations to wanted value with DiffEqFlux?

**URL:** <https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167>\
**Category:** Numerics\
**Tags:** question\
**Created:** [January 29, 2021, 7:15am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167 "2021-01-29T07:15:55Z")\
**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:** [January 29, 2021, 7:15am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/1 "2021-01-29T07:15:55Z")

</div>

Hello,  
I have a population system (bacteria/phages) that I have described with these ODEs that include the logistic factor:

```julia
function sir!(du, u, p, t) # SIRl with logistic term
    μ, κ, φ, ω, η, β = 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
# parms
mu = 0.47 # maximum growth rate susceptible 
kappa = 2.2*10^7 # maximum population density
phi = 10.0^-9 # adsorption rate resident commensal phage
eta = 1.0 # lyse rate resident commensal phage
beta = 50.0 # burst size resident commensal phage
omega = 0.05 # outflow
tmax = 4000.0 # time span 0-tmax
tmax = 4000.0 # time span 0-tmax
r_s0 = 50000.0 # initial susceptible population u[1]
r_i0 = 0.0 # initial infected population u[2]
r_v0 = 1000.0 # initial infectious agent population u[3]
Tp = 1000 # time of infection
# execute
Vp = (mu/phi) * (1-(r_s0/kappa)) - omega/phi
tspan = (0.0, tmax)
u0 = [r_s0, r_i0, 0]
parms = [mu, kappa, phi, omega, eta, beta]
# inoculum
condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += Vp
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u[1]-1
function affect2!(integrator)
  integrator.u[1] = 0
  integrator.u[2] = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u[3]-1
function affect3!(integrator)
  integrator.u[3] = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(SIR!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp])

```

The model includes

1. a time point Tp when inoculation of Vp phages is given
2. a term that accounts for the extintion of a species by setting the equation to zero when the cell count goes below 1.  
With this model, I get a system as:  
 ![treatment_8](https://global.discourse-cdn.com/julialang/original/3X/f/7/f7e33612790aa68d2f698091ef2f7775093ed1cb.png)  
My problem is that I would like to find the Vp and Tp that bring the blue line to zero.  
I have looked into [DiffEqFlux](https://github.com/SciML/DiffEqFlux.jl/blob/master/docs/src/examples/optimization_ode.md) but I don’t know how to implement it.  
In the example given in the link, I can change the `lotka_volterra` function with my `SIR`. I got the initial conditions `u0`, the intermediary points `tspan`, and the parameters `p`. The problem is how do I set the `loss` function. How do I tell `DiffEqFlux` to bring u[1] to zero using u[3]?  
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:** [January 29, 2021, 8:05am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/2 "2021-01-29T08:05:44Z")

</div>

> [@Luigi\_Marongiu](#):
>
> How do I tell `DiffEqFlux` to bring u[1] to zero using u[3]?

What do you mean “using `u[3]`”? You can do a loss function to find what parameters cause `u[1]` to go to zero at a time `T`, but the question you’re asking doesn’t seem to make sense?

---

<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:** [January 29, 2021, 8:17am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/3 "2021-01-29T08:17:39Z")

</div>

If the parameters are `μ, κ, φ, ω, η, β = p`, then no, I am not interested in that. I would like to know how many phages = Vp = u[3] I need to add and at what time point Tp in b1. As you can see in the other [post](https://discourse.julialang.org/t/how-to-bring-a-system-of-differential-equations-to-zero-with-differentialequations/53659/15) (of which the present post is, ideally, a more focused continuation) there is an amount of phage = u[3] that collapses u[1]. As you suggested, how do I train DiffEqFlux to use u[3] as such a parameter? and how do I tell DiffEqFlux to take into account also the time? Thanks

---

<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:** [January 29, 2021, 8:21am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/4 "2021-01-29T08:21:00Z")

</div>

Make `Vp` a parameter.

```julia
parms = [mu, kappa, phi, omega, eta, beta, Vp]
affect1!(integrator) = integrator.u[3] += integrator.p[7]

soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp], sensealg=ForwardDiffSensitivity())

```

should be all you need.

---

<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:** [January 29, 2021, 8:22am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/5 "2021-01-29T08:22:42Z")

</div>

Thanks, I’ll try…

---

<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:** [January 30, 2021, 11:02am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/6 "2021-01-30T11:02:18Z")

</div>

Something weird happened. Until the day before yesterday, when I was introducing 200 million phages, the system went into a cyclic phase and then reached an equilibrium:

 ![treatment_8](https://global.discourse-cdn.com/julialang/original/3X/f/7/f7e33612790aa68d2f698091ef2f7775093ed1cb.png)  
Yesterday I ran this code:

```julia
parms = [mu, kappa, phi, omega, eta, beta, Vp]
# modification
condition1(u, t, integrator) = t==Tϕ
affect1!(integrator) = integrator.u[3] += integrator.p[7]
cb1 = DiscreteCallback(condition1, affect1!)
# extintion
condition2(u, t, integrator) = u[1]-1
function affect2!(integrator)
  integrator.u[1] = 0
  integrator.u[2] = 0
end
cb2 = ContinuousCallback(condition2, affect2!)
condition3(u, t, integrator) = u[3]-1
function affect3!(integrator)
  integrator.u[3] = 0
end
cb3 = ContinuousCallback(condition3, affect3!)
# run
modification = CallbackSet(cb1, cb2, cb3)
prob = ODEProblem(sir!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=modification, tstops=[Tp], sensealg=ForwardDiffSensitivity())

```

and, for the same amount I got the extinction:

 ![before](https://global.discourse-cdn.com/julialang/original/3X/b/8/b8b8373e3da744bfffceff42f3894df8925d1479.png)  
Since the amount inserted is the same, I thought the solver itself has changed. So, I closed Julia’s core and started a new session but the results are the same.  
What happened?  
Can the solver be reversed as it was?  
I think the question I have (found the amount of phage that extinguishes the bacterium) cannot be found mathematically, but from the [other post](https://discourse.julialang.org/t/how-to-bring-a-system-of-differential-equations-to-zero-with-differentialequations/53659/15) (obtained with the solver as it was) it can be seen that there must be an amount that triggers such extinction.  
So, I’ll rephrase the question: is there a computational savvy way to run the sir function in a loop providing different Vp?  
Thanks

---

<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:** [January 30, 2021, 11:35am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/7 "2021-01-30T11:35:50Z")

</div>

> [@Luigi\_Marongiu](#):
>
> Can the solver be reversed as it was?

You can change versions at any time. Are you sure it’s a version thing?

> [@Luigi\_Marongiu](#):
>
> So, I’ll rephrase the question: is there a computational savvy way to run the sir function in a loop providing different Vp?

[https://diffeq.sciml.ai/stable/features/ensemble/](https://diffeq.sciml.ai/stable/features/ensemble/)

---

<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:** [January 30, 2021, 12:56pm UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/8 "2021-01-30T12:56:34Z")

</div>

I don’t think is a solver’s version problem, for I did not install a new version of `DifferentialEquations`: I just added `BoundaryValueDiffEq, Flux, Optim, DiffEqFlux, DiffEqSensitivity`. But how can it be that now the solver brings the system to zero even when using a fresh session but before it did not? I assumed that the code I used changed the solver in a stable manner (if possible)…

---

<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:** [January 30, 2021, 6:46pm UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/9 "2021-01-30T18:46:44Z")

</div>

repro?

---

<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:** [January 31, 2021, 7:50am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/10 "2021-01-31T07:50:34Z")

</div>

Hello, I don’t know what “repro” means but I tried again. If I give the value Vp = 200 million phages, I now get extinction, whereas before that was not the case. Moreover, If I use up to Vp/21 I still get extinction but I don’t with amounts below Vp/22. Standing that I don’t understand why I was getting a different result before running DiffEqFlux, it is still evident that the minimal amount of phages is not Vp but another value that DiffEqFlux did not help finding. I’ll try the ensemble simulations you’ve suggested. Thanks

---

<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:** [February 10, 2021, 8:26pm UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/11 "2021-02-10T20:26:30Z")

</div>

Hello, how do i set the code for the parallel ensemble simulation? I can see that:  
`EnsembleProblem(prob::DEProblem; output_func = (sol, I),-> (sol, false), prob_func = (prob, i, repeat) -> (prob), reduction = (u, data, I) -> (append!(u, data), false)` where `prob` is the one I have set in the first post. Is this verbatim OK? and how do I tell what parameter to check and in what range? The example in the link gives the code `function prob_func(prob, I, repeat) @. [...]` but I don’t understand what it does.

---

<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:** [February 11, 2021, 10:04am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/12 "2021-02-11T10:04:15Z")

</div>

> [@Luigi\_Marongiu](#):
>
> and how do I tell what parameter to check and in what range? The example in the link gives the code `function prob_func(prob, I, repeat) @. [...]` but I don’t understand what it does.

You just say how to change the base problem to new ones. So, remake it with new parameters.

---

<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:** [February 23, 2021, 6:35am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/13 "2021-02-23T06:35:57Z")

</div>

Hello, could you please write an example of how should I set the function? I tried to focalize the problem with [this post](https://discourse.julialang.org/t/how-to-test-multiple-parameters-with-sde/55089/2) but I got no answer… 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:** [February 23, 2021, 8:45am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/14 "2021-02-23T08:45:25Z")

</div>

[https://diffeq.sciml.ai/stable/features/ensemble/#Random-Initial-Conditions](https://diffeq.sciml.ai/stable/features/ensemble/#Random-Initial-Conditions)

---

<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:** [February 23, 2021, 10:32am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/15 "2021-02-23T10:32:17Z")

</div>

I have seen that but I don’t know how to use it…

---

<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:** [February 23, 2021, 12:55pm UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/16 "2021-02-23T12:55:02Z")

</div>

What’s the exact question here? There’s plenty of examples already, so without a simple exact question it’s hard to do much more. It’s literally just a function that spits out the problem to solve for trajectory `i`.

---

<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:** [February 26, 2021, 6:53am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/17 "2021-02-26T06:53:05Z")

</div>

I have been looking at [example 2](https://diffeq.sciml.ai/stable/features/ensemble/) of the vignette since it covers a Lotka-Volterra system. This example has a main function that resembles the one I called `growth` in my [second post](https://discourse.julialang.org/t/how-to-test-multiple-parameters-with-sde/55089). The problem is that the example than sets another function:

```julia
function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
end

```

What is the use of this function in my case?  
I just want to test over u[3] in `growth`. How do I write that down?  
shall I do:

```julia
function g(du,u,p,t)
  du[1] = p[3]*u[1]
  du[2] = p[4]*u[2]
  du[3] = p[5]*u[3]
end

```

or simply

```julia
function g(du,u,p,t)
  du[3] = p[1]*u[3]
end

```

What about the indices?

---

<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:** [February 26, 2021, 10:15am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/18 "2021-02-26T10:15:58Z")

</div>

> [@Luigi\_Marongiu](#):
>
> I just want to test over u[3] in `growth` . How do I write that down?

What do you even mean by this? What are you trying to do mathematically? Do you mean you want to sample over the parameters of the noise? Then yes, make there be a parameter `p` that is used in the diffusion equation.

> [@Luigi\_Marongiu](#):
>
> What about the indices?

That’s completely dependent on the model. Is that the same parameter that is used in the drift equation `f`? Then use the same index. Is it a different one? Then use a different index. It’s literally just the same `p` that you pass in. There’s nothing fancy there.

---

<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:** [February 26, 2021, 10:41am UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/19 "2021-02-26T10:41:07Z")

</div>

I just need to pass different values to u[3] at the selected step. I got:

```julia
condition1(u, t, integrator) = t==Tp
affect1!(integrator) = integrator.u[3] += i

```

I just want to test several combinations of `Tp` and `i`.

---

<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:** [February 26, 2021, 12:00pm UTC](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167/20 "2021-02-26T12:00:27Z")

</div>

```julia
condition1(u, t, integrator) = t==integrator.p[4]
affect1!(integrator) = integrator.u[3] += integrator.p[5]

```

and follow the examples.

[Next page](https://discourse.julialang.org/t/how-to-bring-differential-equations-to-wanted-value-with-diffeqflux/54167.md?page=2)