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?

Improve this question
$\endgroup$
3
  • $\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, 2018 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, 2018 at 1:56
  • $\begingroup$ Do any constants of motion exist? $\endgroup$
    – bbgodfrey
    Oct 3, 2018 at 2:29

2 Answers 2

Reset to default
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.

Improve this answer
$\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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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