I’m using a VectorContinuousCallback to perform an “extinction event”, where for my large 50-equation system, if the values of one species (i’m modeling a bunch of interacting species) goes below 10e-30 then the value of that species goes to zero.
Here’s a mini example of sort of how I’m doing this with a higher threshold:
function dynamics!(du,u,p,t)
α, β, γ, δ = p
du[1] = (α - β * u[2]) * u[1]
du[2] = (-γ + δ * u[1]) * u[2]
end
function condition(out, u , t, integrator)
out[1] = 20-u[1]
out[2] = 20-u[2]
end
function affect!(integrator, idx)
for i in 1:2
if idx == i
integrator.u[i] = 0
end
end
end
cb = VectorContinuousCallback(condition,affect!,2)
u0 = [100, 100]
p = [0.55, 0.028, 0.84, 0.026]
tspan = (0.0,20.0)
prob = ODEProblem{true}(dynamics!,u0,tspan,p)
sol = solve(prob,Tsit5(), callback = cb)
plot(sol)
Which gives me this:
which is what I want.
However, in a few of my simulations in my larger system (with 50 species and where the extinction threshold is 10e–30, I’ve noticed that values that were taken to zero would sometimes pop back up (i.e. to something small like 2.5e-29 or something), but this is messing with my counting of how many extant species I have in the system.
My question is: How do I keep a value at zero once it’s been taken there by a callback? I first thought that I might want to resize the system with deleteat! but that causes problems since all the equations are linked.
For now, within my dynamics! function I put a condition, so it is like this:
function dynamics!(du,u,p,t)
for i in 1:2
if u[i] < 10e-30
du[i] = 0.0
u[i] = 0.0
else
*populate du*
end
end
and I think that worked, but there must be a more elegant way to do this in the solver options?
To be clear, 99% of the time, the value stays at zero, but I need to make sure it does so 100% of the time over millions of sims.
Any thoughts/help much appreciated!!!
Cheers