1

"Mathematica version 13.3.1 for Microsoft Windows (64-bit) (July 24, 2023)"

I have the following Mathematica code:

$Assumptions = {r \[Element] Reals, r >= 0, rh \[Element] Reals, 
   rh > 0, g \[Element] Reals, g > 0, l \[Element] Reals, l > 0};
M[rh_, g_] := (1 + rh^2/l^2) (rh^2 + g^2)^(3/2)/(2 rh^2)
f[r_, rh_, g_] := 1 - (2 M[rh, g] r^2)/(g^2 + r^2)^(3/2) + r^2/l^2
T[rh_, g_] := D[f[r, rh, g], r]/(4 \[Pi]) /. r -> rh
S[rh_, g_] := \[Integral]D[M[rh, g], rh]/T[rh, g] \[DifferentialD]rh
G[rh_, g_] = M[rh, g] - (T[rh, g] S[rh, g]) // Simplify;
{rhc, gc} = 
 Block[{l = 1}, 
    SolveValues[{D[T[rh, g], {rh, 2}] == 0, 
      D[T[rh, g], rh] == 0}, {rh, g}, Reals]][[1]] // N
l = 1;
pplt[g_] := 
 ParametricPlot[Evaluate@{T[rh, g], G[rh, g]}, {rh, 0.001, 10}, 
  PlotRange -> {{0.2547, 0.2555}, {0.0342, 0.0349}}, AspectRatio -> 1]
ip[g_] := Graphics`Mesh`FindIntersections[pplt[g]// DiscretizeGraphics, Graphics`Mesh`AllPoints -> False]
Tp = Table[ip[g], {g, 0.098, 0.0985, 0.0001}]

which gives me intersection points for $g$ near $g_c$. Surprisingly, I get two intersection points for some values of $g$, which is not the case if I try to do it manually:

plot = Manipulate[Show[pplt[g]], {g, 0.0980, 0.0985, 0.0001}]

Why is there such a difference between automated results and manual results?
Can the automation be done differently than the one I used?

codebpr
  • 2,233
  • 1
  • 7
  • 26

1 Answers1

3
  • Since the original code does not work on 13.3.1 in my Windows,here we try to use
Region`Mesh`FindSegmentIntersections

Clear[intersections]; 
intersections[g_] := 
 Region`Mesh`FindSegmentIntersections[
   Cases[Normal@pplt[g], _Line, -1], "ReturnSegmentIndex" -> True, 
   "Ignore" -> {"EndPointsTouching"}][[1, 1]]
intersections /@ Range[0.098, 0.0985, 0.0001]

{{0.255261, 0.0344312}, {0.25518, 0.0344083}, {0.255099, 0.034386}, {0.25502, 0.0343616}, {0.25494, 0.0343376}, {0.25486, 0.0343135}}

cvgmt
  • 72,231
  • 4
  • 75
  • 133
  • Thank you for the insightful answer. I am still wondering about the bug in my original code, I have tried to edit using your technique, as described in this answer: https://mathematica.stackexchange.com/a/289189/78049 I hope the code works now in your system! – codebpr Jan 09 '24 at 04:18