2
$\begingroup$

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:

enter image description here

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!

$\endgroup$
15
  • $\begingroup$ 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. $\endgroup$
    – user64494
    Aug 1, 2023 at 18:16
  • $\begingroup$ @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) . $\endgroup$
    – Cham
    Aug 1, 2023 at 18:20
  • 2
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Aug 1, 2023 at 18:24
  • $\begingroup$ @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?) $\endgroup$
    – Cham
    Aug 1, 2023 at 18:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Aug 1, 2023 at 18:53

1 Answer 1

2
$\begingroup$

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]

enter image description here

solution["Data"]

Mathematica graphics

$\endgroup$
13
  • $\begingroup$ Wow! There are lots of things in there that I don't understand. I'll have to study that solution... $\endgroup$
    – Cham
    Aug 1, 2023 at 21:30
  • $\begingroup$ 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? $\endgroup$
    – Cham
    Aug 2, 2023 at 13:31
  • 1
    $\begingroup$ Yes, that's right. So will (p["Domain"] /. First@solution[q, gamma])[[1, 2]] and Last[m["ValuesOnGrid"] /. First@solution[q, gamma]] respectively. $\endgroup$
    – Michael E2
    Aug 2, 2023 at 13:54
  • 1
    $\begingroup$ @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. $\endgroup$
    – Michael E2
    Aug 2, 2023 at 16:00
  • 1
    $\begingroup$ @Cham My point was that your suggestion to find rMax from m'[r] == 0 leads to this condition. $\endgroup$
    – Michael E2
    Aug 2, 2023 at 21:52

Your Answer

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

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