sol = solve(
sdeprob, LambaEM(), callback = PositiveDomain(save=false),
save_everystep=false, saveat=5.0, seed=42
)
help?> PositiveDomain
search: PositiveDomain CircleDomain TimeDomain position
PositiveDomain(u = nothing; save = true, abstol = nothing, scalefactor = nothing)
Especially in biology and other natural sciences, a desired property of dynamical systems is the positive invariance
of the positive cone, i.e. non-negativity of variables at time t_0 ensures their non-negativity at times t \geq t_0
for which the solution is defined. However, even if a system satisfies this property mathematically it can be
difficult for ODE solvers to ensure it numerically, as these MATLAB examples
(https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html) show.
To deal with this problem, one can specify isoutofdomain=(u,p,t) -> any(x -> x < 0, u) as additional solver option
(https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/), which will reject any step that leads to
non-negative values and reduce the next time step. However, since this approach only rejects steps and hence
calculations might be repeated multiple times until a step is accepted, it can be computationally expensive.
Another approach is taken by a PositiveDomain callback in DiffEqCallbacks.jl, which is inspired by Shampine's et al.
paper about non-negative ODE solutions (https://www.sciencedirect.com/science/article/pii/S0096300304009683). It
reduces the next step by a certain scale factor until the extrapolated value at the next time point is non-negative
with a certain tolerance. Extrapolations are cheap to compute but might be inaccurate, so if a time step is changed it
is additionally reduced by a safety factor of 0.9. Since extrapolated values are only non-negative up to a certain
tolerance and in addition actual calculations might lead to negative values, also any negative values at the current
time point are set to 0. Hence, by this callback non-negative values at any time point are ensured in a
computationally cheap way, but the quality of the solution depends on how accurately extrapolations approximate next
time steps.
Please note, that the system should be defined also outside the positive domain, since even with these approaches,
negative variables might occur during the calculations. Moreover, one should follow Shampine's et al. advice and set
the derivative x'_i of a negative component x_i to \max \{0, f_i(x, t)\}, where t denotes the current time point with
state vector x and f_i is the i-th component of function f in an ODE system x' = f(x, t).
Arguments
=========
• u: A prototype of the state vector of the integrator. A copy of it is saved and extrapolated values are
written to it. If it is not specified, every application of the callback allocates a new copy of the state
vector.
Keyword Arguments
=================
• save: Whether to do the standard saving (applied after the callback).
• abstol: Tolerance up to, which negative extrapolated values are accepted. Element-wise tolerances are
allowed. If it is not specified, every application of the callback uses the current absolute tolerances of
the integrator.
• scalefactor: Factor by which an unaccepted time step is reduced. If it is not specified, time steps are
halved.
References
==========
Shampine, Lawrence F., Skip Thompson, Jacek Kierzenka and G. D. Byrne. Non-negative solutions of ODEs. Applied
Mathematics and Computation 170 (2005): 556-569.
So you just need to set save=false there. I am not sure that’s the right default for the callback but that is how it is right now at least.