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

1 Answer 1

7
$\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$
1
  • $\begingroup$ Now this is fine. Its working now. Thank you very much xzczd. $\endgroup$
    – Emlie
    Commented Feb 8, 2017 at 7:22

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.