I'm doing Numerical integration of ODEs. for a special system that has an always positive coordinate s
and a conjugated momentum ps
.
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 s_new=abs(s)
and ps_new=-ps
around the critical region (for example if(s<1.0e-8)
). The problem is that these are not smooth. Numerical integrators (I use BDF) break down.
One solution that I came across is to use arctan
for abs
and erf
for sign
.
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 df/dt=0