# NDSolve producing message Power::infy

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 and v' == 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?

• 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? Oct 3 '18 at 0:09
• @bbgodfrey I expect r[x] start to decrease at x=-1.5 monotonically to some positive value at r and this is its minimum, then it goes back. v[x] first goes up monotonically to maximum v and then goes down. They are all even functions. Oct 3 '18 at 1:56
• Do any constants of motion exist? Oct 3 '18 at 2:29

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, v[1.5] == 10, v' == 0}, {r, v}, {x, -1.5, 1.5},
Method -> {"Shooting", "StartingInitialConditions" ->
{r == 1.10478, v == 8.5}}]
Row[Plot[#, {x, -1.5, 1.5}, ImageSize -> Medium, PlotRange -> All] &
/@ Through[s[x]]] Because r varies so rapidly near the endpoints, I reran the calculation with WorkingPrecision -> 30. Results are more precise but otherwise unchanged.

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] 