1

I have an ODE (for a phase space density) with a physical interpretation where the values should always be non-negative. Unfortunately, solve_ivp uses large enough timesteps such that the values become negative. Forcing small timesteps would lead to a huge increase in computation time. So I'm looking for other ways to enforce the nonnegativity:

  1. I looked at solve_bvp, but it seems like it can only enforce some linear S*y = 0 and not y > 0
  2. Is there a way to change which steps in the solve_ivp integration are accepted/discarded?
  3. Alternatively, is there a way to adapt the step size when the values are close to zero?

Thanks for any help,

Nik

niklasrb
  • 11
  • 1
  • Use `solve_ivp` with the Radau or LSODA implicit methods. Do you still get negative values? Did you select problem-relevant values for the tolerances `rtol, atol`? Alternatively, replace the density with `y=a*exp(b*u)` with sensibly chosen values for `a` and `b`, then `du/dt = dy/dt/(b*y)`. Then `y` can not become negative. – Lutz Lehmann May 11 '21 at 15:53
  • See for instance https://stackoverflow.com/questions/41648878/replacing-negative-values for earlier discussions of this problem. – Lutz Lehmann May 11 '21 at 16:19

0 Answers0