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!
f[x_, y_] = Ωx
instead off[x_, y_] := Ωx
? Didn't read the rest of the post, but this shouldn't work with:=
. $\endgroup$ – Szabolcs Feb 27 '17 at 12:42f
andg
... 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:54Length@FindRoots2D[{f[x, y], g[x, y]}, {x, -2, 2}, {y, -2, 2}, MaxIterations -> 10000, PlotPoints -> 100, MaxRecursion -> 5]
(perhaps also usingWorkingPrecision -> 40
and settingp=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