6
$\begingroup$

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$.)

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ Have you tried using an ODE library that supports event location, and asking it to stop integration at a time $r(t)=10^{-8}$ or something small but not zero (example in Mathematica, I got $t_0$ to within $10^{-7}$). ODE event location requires a sign change, which $r$ doesn't have at $t_0$, which makes this a non-standard ODE problem. $\endgroup$
    – Kirill
    Jul 1, 2018 at 9:11
  • $\begingroup$ @Kirill: That sounds like what I want. If someone could point me to an open-source library that had such a routine, that would be great. Then I could read it and see what algorithm it uses, and reimplement it myself. What I'm describing, in the terminology you're using, is the fact that I was unsuccessful in rolling my own routine for event location. (I don't own a copy of Mathematica, and it's not open source. Since I'm targeting the browser in javascript, I probably can't just use someone else's software out of the box.) $\endgroup$
    – user28109
    Jul 1, 2018 at 11:46
  • 1
    $\begingroup$ A good adaptive timestep selection method might be another good place to look, because it seems from what you've written you don't have that yet. It'd be much easier than event location, which IIRC is rather tricky to do well. $\endgroup$
    – Kirill
    Jul 1, 2018 at 12:26
  • $\begingroup$ @Kirill: Your answer is very helpful, thanks! It would also be great if you or someone else could point me to information on the theory behind event location. Yes, it does seem tricker than I had imagined when I naively tried to roll my own. In this particular problem, the diffeq solver can actually overshoot into negative values of $r$, which have no physical significance. But that is not guaranteed, and I 've also observed behavior where the solution "bounces" back from the singularity. $\endgroup$
    – user28109
    Jul 1, 2018 at 12:32

2 Answers 2

Reset to default
4
$\begingroup$

Here is an example of this can be done in Julia with an ODE library that can do event location:

using OrdinaryDiffEq

function main()
    function rhs(u, _p, _t)
        _q, r, qd, rd = u
        qdd = -rd*qd/(r*(r-1))
        rdd = -qd^2*(r-1)/2r^3 + rd^2/(2r*(r-1))
        [qd, rd, qdd, rdd]
    end
    cb = ContinuousCallback((t,u,i) -> u[2] - 1e-6, i -> terminate!(i))
    prob = ODEProblem(rhs, [0, 0.9, 9, -1.0540925533894594], (0.0, 1.0))
    sol = solve(prob, Vern9(), reltol=1e-12, abstol=1e-12)
    sol.t[end] - 0.5692099788303082
end

The absolute error in $t_0$ comes out to 7.149836278586008e-14, compared with the value you gave. Like I said in the comment, what makes this problem nontrivial is that you don't have a change of sign in the equation $r(t)=0$, but you can work around that by solving $r(t_1)=\epsilon$, which only differs from $t_0$ by a small error of magnitude $|r(t_1)/r'(t_1)|$.

Also, integrating an ODE next to a singularity should produce lots of extra steps as the solver tries to avoid the singularity itself, so this approach isn't all that efficient, if that matters to you.

Improve this answer
$\endgroup$
3
  • $\begingroup$ If $r$ behaves like $t^{2/3}$ near the singularity, then it's true that the error in $t$ is of order $r(t_1)/r'(t_1)$, but I can't be sure that $r$ always behaves this way. I found this article by Shampine and Thompson, "Event Location for Ordinary Differential Equations," radford.edu/~thompson/webddes/eventsweb.html . It basically says that solvers implement this sort of thing by finding a piecewise polynomial approximation to the solution, and then solving to find a root of the polynomial. This seems suspect near a point where the function's derivative is blowing up. $\endgroup$
    – user28109
    Jul 2, 2018 at 2:32
  • 1
    $\begingroup$ The interpolating polynomial is using Horner's rule for its numerical stability, see the code with the constants here. It allows for BigFloats if you need beyond floating point accuracy. The derivative blow up makes the time steps shorten which then increases the accuracy of the interpolation (and the interpolation is 9th order like the method). $\endgroup$
    – Chris Rackauckas
    Jul 2, 2018 at 5:50
  • 2
    $\begingroup$ DiffEq is already ran as a browser-backend in a webapp, see this app. How this was done is described in this blogpost. If you want help, just come to the DiffEq chatroom. Also, a lot of this should already compile directly to Javascript (WebAssembly) via Charlotte.jl. That would be a cool project if you're up for it. $\endgroup$
    – Chris Rackauckas
    Jul 2, 2018 at 5:54
1
$\begingroup$

Similar to the Julia example, here is some example code using the ode15i function in Octave 4.4. The Octave implementation of ode15i is based on the Sundials IDA solver. IDA uses a family of backward difference formulas of orders one through five. ode15i also includes a basic event location capability.

The Octave code is shown below:

function black_hole_example
y0=[0 9 .9 -1.0540925533894594]';
yp0=[0 0 0 0]';
opts=odeset('Events', @eventFunc, 'RelTol', 1e-12, 'AbsTol', 1e-12);
[y0,yp0] = decic(@odeFunc,0,y0,[1 1 1 1]',yp0,[0 0 0 0]');
tspan=[0 1];
[t,y,te] = ode15i(@odeFunc,tspan,y0,yp0,opts);
err=te - 0.5692099788303082;
fprintf('Final time=%16.12f, final r=%10.2e, time error=%10.2e\n', 
te, y(end,3), err);
figure; plot(t, y(:,3));
end

function f=odeFunc(t,y,yp)
  q=y(1); qp=y(2); r=y(3); rp=y(4);
  qpp=-1/(r*(r-1))*rp*qp;
  rpp=-.5*r^-3*(r-1)*qp^2 + .5/(r*(r-1))*rp^2;
  f=yp-[qp qpp rp rpp]';
end

function [res,isterminal,direction]=eventFunc(t,y,yp)
res=y(3)-1e-9;
isterminal=true;
direction=0;
end

When run in Octave the results are:

Final time=  0.569209978864, final r= 9.99e-010, time error= 3.32e-011

enter image description here

Just as with the Julia ODE solver, a relatively large number of time steps are required to achieve this result but execution time is less than a second on a fairly modest Windows PC.

Although I ran this code in Octave, it will also run in MATLAB with few, if any, changes.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.