Let's define some simple functions
Clear["Global`*"];
x0 = 0; y0 = 0;
x1 = 1/Sqrt[3]; y1 = 0;
x2 = -(x1/2); y2 = 1/2;
x3 = x2; y3 = -y2;
r0 = Sqrt[(x - x0)^2 + (y - y0)^2];
r1 = Sqrt[(x - x1)^2 + (y - y1)^2];
r2 = Sqrt[(x - x2)^2 + (y - y2)^2];
r3 = Sqrt[(x - x3)^2 + (y - y3)^2];
Ω = 1/(3*(1 + β*Sqrt[3]))*(β/r0 + 1/r1 + 1/r2 + 1/r3) + 1/2*(x^2 + y^2);
μ = 0.5;
β = 1/μ - 1;
Ωx = D[Ω, x];
Ωy = D[Ω, y];
Then the module FindRoots2D for finding the intersections of two curves
Options[FindRoots2D] = {WorkingPrecision -> 30, MaxRecursion -> 30};
FindRoots2D[funcs_, {x_, a_, b_}, {y_, c_, d_}, opts___] :=
Module[{fZero, seeds, signs, fy},
fy = Compile[{x, y}, Evaluate[funcs[[2]]]];
fZero =
Cases[Normal[
ContourPlot[
funcs[[1]] == 0, {x, a - (b - a)/97, b + (b - a)/103}, {y,
c - (d - c)/98, d + (d - c)/102},
Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]],
Line[z_] :> z, Infinity];
seeds = Flatten[((signs = Sign[Apply[fy, #1, {1}]];
#1[[
1 + Flatten[
Position[Rest[signs*RotateRight[signs]], -1]]]]) &) /@
fZero, 1];
If[seeds == {}, {},
Select[Union[({x, y} /.
Check[FindRoot[{funcs[[1]],
funcs[[2]]}, {x, #1[[1]]}, {y, #1[[2]]},
Evaluate[FilterRules[{opts}, Options[FindRoot]]]], {x ->
b + 999, y -> d + 999}] &) /@ seeds,
SameTest -> (Norm[#1 - #2] < 10^(-8) &)],
a <= #1[[1]] <= b && c <= #1[[2]] <= d &]]]
Then we apply this module to the first-order derivatives
pts = Chop[FindRoots2D[{Ωx, Ωy}, {x, -2, 2}, {y, -2, 2}, PlotPoints -> 100]];
nps = Length[pts];
Print["N = ", nps]
The final output is
cont = ContourPlot[{Ωx == 0, Ωy == 0}, {x, -2, 2}, {y, -2, 2},
ContourShading -> False, ContourStyle -> {{Thick, Darker[Green]}, {Thick, Blue}},
PlotPoints -> 200, PerformanceGoal :> "Speed",
Epilog -> {AbsolutePointSize[8], Red, Point@pts}, ImageSize -> 500]
We see that there are 9 intersections, the coordinates of which are contained in the list pts
. Personally, I have no idea how Mathematica orders the 9 points inside the list.
My question is the following: how can I automatically order (sort) the points, according to the labels 1, 2, ..., 9? I mean, how can I create a new list pts2
with {{x_1, y_1}, {x_2, y_2}, {x_3, y_3}, ..., {x_9, y_9}}?
Any ideas?
Many thanks in advance!
FindRoots2D[]
. You might want to mention you're still using version 9, as most are on 11 now. $\endgroup$