I'm doing Numerical integration of ODEs. for a special system that has an always positive coordinate
s and a conjugated momentum
At some point the equations become so stiff that
s becomes negative and if I improve tolerances of the numerical integrator, it crashes. To solve this problem, I'd like to approximate the ODEs with
ps_new=-ps around the critical region (for example
The problem is that these are not smooth. Numerical integrators (I use BDF) break down.
It was hard for me to find anything in literature on it. If you know of a better approach, please let me know.
Note, there is conservation of energy, so there exist
f(s,ps) such that