3
$\begingroup$

Let's define a simple system of non-linear equations regarding magnetic dipoles:

Clear["Global`*"];
r1 = Sqrt[(x - 0.5)^2 + y^2];
r2 = Sqrt[(x + 0.5)^2 + y^2];
p1 = 1/r1^3 + λ/r2^3;
q1 = 1/2*(1/r1^3 - λ/r2^3);
Ax = -y*p1;
Ay = x*p1 - q1;
Ω = ω*(x*Ay - y*Ax) + ω^2/2*(x^2 + y^2);

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

ω = 1;
λ = 1;

My goal is to solve the system of nonlinear equations $\Omega_x(x,y) = 0, \Omega_y(x,y) = 0$ and determine the total number of equilibrium points.The total number of the equilibrium points strongly depends on the value of $\lambda$ and it is equal to 3, 5 or 7.

For this purpose, I used the accepted (first) answer of this question. Note that all the other suggested answers also fail to properly work for this particular system of nonlinear equations.

Options[FindRoots2D] = {PlotPoints -> Automatic, 
        MaxRecursion -> Automatic};

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^(-6) &)], 
 a <= #1[[1]] <= b && c <= #1[[2]] <= d &]]]

So

pts = FindRoots2D[{Ωx, Ωy}, {x, -5, 5}, {y, -5, 5}];
ContourPlot[{Ωx == 0, Ωy == 0}, {x, -2, 2}, {y, -2, 2}, 
ContourShading -> False, PlotPoints -> 100, 
PerformanceGoal :> "Quality", 
Epilog -> {AbsolutePointSize[6], Red, Point@pts}]

And the output is the following

enter image description here

As you can see, the code finds five equilibrium points, which is of course wrong! There should be only three equilibrium points located on the $x$ axis, all with $y = 0$.

My question: Is there an efficient way to improve the suggested piece of code so to obtain the correct number of equilibrium points? Of course, any other alternative method for obtaining the equilibrium points will be highly appreciated.

I use version 9.0 of MMA.

$\endgroup$
  • $\begingroup$ FindRoot[] is complaining about accuracy. If you feed FindRoot[] with initial values too far from any actual roots it can give fictional results. $\endgroup$ – Feyre Aug 1 '16 at 9:52
  • $\begingroup$ @Feyre But the problem is that I do not know beforehand the positions of the roots. Sometimes all of them lie on the x axis, while some other times roots have non zero values of y. $\endgroup$ – Vaggelis_Z Aug 1 '16 at 9:56
  • $\begingroup$ Perhaps a modification of the answer here? mathematica.stackexchange.com/questions/54438/… where you ignore any results that give warnings? $\endgroup$ – Feyre Aug 1 '16 at 9:59
  • $\begingroup$ If that solution doesn't work in your case, why not try the other solutions in the thread you linked to? $\endgroup$ – J. M. will be back soon Aug 1 '16 at 12:21
  • $\begingroup$ @J.M. Other solutions present other types of malfunctions for this particular set of equations. $\endgroup$ – Vaggelis_Z Aug 1 '16 at 12:23
3
$\begingroup$

As we figured out in the comments section above, adding PlotPoints->100 as an option to FindRoots2D fixes the problem.

pts = FindRoots2D[{Ωx, Ωy}, {x, -5, 5}, {y, -5, 5}, PlotPoints -> 100]
ContourPlot[{Ωx == 0, Ωy == 0}, {x, -2, 2}, {y, -2, 2}, ContourShading -> False,
  PlotPoints -> 100, PerformanceGoal :> "Quality", Epilog -> {AbsolutePointSize[6], Red, Point@pts}]

enter image description here

$\endgroup$
  • $\begingroup$ Many thanks, once more! $\endgroup$ – Vaggelis_Z Aug 1 '16 at 17:35

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.