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$ Commented Oct 26, 2018 at 18:27
  • $\begingroup$ Can you please tell me where you got this equation from? $\endgroup$
    – ayr
    Commented Jun 12, 2023 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 you have read our privacy policy.

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