In modelling biological communities with ODEs, given that individuals are somewhat finite entities, I am trying to avoid infinitesimal densities which are mathematically correct but make no sense. More specifically, I want to set an extinction threshold at a certain small value (eg 1e-12), to ensure that extinct species are actually extinct and do not come back from the dead later. Currently doing this with a callback (see below). It works but I was wondering if there is a more elegant or efficient way of setting to zero values of the solution under a certain threshold?
My current approach (simplified, this is of course the drug concentration decay example with a large threshold)
function f(du,u,p,t) du = -u du = -u end threshold = 0.2 condition(u,t,integrator) = any( (u .< threshold) .& (u.>0)) function affect!(integrator) integrator.u[(integrator.u .< threshold) .& (integrator.u .> 0)] .= 0 end cb = DiscreteCallback(condition,affect!) u0 = [1, 0.5] prob = ODEProblem(f,u0,(0.0,10.0)) sol = solve(prob,Tsit5(),callback=cb,saveat = 0.1) plot(sol)