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.

  • 2
    $\begingroup$ One possible solution of your ode is w[s]==1/6! $\endgroup$ Oct 26, 2018 at 18:27
  • $\begingroup$ Can you please tell me where you got this equation from? $\endgroup$
    – dtn
    Jun 12 at 8:17

2 Answers 2


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

enter image description here

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!


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.


Your Answer

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

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