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

code image

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.

  • $\begingroup$ Welcome to Mathematica.SE. Are you sure you are posting on the right site? There is nothing in your question making it clear that it is concerned with Mathematica software. $\endgroup$
    – m_goldberg
    Commented Oct 10, 2019 at 13:44
  • $\begingroup$ @m_goldberg I believe this related to NDSolve. $\endgroup$ Commented Oct 11, 2019 at 10:06

2 Answers 2


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.

  • $\begingroup$ I have already seen this code but it is difficult to understand (Actually I am not aware of the command Block and Join), Can anyone help to clear the query in simpler way so that I can do better? I will very much appreciate it if anyone can write the NDSolve part with the plotting of the orbits as in the code already mentioned above. $\endgroup$ Commented Oct 11, 2019 at 10:51

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], 
   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}]]

Figure 1

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[
    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}];


Figure 2


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.