Scipy Odeint Non-negative Solution

Apparently, getting a non-negative solution from an ODE solver is non-trivial. In Matlab, there is the NonNegative option for certain solvers to get a non-negative solution. Is the

Solution 1:

The problem is not merely that you have to avoid square rooting the negative x. The problem is that the "best" way of imposing the constraint still depends on what your system's application is and what behavior you assume to be "reasonable". If your system does not have an equilibrium at 0 then your problem may be ill-posed. What could be the meaning for it to move at non-zero speed into the negative-x domain? If the interpretation is that the solution should stay at zero, then you actually no longer have an ODE system as your intended model: you have a hybrid dynamical system with a non-smooth component, i.e. when the trajectory x(t) hits 0 at t = t_1, it must stay at x(t) for all t > t_1. This can be easily achieved with a proper dynamical systems package such as PyDSTool.

Alternatively, x=0 is a stable equilibrium and you simply need to prevent evaluation of f for x<0. This can also be hacked with event detection.

In either case, event detection at x=0 is tricky when your f is undefined for x<0. There are few standard ODE solvers that can literally be forced to avoid evaluating in a sub-domain under all circumstances, and most event detection will involve evaluations on either side of a boundary. A pragmatic solution is to choose a small number for x below which it is safe (in the context of your application) to declare x = 0. Then make the event detect when x reaches that (which, given that you can control the step size to stay small enough) should prevent x ever being evaluated at a negative value. Then you'd make a condition to make x = 0 after that point, if that is the behavior you want. Again, that's a bit of fuss in scipy/python but you can do it. It's also fairly easy to set up the desired behavior in PyDSTool, which I'd be willing to advise you on if you post in its help forums.