# Orbit followed by a particle around Schwarzschild Black Hole

The following is a equation which describes various possible orbits of a particle around the Schwarzschild black hole spacetime in general relativity. I want to solve it from -Pi to Pi but fails, any suggestions?

d = 6 m; m = 2;

s =
NDSolve[
{u'[ϕ] - Sqrt[2 m u[ϕ]^3 - u[ϕ]^2 + 1/d^2] == 0, u[0] == 0},
u, {ϕ, -0.4, 2.5}]

P1 = PolarPlot[Evaluate[{1/(u[ϕ])} /. s], {ϕ, -0.4, 2.5},AspectRatio->1/GoldenRatio]


• this produces a plot fine for me. You may want to specify an aspect ratio eg, AspectRatio -> 1 to PolarPlot Feb 6 '17 at 18:21
• Yes. Actualy, what I am getting after running the programme is not clear. Feb 6 '17 at 18:24
• I added the plot to the question. Its not clear if you aren't seeing a plot or if you think something is wrong with the plot. Feb 6 '17 at 18:27
• What if we use the range for \phi from -\pi to + \pi. Feb 6 '17 at 18:31
• Can you add a reference for the equation? Feb 7 '17 at 10:22

# Solution 1

The equation can be solved with Method -> "StiffnessSwitching", together with a higher WorkingPrecision:

d = 6 m; m = 2;
s = NDSolve[{u'[ϕ] - Sqrt[2 m u[ϕ]^3 - u[ϕ]^2 + 1/d^2] == 0, u[0] == 0},
u, {ϕ, -Pi, Pi}, Method -> "StiffnessSwitching", WorkingPrecision -> 64]

Plot[Evaluate[u[ϕ] /. s // Re], {ϕ, -Pi, Pi}]


Plot[Evaluate[1/u[ϕ] /. s // Re], {ϕ, -Pi, Pi}]


PolarPlot[Evaluate[1/u[ϕ] /. s // Re], {ϕ, -Pi, Pi}, Exclusions -> ϕ == 0,
PlotRange -> {{-30, 30}, Automatic}]


I added Re because even with WorkingPrecision -> 64 the solution still involves very small imaginary numeric error e.g.

s[[1, 1, -1]][-Pi]
(* -0.0736592722609752787242287643805033631850830491206634648394182567 -
4.668866930386631320069097140978040*10^-31 I *)


# Solution 2

Inspired by Solution 1, I found the problem can be resolved by adding a Re in the equation:

d = 6 m; m = 2;
s = NDSolve[{u'[ϕ] - Re@Sqrt[2 m u[ϕ]^3 - u[ϕ]^2 + 1/d^2] == 0,
u[0] == 0}, u, {ϕ, -Pi, Pi}]


The resulting graph is the same as above so I'd like to omit it here.

# P.S.

I still have a feeling that the equation is wrong, at least incomplete, for example, don't you think the following

d = 6 m; m = 2;
s = NDSolve[{u'[ϕ]^2 == (Re@Sqrt[2 m u[ϕ]^3 - u[ϕ]^2 + 1/d^2])^2,
u[0] == 0}, u, {ϕ, -Pi, Pi}]

PolarPlot[Evaluate[1/u[ϕ] /. s], {ϕ, -Pi, Pi}, Exclusions -> ϕ == 0,
PlotRange -> {{-30, 30}, Automatic}]


seems to be more natural?

• Now this is fine. Its working now. Thank you very much xzczd. Feb 8 '17 at 7:22