# Shell falling into a black hole in Kruskal coordinate

I am modeling the motion of a shell falling into a black hole in Kruskal coordinate with NDSolve. I am sure that the physics behind it is correct, but can't get the desired result i.e. a smooth function defined on {0, rWALL}.

The code:

V[r_] := 1 - 100/r + r^2 - (0.16 (100/r + 1.4375 r^2)^2)/r^2
mfo[r_] := 1 - 100/r + r^2
h = 4.56978
rWALL = 6.03647
k = 13.9282
w[r_] := 1/(2 (1 + 3 h^2) Sqrt[4 + 3 h^2]) ((4 + 6 h^2) ArcTan[(h + 2 r)/Sqrt[4 + 3 h^2]] + h Sqrt[4 + 3 h^2] (2 Log[Abs[-h + r]] - Log[Abs[1 + h^2 + h r + r^2]]))
Twall[r_] := Sqrt[Xwall[r]^2 - Exp[k*w[r]]]
SolX = NDSolve[{Sqrt[(4 fo[r])/(k^2 Exp[k*w[r]])] Sqrt[-V[r]]Sqrt[(Twall'[r]^2 - Xwall'[r]^2)] == 1, Xwall[rWALL] == 1.02227}, Xwall, {r, 0, rWALL} ]


I only got the following warnings:

NDSolve::ndsdtc: The time constraint of 1. seconds was exceeded trying to solve for derivatives, so the system will be treated as a system of differential-algebraic equations. You can use Method->{"EquationSimplification"->"Solve"} to have the system solved as ordinary differential equations.

NDSolve::nlnum: "The function value {-1+(0. +0.00223649\ I)\ Sqrt[0. \-5.31229\Plus[<<2>>]^2]} is not a list of numbers with dimensions {1} \at \!$${r, Xwall[r], \*SuperscriptBox[\"Xwall\", \"\[Prime]\",MultilineFunction->None][r]}$$ = {6.03647,1.02227,0.}."

NDSolve::icfail: unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.


Can someone help me with it? Thanks.

• Can you share more information? What are you modeling? Also, the first error from NDSolve seems strange. Have you set some time constraints that you didn't show in your code, or perhaps globally? – MarcoB Aug 16 '17 at 3:23
• Hi I am modeling the motion of a shell in black hole spacetime in Kruskal coordinate, trying to solve the coordinate X with r (the radius of the shell) as a parameter of it. The code above is all I used to solve this problem, and I didn't set any other constraints. – Zhencheng Aug 16 '17 at 3:31
• I obtain a different error message, NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}. I recommend that you seek ways to eliminate the Sqrt and Abs, which probably are the source of your difficulties. – bbgodfrey Aug 16 '17 at 4:32
• Thanks for your advice. But I only got an interpolating function at one point 'r=rWALL' after removing the Sqrt and Abs... It is not a satisfying result. – Zhencheng Aug 16 '17 at 4:59
• Is Xwall always real? – xzczd Aug 16 '17 at 5:17

At least 3 issues here.

1. As mentioned in the comments, the Abs is causing trouble. If you look carefully, you'll find Abs' term in the equation, which can't be handled properly by NDSolve. This can be easily circumvented by a Abs[a_] :> Sqrt[a^2] replacement.

2. The coefficient of equation is too complicated, which stops NDSolve from transforming the equation to an… Er… explicit ODE i.e. an equation with the form $y'(x)=f(x,y(x))$ in a reasonable amount of time, so the warning ndsdtc is generated. For more information, check this post.

ndsdtc isn't always a problem, because NDSolve will then try to solve the equation with a DAE solver, but the DAE solver of NDSolve is weaker than the ODE solver (at least now) and unfortunately it fails to solve your problem, so icfail is generated.

3. The small imaginary part generated by numeric error is troublesome. Here is a similar problem.

The following is the fixed code:

Clear[V, fo, w, Twall]

Twall[r_] = Sqrt[Xwall[r]^2 - Exp[k w[r]]];
eq = Sqrt[(4 fo[r])/(k^2 Exp[k w[r]])] Sqrt[-V[r]] Sqrt[Twall'[r]^2 - Xwall'[r]^2] == 1;

neweq = Xwall'[r] == Re@Solve[eq, Xwall'[r]][[1, 1, -1]]

V[r_] = 1 - 100/r + r^2 - (0.16 (100/r + 1.4375 r^2)^2)/r^2;
fo[r_] = 1 - 100/r + r^2;
h = 4.56978;
rWALL = 6.03647;
k = 13.9282;
w[r_] = 1/(2 (1 + 3 h^2) Sqrt[
4 + 3 h^2]) ((4 + 6 h^2) ArcTan[(h + 2 r)/Sqrt[4 + 3 h^2]] +
h Sqrt[4 + 3 h^2] (2 Log[Abs[-h + r]] - Log[Abs[1 + h^2 + h r + r^2]])) /.
Abs[a_] :> Sqrt[a^2]

SolX = NDSolveValue[{neweq, Xwall[rWALL] == 1.02227}, Xwall, {r, 0, rWALL}]

ListLinePlot@SolX


The calculation still stops at 4.56978, but this seems to be the nature of model?

# Update

OK, the singularity at $r = h$ seems to be removable. Assuming $r = h$ is a removable singularity, $X_{wall}(r)$ should be smooth around it, so $X_{wall}(h - \epsilon) \approx 2 X_{wall}(h) - X_{wall}(h + \epsilon)$ where $\epsilon$ is a small enough positive number:

solLeft = Block[{eps = 10^-3},
NDSolveValue[Rationalize[{neweq, Xwall[h - eps] == 2 SolX[h] - SolX[h + eps]}, 0],
Xwall, {r, 0, h}, WorkingPrecision -> 32]]

{{rL, rR}} = solLeft["Domain"]

Plot[solLeft[r], {r, rL, rR}, PlotRange -> All]~Show~Plot[SolX[r], {r, h, rWALL}]
`

This time the singularity at $r=0$ causes trouble so I fail to reach $r = 0$, but I think $r=0.08$ is already acceptable?

• Thanks! It is good to see a result. I think the calculation stops at that point because there exist a singularity in the function w[r] at that point. The value of w[r] at that point should be negative infinity. But one of the purpose of this calculation is to remove the influence of that singularity and get a solution defined from r=0 to r=rWALL... Now it doesn't seem very good. Do you know any method that can remove the influence of the singularity? – Zhencheng Aug 16 '17 at 6:17
• @Zhencheng Er… what does "remove the influence of the singularity" mean? – xzczd Aug 16 '17 at 6:25
• Well, I mean I expect to get a solution as a smooth function defined on [0,rWALL], but not only on this small region. It seems that the singularity stops the calculation to move on. – Zhencheng Aug 16 '17 at 6:33
• Thanks a lot for your help. But in physics, what I am trying to do is to express the equation of motion of a shell falling into a black hole in Kruksal coordinate. In this coordinate system, there is no singularity in coordinates. Everything should be smooth. So I am really confused here. – Zhencheng Aug 16 '17 at 7:01
• Then why not look for a reformulation (change of variables or coordinates) that won't exhibit the singularity? – J. M.'s technical difficulties Aug 16 '17 at 7:27