For the MacKay (Henon action-angle transformation) Map (Mackay 1983):
$x_{n+1}=q-\alpha x_n-x_n^2-y_n \\ y_{n+1}=x_n$
with $q, \alpha \in \mathbb{R}$ real parameters I would like to calculate the points which compose periodic orbits of period $m$. With a collaborator we wrote the following piece of code which calculates those pairs of points.
periodic[period_, initvalx_, initvaly_] :=
With[{q = 2.382163 - 2*0.8390658634208773 - (0.8390658634208773)^2},
Module[{x, y, val, valx, valy, init, eqn},
x = Table[Symbol["x" <> ToString@i], {i, period}];
y = Table[Symbol["y" <> ToString@i], {i, period}];
val = Join[x, y]; init = Array[0 &, {2*period, 2}];
init[[1]][[2]] = initvalx; init[[1 + period]][[2]] = initvaly;
Do[init[[i]][[1]] = x[[i]];
init[[i + period]][[1]] = y[[i]], {i, period}];
Do[init[[i + 1]][[2]] =
q - 2*0.8390658634208773*init[[i]][[2]] -
init[[i + period]][[2]] - init[[i]][[2]]^2;
init[[i + period + 1]][[2]] = init[[i]][[2]], {i, period - 1}];
eqn = Array[0 &, 2*period];
eqn[[1]] =
q - 2*0.8390658634208773*x[[1]] - y[[1]] - x[[1]]^2 - x[[2]] == 0;
eqn[[period + 1]] = x[[1]] - y[[2]] == 0;
Do[eqn[[i + 1]] =
q - 2*0.8390658634208773*x[[i + 1]] - y[[i + 1]] - x[[i + 1]]^2 -
x[[i + 2]] == 0;
eqn[[i + period + 1]] = x[[i + 1]] - y[[i + 2]] == 0, {i,
period - 2}];
eqn[[period]] =
q - 2*0.8390658634208773*x[[period]] - y[[period]] -
x[[period]]^2 - x[[1]] == 0;
eqn[[2*period]] = x[[period]] - y[[1]] == 0;
sol = FindRoot[eqn, init]; val = val /. sol;
valx = Array[0 &, {period}]; valy = Array[0 &, {period}];
Do[valx[[i]] = val[[i]];
valy[[i]] = val[[i + period]], {i, period}];
ListPlot[Table[{valx[[i]], valy[[i]]}, {i, period}]]]]
Although it works very well, we face the following issues:
For a periodic orbit of period $m=5$ we get the following points:
{x1 -> 0.114092, x2 -> -0.31857, x3 -> 0.319023, x4 -> -0.31857, x5 -> 0.114092, y1 -> 0.114092, y2 -> 0.114092, y3 -> -0.31857, y4 -> 0.319023, y5 -> -0.31857}and the correspoding listplot is the following: (cannot upload due to some Imgur error)
but for a period of $m=10$ for the same initial conditions we get:
{x1 -> 0.244737, x2 -> -0.375863, x3 -> 0.244737, x4 -> -0.0947348,
x5 -> -0.0947348, x6 -> 0.244737, x7 -> -0.375863, x8 -> 0.244737,
x9 -> -0.0947348, x10 -> -0.0947348, y1 -> -0.0947348,
y2 -> 0.244737, y3 -> -0.375863, y4 -> 0.244737, y5 -> -0.0947348,
y6 -> -0.0947348, y7 -> 0.244737, y8 -> -0.375863, y9 -> 0.244737,
y10 -> -0.0947348}
Which are the same points repeating themselves after $(x_5,y_5)$ with list plot of five and not ten points, as it should be.
I am guessing this has to do with FindRoot function of Mathematica, as apparently the corresponding non-linear system of equations ($a \times 5$ in this particular case) accepts this set of solutions multiple times. Therefore that is not an actual period $m=10$ but a period $m=5$.
How could it be possible to tell Mathematica to not accept this as a solution and show only the period $m=5$? In this case it is easy to see by eye, but when I search for e.g. period $m=232$ this is not easy to deal with. I am thinking it would require a logical condition, but I am not yet able to implement something like that.
- Also, is it possible to tell Mathematica to show the results as pairs of $(x_i,y_i)$ and not in the way it does (first all $x_i$ and then all $y_i$)? Because this way it is hard to connect each pair, especially when the number of periodic points is large.
Thanks!
