How to modulate the resolution of Julia DifferentialEquations solver

Specific Domains Numerics
,
1

Hello,
I have this model with delayed differential equations:

using DifferentialEquations
# function
function baseInfect!(du, u, p, t)
    μ, κ, φ, ω, η, β = p
    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
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

t_in  = 1000     # injection time
gran  = 1000        # granularity of the ODE

# 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=[gran])

I have set t_in to determine the time of modification of the equations and gran to determine the “granularity” of the model via tstops. In fact, the model works only if t_in and gran are the same. I reckon the problem is that if the solver has a too large “granularity” (I don’t know the exact term), it might not kick on the modification. But the model does not work even if gran goes to 1.
In other words, I am having this sort of model:

ModelTherapy
ModelTherapy1920×1440 107 KB

only when t==t_in is the same as tstops=[gran], otherwise I get this:

ModelTherapy
ModelTherapy1920×1440 69.4 KB

How can I module the solver?
Thank you

2

Actually it is a similar thing in here:

  • Simply use ContinuousCallback 's and forget about tstops
3

Of course you need a tstop at a point if a DiscreteCallback is asking if a value is hit exactly. I don’t quite get the question. You can use a PresetTimeCallback if you want that part to be done automatically.

4

I ran it without tstops but in this case I get the second graph even with t=1000…
But according to this post just using tstops=t_in, works just fine. (The present post was a special case of that post…)

5

Here, we are saying that if the discrete callback condition cheks equality namely time then this preset time should be in tstops so that solver steps into this time exactly. What purpose gran has is another thing. I advise you to make it continuos callback (like t-t_in in condition) as @MatFi suggests and you don’t need to add the time, t_in, in tsops.

1 Like
6

But if you can specify the time explicitly, do that. It’s more efficient than a continuous callback and a lot less involved.

1 Like
7

Thank you!