1
$\begingroup$

"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?

$\endgroup$

1 Answer 1

3
$\begingroup$
  • 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}}

$\endgroup$
1
  • $\begingroup$ 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: mathematica.stackexchange.com/a/289189/78049 I hope the code works now in your system! $\endgroup$
    – codebpr
    Jan 9 at 4:18

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.