I'm working on a fun open-source project to calculate the trajectories of objects near black holes. This is obviously not the first time anyone has done this sort of thing, but I have some design goals that have led me to roll my own solution. (These considerations are spelled out in the project's README file, but briefly, I want to be able to construct browser-based educational apps, to model motion that crosses the horizon or hits the singularity, and to simulate things like rocket ships, for which the motion is not a geodesic.) Because of these considerations, I've written my own Runge-Kutta routine.
Basically everything is working fine and passes various tests, except for the following. I want to be able to simulate motion that terminates at a singularity. When I do this, I want to be able to precisely determine the time (or technically the affine parameter) at which it terminates. So in the language of ODEs, I have a vector ODE of the form
$$y'=f(t,y),$$
in which the function $f$ misbehaves at a known surface, and I want to determine accurately the time $t_0$ at which I hit the surface. "Accurately" would hopefully mean something like the precision with which a 4th-order Runge-Kutta can determine $y$ for a given $t$, i.e., with a reasonable step size I want to be able to get $t_0$ to within rounding errors. There is a paper by Lewis and Kwan that discusses something similar. They say they used software called odepack, but nothing about what routine they used or how they managed the singularity. They quote errors on the order of $10^{-4}$, which is good enough for their purposes but would be sad for a general-purpose library like the one I'm writing.
Can anyone give me pointers on how to proceed here, point me to references, or suggest keywords to search on? I'm not a whiz at numerical methods, although I took an undergraduate course a long time ago. I'm thinking that my ODE is not stiff. I ran across information about things called singular solutions, but the usage of this term doesn't seem to match up very well with my situation. What I'm dealing with is what a relativist would call an incomplete geodesic.
The following is a description of what I've tried.
I tried reducing the step size according to a seemingly reasonable rule as I approached the singularity. This sort of worked, but basically I could never get it to behave better than the equivalent of first order, i.e., the relative error in $t$ seemed to scale like the step size $h$, so that I could get down to errors of about $10^{-6}$ with reasonable computational effort, but nowhere near machine precision.
I also tried doing a change of variables. For solutions that approach the singularity, I had estimates that one of the components of the $y$ vector (it's actually the Schwarzschild radius $r$) should behave as $r\sim (t_0-t)^{2/3}$. So I tried a change of variable to $u=r^{3/2}$. This sort of helped, because with these new coordinates, $u'$ is nearly constant (should be exactly constant for trajectories that free fall from rest at infinity). However, I'm still not having much luck getting an algorithm to determine $t_0$ precisely. What tends to happen is that my techniques for predicting when I'm going to hit the singularity are a little off, so I overshoot, and I start getting garbage results, such as the particle flying back out of the black hole.
Is there some algorithm implemented in open-source software that I could try out, and possibly just reimplement? If it's of interest to see the equations concretely, the simplest example I can make up, omitting some degrees of freedom, is
$$q''=-r^{-1}(r-1)^{-1} r'q'$$
$$r''=-\frac{1}{2}r^{-3}(r-1)(q')^2+\frac{1}{2}r^{-1}(r-1)^{-1}(r')^2,$$
where $q$ is actually the Schwarzschild time coordinate, and the primes indicate differentiation with respect to proper time $t$. An example of initial conditions of interest would be $q=0$, $r=0.9$, $q'=9$, $r'=-1.0540925533894594$, for which there is an exact solution giving $t_0=0.5692099788303082$. (Physically, these conditions represent someone who has fallen into the black hole from rest at infinity, starting the simulation just inside the event horizon at $r=1$.)