How many equilibrium points exist?

Question

Here comes the definition of a simple potential function

Clear["Global`*"];
r1 = Sqrt[(x + μ)^2 + y^2];
r2 = Sqrt[(x - 1 + μ)^2 + y^2];
Ω = (1 - μ)/r1 + μ/r2 - p/(2*c^4)*((1 - μ)^3/r1^3 + μ^3/r2^3) + 1/2*(x^2 + y^2);

μ = 1/2;
p = 0.01;
c = 1;

Ωx = D[Ω, x];
Ωy = D[Ω, y];

f[x_, y_] := Ωx
g[x_, y_] := Ωy

Now I want to determine how many equilibrium points exist when $p = 0.01$ and also when $p = 0.9$. Let's start with $p = 0.01$.

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} /. 
     FindRoot[{funcs[[1]], 
       funcs[[2]]}, {x, #1[[1]]}, {y, #1[[2]]}, 
      Evaluate[FilterRules[{opts}, Options[FindRoot]]]] &) /@ 
  seeds, SameTest -> (Norm[#1 - #2] < 10^(-8) &)], 
 a <= #1[[1]] <= b && c <= #1[[2]] <= d &]]]

pts = FindRoots2D[{Ωx, Ωy}, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 500]
nps = Length[pts];
Print["N = ", nps]

The code reports 27 equilibrium points, while I strongly believe that this result is false.

Let's see what the contours have to say

cont0 = ContourPlot[{Ωx == 0, Ωy == 0}, {x, -2, 2}, {y, -2, 2}, 
ContourShading -> False, ContourStyle -> {{Thick, Darker[Green]}, {Thick, Blue}}, 
PlotPoints -> 150, PerformanceGoal :> "Speed", 
Epilog -> {AbsolutePointSize[8], Red, Point@pts}, ImageSize -> 600]

I believe that the correct answer is 13 equilibrium points.

Similarly for $p = 0.9$ the code reports 3 equilibrium points.

However also in this case I believe that there is only one equilibrium point at (0,0), while the other two reported points are false.

So, my question: What is the correct number of equilibrium points for $p = 0.01$ and for $p = 0.9$? Needless to say that any proposition of another alternative way of correctly computing the number of equilibrium points is more than welcome!

I am using v9.0 of Mathematica.

Many thanks in advance!

Answer 1

I noticed a lot of FindRoot::cvmit errors when you run FindRoots2D. I wrapped a Check around the FindRoot in FindRoots2D that effectively excludes these nonconvergent points (by returning points that are out of bounds by 999 and therefore discarded). This gives the answers you suspect.

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 &]]]

pts = FindRoots2D[{Ωx, Ωy}, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 500];
nps = Length[pts];
Print["N = ", nps]

(* 13 *)