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.