# Basins of attraction of equilibrium points

The Henon-Heiles potential is the following

V = 1/2*(x^2 + y^2) - y*(1/3*y^2 - x^2);


which has four equilibrium points

Vx = D[V, x];
Vy = D[V, y];
sol = Solve[{Vx == 0, Vy == 0}, {x, y}]


{{x -> 0, y -> 0}, {x -> 0, y -> 1}, {x -> -(Sqrt[3]/2), y -> -(1/2)}, {x -> Sqrt[3]/2, y -> -(1/2)}}

EQs = {{0, 0}, {0, 1}, {-(Sqrt[3]/2), -(1/2)}, {Sqrt[3]/2, -(1/2)}};
L0 = ListPlot[EQs, PlotStyle -> {Blue, PointSize[0.02]},
Frame -> True, Axes -> False, PlotRange -> {{-2, 2}, {-2, 2}},
AspectRatio -> 1]


In this previous question Post 1 we see how using the Newton's iterative method we can find the basins of attraction in the complex plane.

Let's create a uniform square grid in the (x,y)-plane.

data = Flatten[Table[{i, j}, {i, -2, 2, 0.01}, {j, -2, 2, 0.01}], 1];


Now I would like to do the following: Use again the Newton's iterative process in order to see to which equilibrium point each (x,y) initial condition of the data list is attracted to.

In our case the Newton's algorithm takes the form

$x_n = x_{n-1} - \left(\frac{V_xV_{yy} - V_yV_{xy}}{V_{yy}V_{xx} - V_{xy}^2}\right)_{(x_{n-1},y_{n-1})}$

$y_n = y_{n-1} + \left(\frac{V_xV_{yx} - V_yV_{xx}}{V_{yy}V_{xx} - V_{xy}^2}\right)_{(x_{n-1},y_{n-1})}$,

where $V_x,V_y,V_{xx},V_{xy},V_{yy}$ are the partial derivatives of $V(x,y).$

Any ideas how to obtain this using Mathematica?

NOTE: As far as I know, the task of obtaining the basins of attractions to equilibrium points in real (not complex) functions has not been studied before. So, whoever finds a way to do this using Mathematica will add his/her name (or nickname) in the acknowledgments of the resulted research paper.

• I'm quite sure that the basins of attraction of Newton's method in two real variables have been well studied; there's an elementary discussion on pages 486-487 of Gil Strang's calculus text. If you're interested in studying the basins of attraction of the smooth dynamical system with a given potential, however, that's quite a different thing. Only one of your equilibrium points is attractive for that system. Dec 7 '15 at 15:55
• @MarkMcClure For a given potential, such as that corresponding to the restricted three-body problem, I think there are no studies on the basins of attraction of the equilibrium points. What do you mean that only one equilibrium point is attractive? There are basins of attraction for all of them. Dec 7 '15 at 16:13
• As you have shown, the algebraic system $V_x=0,V_y=0$ has four numeric solutions and, if we apply Newton's method to that system, each of those solutions has an interesting basin of attraction. But that is a different question from whether those equilibrium points are attractive under a system of differential equations like $x'=V_x,y'=V_y$. Dec 7 '15 at 16:29
• @MarkMcClure You are absolutely right! I know that in the Henon-Heiles system only the central point (0,0) is attrcative. But how do we know this? Is there a formula showing that (0,0) is attractive, while the other three are not? Dec 7 '15 at 16:37
• For the system $x'=V_x,y'=V_y$, an equilibrium will be attractive iff it's a local minimum of the potential function. Before proceeding further, though, you should make sure that's the system you're interested in. Perhaps, you're really interested in a second order system, which is likely in celestial mechanics. That's a more complicated question, as you've got to deal with initial conditions on $x'$ and $y'$ as well. Dec 7 '15 at 16:43

I think there maybe be a sign error in the maths, bear with me.

Given your

    V[x_, y_] := 1/2 (x^2 + y^2) - y (1/3 y^2 - x^2)


we can find the partial derivatives:

    Vx = D[V[x, y], x]
Vy = D[V[x, y], y]
Vxx = D[V[x, y], {x, 2}]
Vyy = D[V[x, y], {y, 2}]
Vxy = D[V[x, y], x, y]
Vyx = D[V[x, y], y, x]

sol = Solve[{Vx == 0, Vy == 0}, {x, y}]


Then your function for Newton's algorithm is simply:

    newton[{x_, y_}] := {x, y} - {Simplify[(Vx Vyy - Vy Vxy)/(Vyy Vxx - Vxy^2)], -Simplify[(Vx Vyx - Vy Vxx)/(Vyy Vxx - Vxy^2)]}

newton[{x_, y_}] := {x, y} - {(
x (-1 + 2 x^2 + 2 y + 2 y^2))/(-1 + 4 x^2 +
4 y^2), ((1 + 2 y) (x^2 + (-1 + y) y))/(-1 + 4 x^2 + 4 y^2)}


This is where I think a sign is wrong, note the additional minus before the second Simplify. I believe the minus in your formula for yn should be a plus.

Using the function FixedPoint which you should now be familiar with you can find the coordinate to which points gravitate:

    tab = ParallelTable[FixedPoint[newton, {i, j}], {j, -2, 2, 0.003}, {i, -2, 2, 0.003}];


Then it is simply a case of correlating the coordinate to a solution:

    rules = Rule @@@ Transpose[{sol[[;; , ;; , 2]],Range[Length[sol]]}]
newtab = Map[First@Nearest[rules, #] &, tab, {2}]


Now we simply have to plot this:

    ArrayPlot[newtab, ColorFunction -> "Rainbow", DataReversed -> True]


Which we can see lines up well with your expectations:

    ListPlot[List /@ sol[[;; , ;; , 2]],
PlotStyle -> Evaluate[ColorData["Rainbow"] /@ Rescale[Range[4]]],
Frame -> True, Axes -> False, PlotRange -> {{-2, 2}, {-2, 2}},
AspectRatio -> 1]


Edit: Data as a list use FixedPointList rather than FixedPoint. Then take both the Length (to get the number of iterations required) and the Last to get the final solution.

    tab = ParallelTable[{{i, j},
Sequence @@ Through[{Length, Last}[
FixedPointList[newton, {i, j}]]]}, {j, -2, 2, 0.3}, {i, -2, 2,
0.3}];
tab[[;; , ;; , 3]] = Map[First@Nearest[rules, #[[3]]] &, tab, {2}];

dataList = Flatten /@ Flatten[tab, 1];
TableForm[dataList,
TableHeadings -> {None, {"x0", "y0", "iters", "sol"}}]


Edit2: Second Potential Using

    μ = 1/5;
V[x_, y_] := ((1 - μ)/Sqrt[(x + μ)^2 + y^2]) + μ/Sqrt[(x + μ - 1)^2 + y^2] + 1/2*(x^2 + y^2);


Be sure to obtain the numerical solutions: sol = N /@ Solve[{Vx == 0, Vy == 0}, {x, y}]

And limit the number of iterations:

    tab = ParallelTable[
FixedPoint[newton, {i, j}, 100], {j, -2, 2, 0.013}, {i, -2, 2, 0.013}];


Everything else is the same and gives (now using ListDensityPlot):

• Is this method applicable for any type of potential V(x,y)? I would like to have a data list containing the following: x0,y0,iters, where x0, y0 are the initial conditions, iters is the number of Newton's iteration and n is either 1,2,3, or 4 if we number the four solutions. Is this doable? Dec 7 '15 at 12:05
• Sure, see my edit for a data table. I don't see why the general method shouldn't work for other potentials. Dec 7 '15 at 12:14
• Why you define newton twice? And later on rules can not be recognized. Dec 7 '15 at 12:14
• I only define newton twice to show how I obtained the expressions in the second {} (ie from the Simplify of the other stuff, in reality I would only use the second definition. Sorry I forgot to past the code for rules it's there now. It just provides a mapping from coordinate to solution index. Dec 7 '15 at 12:17
• Many many thanks! Do you want to include your name (or nickname) in the acknowledgments of the resulted paper? Dec 7 '15 at 13:55