Hi numerical ODE solvers enthusiasts, we’ve got a problem here that we think necessitates your kind input Here are two related questions:
The following ODE system represents biological species biomasses
B, and it is theoretically guaranteed that species with null biomass
B = 0 have a constant, null trajectory.
However, and depending on how we parametrize
DifferentialEquations’s solver, we get “zombie” species rising out from the dead to nonzero biomasses even if their initial state is zero. We are not sure why.
using DifferentialEquations """ ODE system defining the evolution of the biomasses of my 3-species system. Species 2 and 3 feed on species 1. """ function dBdt!(dB, B, p, t) println("t: $t, B: $B") # <- Sometimes called with B > 0 : why? # @assert B == 0 # <- This would fail, why? r, K, B0, x, y2, y3 = p # unpack parameters # Compute feeding rates F21 = (x * y2 * B * B^2) / (B0^2 + B^2) F31 = (x * y3 * B * B^2) / (B0^2 + B^2) # Compute growth rates dB = r * B * (1 - B / K) - F21 - F31 dB = F21 - x * B dB = F31 - x * B # We have (B = 0) ⇒ (dB = 0), # but even (B = 0) is not always true, why? end # Second species biomass is set to zero, expected to stay at zero. B0 = [0.1, 0.0, 0.2] p = (1, 1, 0.5, 0.314, 7.992, 8) problem = ODEProblem(dBdt!, B0, (0, 100_000), p) no_zombies = solve(problem, reltol=1e-6) # All B values are zero. zombies! = solve(problem, reltol=1e-3) # Some B values are becoming positive.
Sneaking into the above example, we realized that
dBdt! is sometimes called with non-zero values of
B. As a consequence, guaranteeing that
B = 0 ⇒ dB = 0 is not enough to prevent the apparition of zombies.
Here is my first question : Why would
solve ever call
dBdt! with non-zero values of
B in the above example?
Now, here is more general context:
dBdt! is actually generated from user input (foodweb, species capabilities, etc.). As a consequence, even though zombies did not appear in the particular above example when
reltol = 1e-6, there is no way we can calculate generic, adequate values for
In fact, the real need we have is exactly the one described in this other thread : Given a positive extinction threshold
xt > 0, trigger a callback as soon as
B[i] < xt is detected and force it to
B[i] = 0 instead… and forever.
In this other thread, @ChrisRackauckas suggested that using callbacks from
DiffEqCallbacks to this end was a fine approach. However, the example above demonstrates that hard-setting
B[i] = 0 is not a guarantee that the trajectory of
B[i] will stick to
So here is my second question : Are
DiffEqCallbacks really the right tool for the job?
If yes, then how to use them properly and avoid zombie species?
If not, should we rely on other approaches, like stopping simulation whenever
B[i] < 0 is detected then re-generating
dBdt! with fewer variables?