# Rotation Number of the Henon Map

I would like to plot the rotation number of the last KAM of the Henon Map over several inital conditions. This means that for the 2D map:

$\qquad x_{i+1}=x_i \cos \alpha - \sin \alpha (y_i-x_i^2 ),\\ \qquad y_{i+1}=x_i \sin \alpha + \cos \alpha (y_i-x_i^2)$

I would like to calculate the number: $\sigma = \lim_{N \to \infty} \frac{\sum_{i}^N \sigma _i}{N}$ So I have written the following:

x = Array[0 &, {10000}];
y = Array[0 &, {10000}];
σ = Array[0 &, {10000}];
S = 0
x[] = 0.1;
y[] = x[];
α = 0.4*Pi;
Do[x[[i + 1]] = x[[i]]*Cos[α] - Sin[α]*(y[[i]] - x[[i]]^2);
y[[i + 1]] =
x[[i]]*Sin[α] + Cos[α]*(y[[i]] - x[[i]]^2);σ[[i]] =
ArcCos[y[[i + 1]]/(Sqrt[x[[i + 1]]^2 + y[[i + 1]]^2])] -
ArcCos[y[[i]]/(Sqrt[x[[i]]^2 + y[[i]]^2])], {i, 9999}]
ListPlot[Table[{x[[i]], y[[i]]}, {i, 1, 10000}]]
ListPlot[Table[σ[[i]], {i, 1, 9999}]]
Do[S = (S + σ[[i]])/10^4, {i, 9999}]


which works just fine for one initial condition $x_1=y_1$. How can I now add more initial conditions in the loop and calculate the $\sigma$ for each one of these so that I make a plot $\sigma$ versus, lets say $x_i$ (or $y_i$ does not matter since in my case they are equal)?

It is important that all the pairs of initial conditions are between 0.1 and 0.7

• Why not use S = Mean[σ] in place of Do[S = (S + σ[[i]])/10^4, {i, 9999}]? Unless I misunderstand your purpose, your calculation of S is not even correct. (you shouldn't be dividing by 10^4 at every iteration) – m_goldberg Feb 20 '17 at 15:29

This is not the best way to fix your code. I don't have time to rewrite it in functional style. But this will let you generate your curves and limit approximations.

The key to running your process with different initial values (and also allowing for specifying the number of iterations) is to define a function taking those parameters as arguments.

henon[init_, n_] :=
With[{α = 0.4*Pi},
Module[{x, y, σ},
x = ConstantArray[0, n];
y = ConstantArray[0, n];
σ = ConstantArray[0, n];
x[] = y[] = init;
σ[] = Null;
Do[
x[[i + 1]] = x[[i]]*Cos[α] - Sin[α]*(y[[i]] - x[[i]]^2);
y[[i + 1]] = x[[i]]*Sin[α] + Cos[α]*(y[[i]] - x[[i]]^2);
σ[[i + 1]] =
ArcCos[y[[i + 1]]/(Sqrt[x[[i + 1]]^2 + y[[i + 1]]^2])] -
ArcCos[y[[i]]/(Sqrt[x[[i]]^2 + y[[i]]^2])],
{i, n - 1}];
{x, y, σ}]]


With this definition, your example can be reproduced as follows.

The xy plot.

h[.1] = henon[.1, 1000];

ListPlot[Transpose[h[.1][[1 ;; 2]]]] The σ plot

ListPlot[h[.1][]] S = Mean[Rest[h[.1][]]]


-0.000639365

To run the process for many initial values is easy.

With[{n = 1000}, Do[h[i] = henon[i, n], {i, .1, .6, .1}]]


The approximations to the σ limits is given by

Table[Mean[Rest[h[i][]]], {i, .1, .6, .1}]

{-0.000639365, -0.000273044, 0.0000968227, 0.000195651,
-0.000465987, -0.000645517}


Note: I do not show h[.7] because that evaluation produces numeric errors. I do not have time to explore the problem which is probably more mathematical than Mathematica-ish.

### Update

This is to answer an issue raised by the OP in a comment made below.

The last Do-loop above produces six indexed variables h[.1], ..., h[.6], each storing a 3 x n array containing all the data needed to plot x, y, and σ in any combination.

To plot σ against x for the initial condition x = .3, you simply need to write

ListPlot[Transpose[h[.3][[{1, 3}]]]]


because data for initial condition x = .3 is assigned to h[.3]; the x-values are in row 1 and the σ-values are in row 3. In this case, you would get • Thank you for taking the time to write this piece of code. I will try it and let you know how it works. Btw, the overflow at $.7$ is expected, since this is a chaotic map. Thanks! – Mitscaype Feb 20 '17 at 18:57
• Ok, everything seems fine, but still how can I plot the $\sigma$ versus each $x_i$ initial condition? – Mitscaype Feb 20 '17 at 20:27
• @Mitscaype. The vector h[i][]] contains the σ data for initial condition $x_i$. Look at how I made the σ plot for $x_1$. – m_goldberg Feb 20 '17 at 21:51
• Yeap I can see this. But I need to plot the sigma versus not only $x_1$ but all $x_i$. For every initial condition I will get a different sigma, like you did for $x_1=0.1$. So what I really want is a plot of those sigmas on the vertical axes versus each $x_i$ on the horizontal axes. But I do not know how to extract those $x_i$ from your code :/ – Mitscaype Feb 20 '17 at 21:57
• @Mitscaype. See the update I added to my answer. – m_goldberg Feb 20 '17 at 22:23