Black hole orbit in Mathematica

Question

I am currently investigating the motion of a particle around a black hole. The Lagrangian for the system is

$$ \mathcal{L}=F(r)\dot{t}^2-\frac{\dot{r}^2}{F(r)}-(r\dot{\phi})^2 $$

Where $ F(r) = 1 - 1/r $. Dots represent differentiation wrt a parameter $ s $. This gives constants of motion:

$$ F(r)\dot{t} = E $$ $$ r^2\dot{\phi} = L $$

The Lagrangian can be re-written in terms of $u=1/r$ and becomes

F[u_] := 1 - u;
L[u_, t_, p_] := F[u]*D[t, s]^2 - D[u, s]^2/(F[u]*u^4) - (D[p, s]/u)^2;

Solving the Euler-Lagrange equation can be done in one line as

Solve[D[D[L[u[s], t[s], p[s]], u'[s]], s] == D[L[u[s], t[s], p[s]], u], u''[s]]

Which gives $$ \ddot u = \frac{(5 u-4) \dot u^2}{(u-1)u} $$ Solving numerically in Mathematica is done below (u has been denoted w here)

 solu = NDSolve[{w''[s] == ((-4 + 5 w[s])(w'[s])^2)/((-1 + w[s]) w[s]), w'[0] == 0, w[0] == 1/6}, {w, w'}, {s, 0, 10}]

but this gives a constant w (so u) for the whole range. Why is this? Obviously the particle is going to move if we put it around a black hole.

Answer 1

If you change the intial condition to w'[0]=.01 you get a variing solution

W = NDSolveValue[{w''[s] == ((-4 + 5 w[s]) w'[s]^2)/((-1 + w[s]) w[s]), 
w'[0] == 1/100, (* small initial velocity*)
w[0] == 1 /6}, w , {s, 0, 10}]

Plot[W[s], {s, 0, 10}, PlotRange -> {0, 1.1}, AxesLabel -> {s, w[s]}]

Decreasing w'[0] furthermore moves the significant change of w to greater values of s. The limit case w'[0]->0 moves the step to infinity. That's what NDSolve calculated in your numerical solution w[s]==1/6~constant!

Answer 2

I wonder if the problem lies in the numerical solution of the differential equation. With all due respect, I question the differential equation itself you obtained. Because this might tell a different story:

<< VariationalMethods`
Clear[F, ℒ]

F[r_] := 1 - 1/r;
ℒ[r_, t_, ϕ_][s_] := F[r[s]] (t'[s])^2 - (r'[s])^2/F[r[s]] - (r[s] ϕ'[s])^2
EulerEquations[ℒ[r, t, ϕ][s], {r[s], t[s], ϕ[s]}, s] // FullSimplify

$$ \begin{align} r \left(\frac{2 \ddot r}{r-1}-2 \dot\phi^2\right) + \frac{\dot t^2}{r^2} &= \frac{\dot r^2}{(r-1)^2}, \\ 0 = \frac{\dot r \dot t}{r} + (r-1) \ddot t &= r\frac{\mathrm d}{\mathrm ds}\left(F(r)\dot t\right), \\ 0 = r \left(2 \dot r \dot\phi + r \ddot\phi\right) &= \frac{\mathrm d}{\mathrm ds}\left(r^2\dot\phi\right). \end{align} $$

The last two equations establish the constants of motion, indeed, as you said. By using the constants of motion, however, the first equation becomes

$$ \ddot r = \frac{2L^2(r-1)^2+r^3(\dot r^2-E^2)}{2r^4(r-1)}, $$

or in terms of $ u = 1/r $

$$ \ddot u = \frac{u^5 \left[2 L^2 (u-1)^2 u-E^2\right] + (5u-4) \dot u^2}{2(u-1) u}. $$

It is this equation that is worthy of solving, I suppose.