4
$\begingroup$

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]

enter image description here

$\endgroup$
  • $\begingroup$ this produces a plot fine for me. You may want to specify an aspect ratio eg, AspectRatio -> 1 to PolarPlot $\endgroup$ – george2079 Feb 6 '17 at 18:21
  • $\begingroup$ Yes. Actualy, what I am getting after running the programme is not clear. $\endgroup$ – Emlie Feb 6 '17 at 18:24
  • 1
    $\begingroup$ 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. $\endgroup$ – george2079 Feb 6 '17 at 18:27
  • $\begingroup$ What if we use the range for \phi from -\pi to + \pi. $\endgroup$ – Emlie Feb 6 '17 at 18:31
  • 1
    $\begingroup$ Can you add a reference for the equation? $\endgroup$ – xzczd Feb 7 '17 at 10:22
6
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

seems to be more natural?

$\endgroup$
  • $\begingroup$ Now this is fine. Its working now. Thank you very much xzczd. $\endgroup$ – Emlie Feb 8 '17 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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