3
$\begingroup$

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]

enter image description here

I believe that the correct answer is 13 equilibrium points.

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

enter image description here

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!

|improve this question
$\endgroup$
  • $\begingroup$ Did you mean to use f[x_, y_] = Ωx instead of f[x_, y_] := Ωx? Didn't read the rest of the post, but this shouldn't work with :=. $\endgroup$ – Szabolcs Feb 27 '17 at 12:42
  • 1
    $\begingroup$ Yes, but it seems that you don't actually use f and g ... Looking at higher resolution contour plots, it appears that there are an infinite number of equilibrium points which all lay on a circle. (This may not be exactly true, but it will throw off th numerics) $\endgroup$ – Szabolcs Feb 27 '17 at 12:54
  • 1
    $\begingroup$ @Szabolcs Well, I still believe that this (infinite points) is actually a numerical artefact. Take $p = 0.2$ and inspect the contours. You will see that there are four points in this area. Now as you tend to $p = 0$ the four points come close and finally for $p = 0$ they disappear. $\endgroup$ – Vaggelis_Z Feb 27 '17 at 13:00
  • 1
    $\begingroup$ Yes, I agree that it must be a numerical artefact. It is also what makes the problem numerically difficult. $\endgroup$ – Szabolcs Feb 27 '17 at 13:04
  • 1
    $\begingroup$ It seems to be possible to brute force the difficult case using Length@FindRoots2D[{f[x, y], g[x, y]}, {x, -2, 2}, {y, -2, 2}, MaxIterations -> 10000, PlotPoints -> 100, MaxRecursion -> 5] (perhaps also using WorkingPrecision -> 40 and setting p=1/100 exactly). This returns 13. It does not prove that this answer is correct, but there are indications that it is, such as your observations about the contour plot when p is increased, looking at the numerical value of the result as WorkingPrecision is increased further, and that the FindRoot errors are gone. $\endgroup$ – Szabolcs Feb 27 '17 at 13:10
2
$\begingroup$

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 *)
|improve this answer
$\endgroup$
  • $\begingroup$ That is simply marvelous! Many thanks :) $\endgroup$ – Vaggelis_Z Feb 27 '17 at 13:14
  • $\begingroup$ This is obviously a quick hack. Maybe a better approach would be to test all the putative roots whether they actually satisfy the equations? $\endgroup$ – Chris K Feb 27 '17 at 13:19
  • $\begingroup$ This looks very interesting. The FindRoots2D module creates the list pts of possible roots. Then how can we easily filter them and see which of them are actually satisfy the system of equations? $\endgroup$ – Vaggelis_Z Feb 27 '17 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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