Geodesics in 2+1 Schwarzschild metric

Question

Would appreciate if someone could explain the query that, I am trying to understand the code for timelike geodesic in 2+1 Schwarzschild $(\theta=\frac{\pi}{2})$ where we have taken $M=1$. From this how to get orbit plots? Is this possible to draw orbit without involving effective potential and angular momentum? Here I have taken $\dot{t}, \dot{r}, \dot{\phi}$ such way initial condition that satisfy the timelike condition

$(-\dot{t}^{2}+ \dot{r}^{2}+\dot{\phi}^{2}=-1)$. We are currently using NDSolve. $$ds^{2}=-(1-\frac{2M}{r})dt^{2}+(1-\frac{2M}{r})^{-1}dr^{2}+r^{2}d \phi^{2}$$

Geodesic equation: $$\ddot{t}+\frac{2 \dot{t} \dot{r}}{(-2+r)r}=0$$ $$\ddot{r}+\frac{(-2+r)\dot{t}^{2}}{r^{3}}+\frac{\dot{r}^{2}}{2r-r^{2}}+(2-r)\dot{\phi}^{2}=0$$ $$\ddot{\phi}+\frac{2 \dot{r} \dot{\phi}}{r}=0$$

Initial Condition: $$t[0]=1, \dot{t}[0]=2$$ $$r[0]=1, \dot{r}[0]=\sqrt{3}$$

$$\phi[0]=2, \dot{\phi}[0]=0$$

Reference: Introduction to Black Hole physics, Valeri P. Frolov and Andrei Zelnikov, Chapter-7, section-7.2.3 (pg.no:192). We are not exactly looking at this picture, similar to this but more concessions on trajectories.

Answer 1

Generally we do NDSolve for geodesics which are 2nd order ODEs. Here is a link of answer given by Alex Trounev Interpreting Mathematica code on black holes

Please look at it.

Answer 2

We use part of the code from KerrOrbitGRProject and data to generate first picture on this page:

M = 1; a = 0; max\[Tau] = 1750; 
eq = {Derivative[2][t][\[Tau]] == -((2*Derivative[1][r][\[Tau]]*Derivative[1][t][\[Tau]])/
       ((-2 + r[\[Tau]])*r[\[Tau]])), Derivative[2][r][\[Tau]] == 
     ((-r[\[Tau]]^2)*Derivative[1][r][\[Tau]]^2 + (-2 + r[\[Tau]])^2*(Derivative[1][t][\[Tau]]^2 - 
         r[\[Tau]]^3*(Derivative[1][\[Theta]][\[Tau]]^2 + Sin[\[Theta][\[Tau]]]^2*Derivative[1][\[CurlyPhi]][\[Tau]]^2)))/
      ((2 - r[\[Tau]])*r[\[Tau]]^3), Derivative[2][\[Theta]][\[Tau]] == 
     -((2*Derivative[1][r][\[Tau]]*Derivative[1][\[Theta]][\[Tau]])/r[\[Tau]]) + Cos[\[Theta][\[Tau]]]*Sin[\[Theta][\[Tau]]]*
       Derivative[1][\[CurlyPhi]][\[Tau]]^2, Derivative[2][\[CurlyPhi]][\[Tau]] == 
     -((2*(Derivative[1][r][\[Tau]] + Cot[\[Theta][\[Tau]]]*r[\[Tau]]*Derivative[1][\[Theta]][\[Tau]])*Derivative[1][\[CurlyPhi]][\[Tau]])/
       r[\[Tau]]), Derivative[1][t][0] == 1.384573420066829, Derivative[1][r][0] == 0, 
    Derivative[1][\[Theta]][0] == 0, Derivative[1][\[CurlyPhi]][0] == 0.088, t[0] == 0, r[0] == 6.5, 
    \[Theta][0] == Pi/2, \[CurlyPhi][0] == 0}; 
coords = {t, r, \[Theta], \[CurlyPhi]};
soln = NDSolve[eq, coords, {\[Tau], 0, max\[Tau]}];

Now we visualize trajectory together with BH

 bg = Show[
  Graphics[{Yellow, Disk[], Red, Circle[], Black, Disk[{0, 0}, .8]}, 
   Background -> Blue], 
  ParametricPlot[
   Evaluate[{r[\[Tau]] Sin[\[Theta][\[Tau]]] Cos[\[CurlyPhi][\[Tau]]] \
/. soln[[1]], 
     r[\[Tau]] Sin[\[Theta][\[Tau]]] Sin[\[CurlyPhi][\[Tau]]] /. 
      soln[[1]]}], {\[Tau], 0, max\[Tau]}, PlotStyle -> {LightBlue}]]

Now we prepare frames for animation

lst = Table[
   Evaluate[{r[\[Tau]] Sin[\[Theta][\[Tau]]] Cos[\[CurlyPhi][\[Tau]]] \
/. soln[[1]], 
     r[\[Tau]] Sin[\[Theta][\[Tau]]] Sin[\[CurlyPhi][\[Tau]]] /. 
      soln[[1]]}], {\[Tau], 0, max\[Tau], 1}];

frames = Table[
   Show[bg, 
    ListPlot[Take[lst, {i, i + 10}], 
     PlotRange -> {{-15, 15}, {-15, 15}}, Frame -> True, 
     AspectRatio -> Automatic, 
     PlotStyle -> {Yellow, PointSize[Medium]}]], {i, 1, 
    Length[lst] - 10, 20}];

ListAnimate[frames]