Shell falling into a black hole in Kruskal coordinate

Question

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.

Answer 1

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?