Runge-Kutta Approach to Relativistic Kepler Problem

Question

I am attempting to solve the relativistic Kepler problem using fourth order Runge-Kutta, however I am having trouble and am unsure why I simply cannot get the correct orbits no matter how I adjust the initial values. My orbits are not even incorrect but they are nonsensical, I was wondering if there are issues with the Runge-Kutta algorithm when solving such problems? The geodesic equations in Schwarzschild spacetime (in geometric units) are:

\begin{align} &\frac{d^2r}{d\tau^2} = \frac{M}{r(r - 2M)} \left(\frac{dr}{d\tau}\right)^2-\frac{M(r - 2M)}{r^3} \left(\frac{dt}{d\tau}\right)^2 + (r - 2M) \left(\frac{d\phi}{d\tau}\right)^2 \\[0.5em] &\frac{d^2t}{d\tau^2} = - \frac{2M}{r(r - 2M)} \frac{dr}{d\tau} \frac{dt}{d\tau} \\[0.5em] &\frac{d^2\phi}{d\tau^2} = - \frac{2}{r} \frac{dr}{d\tau} \frac{d\phi}{d\tau} \end{align}

We can convert the above to a system of 6 1st order ODE's by setting: \begin{align} y_1 = r \\[0.5em] y_2 = \frac{dr}{d\tau} \\[0.5em] y_3 = \phi \\[0.5em] y_4 = \frac{d\phi}{d\tau} \\[0.5em] y_5 = t \\[0.5em] y_6 = \frac{dt}{d\tau} \end{align}

So then, we obtain 6 1st order ODEs which can be solved using the fourth order Runge - Kutta. For example, I tested my code on mercury with the following parameters (geometrized units): \begin{align} &\text{Mass of sun: } 1474 m \\[0.5em] &\text{Mecury radius at perihelion: } 46.002\times 10^9 m \end{align}

I am not sure about the initial values for the Runge Kutta so I tried a range such as $$ [y_1(0), y_2(0), y_3(0), y_4(0), y_5(0), y_6(0)] = [46.002\times 10^9, 0, 0, 1\times 10^{-16}, 0, 1] $$ Despite the range of initial values i have tried, I cannot seem to get even a decent orbit.

Is there a problem with using Runge-Kutta to solve these equations or should there be a more deterministic way to obtain the initial values?

Answer 1

I have tested this model with RK method of order 8 to understand how it depends on input error and step size as well. First we should note, that relativistic model has hidden suggestion in the form $c=1$ where c is the speed of light. In this case the unit of time and coordinate is the same, for instance, meter. Second, we take Mercury semimajor axis, perihelion, eccentricity, period in days and maximal orbital speed as

Mmercury = 3.3010 10^23; RAmercury = 5.7909227 10^10; RPmercury = 4.6002 10^10; Emercury = 0.20563593; Pmercury = 87.97; vp=58983;

We use Mercury perihelion as a scale of length, therefore in this case we have parameter M of the model and initial condition as follows

c0 = 299792458.; GM = 1.32712 10^(20); mS = GM/c0^2; M=mS/RPmercury ;

$$t(0)= 0, r(0) = 1, \phi(0) = 0, t'(0) = u_1, r'(0) = 0, \phi'(0)= vp/c0$$ Here u1 is computed as a root of equation $$u.g.u=-1, u=(u1,0,0,vp/c0)$$ g is the metric tensor. We have calculated, that $u1=1.0000000514518859$. Also we define time scale $t_0=RPmercury/c0=153.446 s$, hence in equations we should normalize $t$ on $t_0$. Finally we compute orbit as shown in Figure 1

According this calculation the period is about $49569*t0=7.60617*10^6$, while astronomical data says that period is about $Pmercury *24*3600=7.60061*10^6$. Therefore we have discrepancies in time period, but it is not due to relativistic dilation , but due to error in definition of input parameter and computation method.