Hello, I don’t know if this is more a math question than computational one, but since in the end the implementation with Julia is what I seek, I am posting it here.
I am using DifferentialEquations to model some data that looks like this:
It can be seen that the data follows an oscillatory path but the oscillations are reducing themself until the values are stable.
My question is: is there a way (easy) to determine when the oscillation ceases (or more properly: when the oscillation goes beyond a threshold)?
For example, the figure shows that at t=3000h the values are virtually stable. Is it possible to extract that t datum?
Thank you
You could subtract the final value from all of the data then look backwards from the end until you find the first value that is above some threshold.
Yes, but this involves a bit of data analysis. Is there instead some julia function (perhaps within DifferentialEquations) that does that directly?
Are you looking to automatically halt the integration when it hits a steady state?
https://docs.juliadiffeq.org/latest/types/steady_state_types/
https://docs.juliadiffeq.org/latest/solvers/steady_state_solve/
https://docs.juliadiffeq.org/latest/features/callback_library/#TerminateSteadyState-1
4 Likes
It’s a bit unclear how you want to do this. Do you want to get the information from the model of the system, from the solution of the model, or from observed data?
2 Likes
Essentially I would like to know at what time (h) the data has became stable enough. I would say I am seeking the solution of the model…
Did you try the steady state stuff from above?
2 Likes
The approach suggested by Sanders seems like it would require a single line of code, is that really too much “data analysis”?
3 Likes
not yet but i will…
Hello,
this is exactly what I am looking for but I am a bit lost in how to define the problem.
If I have the case:
function dynDelay!(du, u, h, p, t)
μ, κ, ϕ, ω, β, ρ, τ = p
#=
du[1] = susceptible
du[2] = infected
du[3] = phages
=#
history = h(p, t - τ)
delay_lysis = ϕ * history[1] * history[3] * exp(-1 * ω * τ)
du[1] = (μ * u[1]) * (1 - (u[1]+u[2])/κ) - (ϕ * u[1] * u[3]) - (ω * u[1])
du[2] = (ϕ * u[1] * u[3]) - delay_lysis - (ω * u[2])
du[3] = (β * delay_lysis) - (ρ * ϕ * u[1] * u[3]) - (ω * u[3])
end
# set constants
mu = 0.16 # maximum growth rate susceptible strain
nu = 0.12 # maximum growth rate resistant strain
kappa = 2.2e7 # maximum population density
phi = 1.0e-9 # adsorption rate
omega = 0.05 # outflow
beta = 50.0 # burst size
rho = 1.0 # reinfection rate
tau = 3.62 # latency time
tmax = 4000.0 # time span 0-tmax
s0 = 50000.0 # initial susceptible population
i0 = 0.0 # initial infected population
v0 = 80.0 # initial phage population
r0 = 500.0 # initial resistant population
#### ONE BACTERIUM ONE PHAGE ####
# set parameters
h(t, p) = [s0, i0, v0]
tspan = (0.0, tmax)
u0 = [s0, i0, v0]
parms = [mu, kappa, phi, omega, beta, rho, tau]
# instantiate model
prob = DDEProblem(dynDelay!, u0, h, tspan, p=parms; constant_lags = [tau])
algt = MethodOfSteps(Tsit5())
soln = solve(prob, algt)
how would I implement the auto-termination at steady state (which I reckon should give a value between 2500 and 3000 h)?
Thank you.
soln = solve(prob, algt, callback = TerminateSteadyState())