# How to remove an apparent singularity in NDSolve?

I'm trying to numerically solve two coupled differential equations using NDSolve. The following code works, except that I had to remove a starting point singularity at $$r = 0$$ by adding an arbitrary 0.01 value in a function (see the comment, in the code):

Clear["Global*"]

rmax = 1;

rho[p_, gamma_] := p^(1/gamma) + p/(gamma - 1)
TOV[p_, m_, r_, gamma_] := - (m + 4 Pi p r^3) (rho[p, gamma] + p)/((r + 0.01) (r + 0.01 - 2 m)) (* Notice the 0.01 to remove a singularity *)

solution[q_, gamma_] := NDSolve[{
p'[r] == TOV[p[r], m[r], r, gamma], (* Equ. 1 *)
m'[r] == 4 Pi rho[p[r], gamma] r^2, (* Equ. 2 *)
p[0] == q, (* Central pressure *)
m[0] == 0, (* Mass at center *)
WhenEvent[(p[r] < 0 || 2 m[r] >= r), "StopIntegration"]
}, {p, m}, {r, 0, rmax},
Method -> "StiffnessSwitching",
MaxSteps -> Automatic
]

pressure[q_, gamma_] := Plot[Evaluate[p[r] /. solution[q, gamma]], {r, 0, rmax}, PlotRange -> All]

Manipulate[Show[pressure[q, gamma], PlotRange -> {{0, rmax}, {0, q}}],
{{q, 0.5, Style["Central pressure", Italic, 10]}, 0, 1, 0.01},
{{gamma, 4/3, Style["gamma", Italic, 10]}, 1.01, 5/3 - 0.01, 0.01}
]

Preview of what this code is doing:

Yet, there should be a way to integrate both equations from $$r = 0$$ and up (without adding a pesky 0.01 value), since there's a finite limit at that starting point (i.e since $$m(0) = 0$$). So how should I tell Mathematica to evaluate that limit, without adding the 0.01 arbitrary value?

EDIT: Micheal's trick of the PieceWise function works for the code above, but NOT for the simpler version below (classical non-relativistic version):

Clear["Global*"]

rmax = 1;

rho[p_, gamma_] := p^(1/gamma)
classical[p_, m_, r_, gamma_] := - m rho[p, gamma]/r^2

solution[q_, gamma_] := NDSolve[{
p'[r] == PieceWise[{{0, r <= 0}, {classical[p[r], m[r], r, gamma], r > 0}}], (* Equ. 1 *)
m'[r] == 4 Pi rho[p[r], gamma] r^2, (* Equ. 2 *)
p[0] == q, (* Central pressure *)
m[0] == 0, (* Mass at center *)
WhenEvent[(p[r] < 0 || 2 m[r] >= r), "StopIntegration"]
}, {p, m}, {r, 0, rmax},
Method -> "StiffnessSwitching",
MaxSteps -> Automatic
]

pressure[q_, gamma_] := Plot[Evaluate[p[r] /. solution[q, gamma]], {r, 0, rmax}, PlotRange -> All]

Manipulate[Show[pressure[q, gamma], PlotRange -> {{0, rmax}, {0, q}}],
{{q, 0.5, Style["Central pressure", Italic, 10]}, 0, 1, 0.01},
{{gamma, 4/3, Style["gamma", Italic, 10]}, 1.01, 5/3 - 0.01, 0.01}
]

Changing r^2 to (r + 0.001)^2 solves the issue but it's an hack. Now, I don't understand what's going on here!

• Sorry, rho[p_, gamma_] := p^(1/gamma) + p/(gamma - 1); Limit[-(m + 4 Pi p r^3) (rho[p, gamma] + p)/((r) (r - 2 m)), r -> 0] produces Indeterminate. Aug 1, 2023 at 18:16
• @user64494, at $r = 0$, we have $m(0) = 0$, $p(0) = q < \infty$ so $\rho(0) = \text{finite value}$. Then (rho + p)is finite, and the factor (m + 4 Pi p r^3) goes to 0 faster than the denumerator r (r - 2 m) .
– Cham
Aug 1, 2023 at 18:20
• I usually do it with Piecewise: p'[r] == Piecewise[{{TOV[p[r], m[r], r, gamma], r != 0}}, 0], once I've found the value at the singularity. Aug 1, 2023 at 18:24
• @MichaelE2, thanks. This is what I was about to do, but I'm wondering if there's a "better" way with Mathematica, when using NDSolve (maybe with a condition, in the WhenEvent?)
– Cham
Aug 1, 2023 at 18:26
• It looks like p decays quickly. When the pressure goes negative, it is raised to the 1/gamma before the WhenEvent. Maybe change rho[p_, gamma_] := Abs[p]^(1/gamma) + p/(gamma - 1), since p is supposed to be positive. Aug 1, 2023 at 18:53

Here's summary in code of the better suggestions I made or improvements/alternatives to them:

rmax = 1;

rho // ClearAll;
rho[p_?NumericQ, gamma_] := (* clip when p<0 *)
If[p < 0, 0., p^(1/gamma) + p/(gamma - 1)];
TOV // ClearAll;
TOV[p_, m_, r_,
gamma_] := -(m + 4 Pi p r^3) (rho[p, gamma] +
p)/((r(*+0.01*)) (r - 2 m));

solution // ClearAll;
solution[q_, gamma_] := (solution["Data"] = {}; (* for saving data
since getting data out was raised in comments *)
With[{
ndsol = NDSolve[{
p'[r] ==
Piecewise[{{TOV[p[r], m[r], r, gamma], r != 0}}, 0],(*Equ.1*)
m'[r] == 4 Pi rho[p[r], gamma] r^2,(*Equ.2*)
p[0] == q,(*Central pressure*)
m[0] == 0,(*Mass at center*)
WhenEvent[p[r] < 0.001,  (* from a comment *)
solution["Data"] = {"stop" -> "p < 0.001", "m@Rmax" -> m[r],
"Rmax" -> r};
"StopIntegration"],
WhenEvent[ (* event shouldn't happen b/c of prev. event *)
(Developer`MachineComplexQ[p[r]] || p[r] < 0 || 2 m[r] >= r),
solution["Data"] = {"stop" -> "p negative-ish"};
"StopIntegration",
"LocationMethod" -> "StepEnd"],
WhenEvent[r > 1 && p[r] == 0, (* stop infinite integration case *)
solution["Data"] = {"stop" -> "p = 0",
"m@Rmax" -> Indeterminate, "Rmax" -> Indeterminate};
"StopIntegration"]
}, {p, m}, {r, 0, rmax; Infinity},
Method -> "StiffnessSwitching", MaxSteps -> Automatic]
},
If[("stop" /. solution["Data"]) === "p negative-ish",
solution["Data"] = Join[
solution["Data"],
{"m@Rmax" -> Last[m@"ValuesOnGrid" /. First@ndsol],
"Rmax" -> First@Last[m@"Grid" /. First@ndsol]}]
];
solution["Data"] =
Association@AppendTo[solution["Data"], "solution" -> ndsol]];
solution["Data"]["solution"]
);
pressure[q_, gamma_] := With[{ndsol = solution[q, gamma]},
Plot[
Evaluate[p[r] /. ndsol],
Evaluate@Flatten@{r, p@"Domain" /. ndsol},
PlotRange -> All,
PlotLabel -> Row[{
solution["Data"]["stop"], ", ",
Subscript[r, "max"], "=", solution["Data"]["Rmax"], ", ",
m@Subscript[r, "max"], "=", solution["Data"]["m@Rmax"]
}],
ImageSize -> 400]
];

Manipulate[
pressure[q, gamma],
{{q, 0.5, Style["Central pressure", Italic, 10]}, 0, 1, 0.01},
{{gamma, 4/3, Style["gamma", Italic, 10]}, 1.01, 5/3 - 0.01, 0.01},
TrackedSymbols :> Manipulate]

solution["Data"]

• Wow! There are lots of things in there that I don't understand. I'll have to study that solution...
– Cham
Aug 1, 2023 at 21:30
• Apparently, the Max values of rand mcan be found with Rmax[q_, gamma_] := (p/.solution[q, gamma])[[1]][[1]][[1]][[2]]and Mmax[q_, gamma_] := m[Rmax[q, gamma]]/.solution[q, gamma][[1]], is that right?
– Cham
Aug 2, 2023 at 13:31
• Yes, that's right. So will (p["Domain"] /. First@solution[q, gamma])[[1, 2]] and Last[m["ValuesOnGrid"] /. First@solution[q, gamma]] respectively. Aug 2, 2023 at 13:54
• @Cham real^noninteger yields complex when the base is negative. There are other functions such that f[real] can be complex, like Log[negative] and inverse trig. But the only place I saw it in your original code was in rho[] with p^(1/gamma). Once the complex numbers creep in, they propagate. Here 0^3 with round-off error yields complex: (-$MachineEpsilon)^(3 + 3$MachineEpsilon). For p^(1/gamma), I suggested Abs[p]^(1/gamma) originally, also tried Re[p]^(1/gamma), and settled for If[..] above. A somewhat higher AccuracyGoal sometimes can help keep p from going negative. Aug 2, 2023 at 16:00
• @Cham My point was that your suggestion to find rMax from m'[r] == 0 leads to this condition. Aug 2, 2023 at 21:52