2
$\begingroup$

I am trying to solve geodesic equations in some 3D black hole spacetime. It is a coupled ODE system with boundary conditions. Due to the symmetry of the spacetime, I expect the solutions to be even functions, with r'[0] == 0 and v'[0] == 0.

Here's my code:

f[r_, v_] := r^2 - 1/2 Tanh[v] - 1/2

NDSolve[
  {0 == D[(r[x]^2 + 2 r'[x] v'[x] - f[r[x], v[x]] v'[x]^2)/r[x]^4, x], 
   r[x]^2 - r[x]^2 v'[x]^2 - r[x] v''[x] + 2 v'[x] r'[x] == 0, 
   r[-1.5] == 100, r[1.5] == 100, v[-1.5] == 10, v[1.5] == 10}, 
  {r, v}, {x, -1.5, 1.5}]

But I don't get results, but only the following messages:

Power::infy: Infinite expression 1/0. encountered.
Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity ComplexInfinity encountered.
Power::infy: Infinite expression 1/0. encountered.

Can someone point out what is wrong with my code?

$\endgroup$
  • $\begingroup$ Probably, an trial solution passes through r == 0 as NDSolve attempts to match the boundary conditions. Qualitatively, what do you expect the solutions to look like? $\endgroup$ – bbgodfrey Oct 3 '18 at 0:09
  • $\begingroup$ @bbgodfrey I expect r[x] start to decrease at x=-1.5 monotonically to some positive value at r[0] and this is its minimum, then it goes back. v[x] first goes up monotonically to maximum v[0] and then goes down. They are all even functions. $\endgroup$ – Zhencheng Oct 3 '18 at 1:56
  • $\begingroup$ Do any constants of motion exist? $\endgroup$ – bbgodfrey Oct 3 '18 at 2:29
2
$\begingroup$

The solution is very sensitive to initial conditions. Here, we take advantage of the symmetry described in the question to integrate from x == 0.

s = NDSolveValue[{0 == D[(r[x]^2 + 2 r'[x] v'[x] - f[r[x], v[x]] v'[x]^2)/r[x]^4, x], 
    r[x]^2 - r[x]^2 v'[x]^2 - r[x] v''[x] + 2 v'[x] r'[x] == 0, 
    r[1.5] == 100, r'[0] == 0, v[1.5] == 10, v'[0] == 0}, {r, v}, {x, -1.5, 1.5}, 
    Method -> {"Shooting", "StartingInitialConditions" -> 
        {r[0] == 1.10478, v[0] == 8.5}}]
Row[Plot[#, {x, -1.5, 1.5}, ImageSize -> Medium, PlotRange -> All] &
    /@ Through[s[x]]]

enter image description here

Because r varies so rapidly near the endpoints, I reran the calculation with WorkingPrecision -> 30. Results are more precise but otherwise unchanged.

$\endgroup$
0
$\begingroup$

There are two types of solutions in this model (maybe more). One indicated the author and found bbgodfrey. I will show the second solution, apparently connected with the fall to the center

f[r_, v_] := r^2 - 1/2 Tanh[v] - 1/2

plot[a_, b_, x1_] := 
 Block[{A = a, B = b, x0 = 10^-6, xm = x1}, 
  sol = NDSolve[{A^2 == (r[x]^2 + 2 r'[x] v'[x] - 
         f[r[x], v[x]] v'[x]^2)/r[x]^4, 
     r[x]^2 - r[x]^2 v'[x]^2 - r[x] v''[x] + 2 v'[x] r'[x] == 0, 
     r[x0] == 1/A, v'[x0] == 2*x0/A, v[x0] == B}, {r, v}, {x, 0, xm}, 
    WorkingPrecision -> 20, MaxSteps -> 10^6]; {Plot[
    Evaluate[r[x] /. First[sol]], {x, x0, xm}, PlotRange -> All, 
    AxesLabel -> {"x", "r"}], 
   Plot[Evaluate[v[x] /. First[sol]], {x, x0, xm}, PlotRange -> All, 
    AxesLabel -> {"x", "v"}]}]

Solution of the first type

plot[1001/1000, 1, 1990/1000]

Second type solution

plot[11/10, 1, 1990/100]

fig1

$\endgroup$

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.