How to get numerical solutions for Geodesic path using numerical Metric in General Relativity? I have question regarding numerical method to calculate geodesic path, when the analytic form of the metric is unknown, please.
Suppose that we have a black-box at hand that can show us a measurement value for the metric $ds$ line element every time we move a short distance in space.
For space coordinates $(x,y,z)$, the geodesic path from $[x_1, y_1, z_1]$ to $[z_2, y_2, z_2]$ can be 'guessed' by for instance starting first with a linear interpolation along all the 3 coordinate bases, then iterating the guess over and over while trying to minimize the sum over all ds.
But when one of the variables is $t$, I am unsure if this changes anything, or places any new restrictions on the algorithms normally usable for $(x,y,z)$?
For example, suppose that the black-box is secretly (unbeknownst to us) outputting:
$$ds^2 = -\left(1-\frac{2Gm}{r}\right)dt^2 + \frac{dr^2}{\left(1-\frac{2Gm}{r}\right)} + r^2 d \Omega^2$$
Are there any special restrictions imposed by physics on minimizing sum of ds over a path from $[t_1, r_1, \theta_1, \phi_1]$ to $[t_2, r_2, \theta_2, \phi_2]$? In particular, can $t_1$ and $t_2$ be any values, or only a restricted set of values are allowed in this case?
Thank you.
 A: When you have a second-order differential equation, you need initial conditions (coordinates and velocities) and then you can use numerical methods to evolve this forward or backward in the evolution parameter.
Conceptually the easiest way would be to numerically integrate:
$$ x(\tau + \Delta \tau) \approx x(\tau) + v(\tau) \Delta \tau, $$
$$ v(\tau + \Delta \tau) \approx v(\tau) + F(x(\tau), v(\tau)) \Delta \tau $$
where $F$ is the function that defines your differential equation:
$$ \ddot{x} = F(x, \dot{x}). $$
In practice, accumulating precision loss comes from the fact that in the method above we've taken $F$ at point $\tau$ and not somewhere in between of $\tau$ and $\tau + \Delta \tau$. A better method that avoids much of this loss is Runge-Kutta.
The geodesic equation is precisely a second-order differential equation, so all of the above applies:
$$ \ddot{x}^{\mu} = - \Gamma_{\alpha \beta}^{\mu} \dot{x}^{\alpha} \dot{x}^{\beta}. $$
You can take this exact form and reconstruct the parametric solution using numerical methods, but there's a better way. The evolution parameter in the form of the geodesic equation above is the proper time $\tau$. Hovewer, for most practical implications, we're interested in how the trajectory of the test particle depends on the coordinate time $t = x^0$.
You could first solve the equation numerically for $x^{\mu}(\tau)$ and then construct the function $x^k(x^0)$ using interpolation, however, in situations where proper and coordinate times are greatly out of sync (such as e.g. in the vicinity of the event horizon), you will lose precision.
A better way is to change the form of the geodesic equation first (on a whiteboard or on a piece of paper), to change its form such that the evolution parameter is $t = x^0$ and not $\tau$. This is not super hard to do, but it is a slightly nontrivial calculation.
Here I'll draft the strategy for you, but I expect you to do the calculation yourself (feel free to ask a homework-and-exercise question if you are stuck).
The calculation roughly goes as follows. First, note that for any quantity $Q$,
$$ \frac{d Q}{d \tau} = \frac{d Q}{dt} \frac{dt}{d \tau}. $$
The quantity
$$ \gamma = \frac{dt}{d \tau} $$
will obviously play a major role in the calculation, along with its first derivative
$$ \frac{d \gamma}{d \tau} = \gamma \frac{d \gamma}{dt}, $$
as is evident from the zero'th component geodesic equation (with $\ddot{x^0}$).
We now rewrite the spatial components of the geodesic equation such that only the coordinate time derivatives, $\gamma$ and $\dot{\gamma}$ are allowed.
Finally, we use the time component of the geodesic equation to fix a relationship between $\gamma$ and $\dot{\gamma}$ and plug it back in the spatial components. This will eliminate both $\gamma$ and $\dot{\gamma}$ from the equation, leaving us with the desired form of the geodesic equation on $x^k(t)$.
Then use Runge-Kutta or your other favorite numerical method to solve that and reconstruct $x^k(t)$. You'll need to know the components of the affine connection, but those are trivially computable from the metric:
$$ \Gamma_{\alpha \beta}^{\mu} = \frac{g^{\mu \nu}}{2} \left( \frac{\partial g_{\nu \beta}}{\partial x^{\alpha}} + \frac{\partial g_{\alpha \nu}}{\partial x^{\beta}} - \frac{\partial g_{\alpha \beta}}{\partial x^{\nu}}\right), $$
and you can take the partial derivatives in the formula above numerically by approximating them with finite differences.
