I have an ODE system of 60 equations and run multiple simulations with different sets of parameters. I want to prevent the outputs from taking any negative values since it is not physically plausible. I currently use something like that, based on this topic: Setting really small values in ODE to zero

condition0(u,t,integrator) = any( (u .< 0.0))
function affect0!(integrator)
integrator.u[(integrator.u .< 0.0)] .= 0;
end
cb0 = DiscreteCallback(condition0,affect0!)

However, that way the computational time increases significantly.

Without using the callback function

0.831594 seconds (845.83 k allocations: 51.515 MiB, 8.17% gc time, 92.77% compilation time: 25% of which was recompilation)

Using the callback function

16.886240 seconds (40.82 M allocations: 2.280 GiB, 16.01% gc time, 11.04% compilation time: 11% of which was recompilation)

What would be a better way to tackle this issue? I have around 160 parameters so it is difficult to know which parts of the space cause these irregularities.

Iâ€™m just guessing here, but could you add a barrier function that prevents the solution from going negative?

e.g. instead of du/dt = f(u,t), use

\frac{du}{dt} = f(u,t) + \mathrm{barrier}.(u)

where barrier.(u) is an elementwise application of something like barrier(x) = ÎĽ/x^2 (for some strength ÎĽ) that pushes the solutions away from zero?

Or alternatively multiply f(u,t) * barrier.(u) for some function like barrier(x) = exp(-ÎĽ/x) that goes smoothly to zero as u \to 0^+ but \to 1 for large u.

That is, if negative u are â€śnot physically plausibleâ€ť, add something to your physics that enforces this.

One important point is that it is isnâ€™t enough to simply set u = 0 when it goes non-positive â€” you should also make sure your right-hand side is 0 at this point, so that the equations donâ€™t try to keep pushing u to the negative regime. i.e. your right-hand side f(u,t) should be replaced with

F = f(u,t)
return @. ifelse(u <= 0, max(0, F), F)

Of course, this introduces a discontinuity into your right-hand-side, but with continuous callbacks this is manageable I guess. I havenâ€™t read Shampineâ€™s paper closely, but it seems to be about refining the handling of this discontinuity.

Somewhat equivalent to adding a barrier function is reparametrizing variable to make negative u impossible. For example u = e^w and finding the appropriate new form for the differential equation in w.

That turns du/dt = f(u,t) into dw/dt = f(e^w, t) / e^w for a single variable (and the corresponding elementwise thing for multiple variables).

It looks like you run into a problem, however, if the system really â€śwantsâ€ť to go negative: that is, if f(u,t) < 0 as u \to 0^+. In this case, your w equations will have a divergence where it wants to go to w \to -\infty faster and faster (the right-hand side will diverge exponentially to -\infty as w gets more negative).

I suspect that this diverging solution will cause problems with your numerical ODE solver.

Doing both a barrier function and reparametrization might be a good solution. The barrier keeping the solution from -\infty and the reparametrization preventing discrete/edge effects from pulling u negative.

However, I noted that my system is indeed keep pushing to negative, and while the u stays close to 0, there a small negative perturbations that cause other u to exceed other physical limits (e.g., exceed mass conservation). I have not tried yet the combinations mentioned earlier by @Dan and @stevengj.

Oh, brilliant, that implements exactly the Shampine method I linked! It also mentions the \mathrm{ifelse}(u \le 0, \max(0, f(u,t)), f(u,t)) modification.