# Numerical search for periodic orbits

Suppose that $\phi$ is a map defined in all of $R^n$, with image also in $R^n$. Is there a known good algorithm to find the periodic points of $\phi$ of order $m$, (or fixed points of $\phi^m$) besides:

1. Applying Newton's method in some (random) initial $m$ points $x_0,x_1, ..., x_{m-1}$ and seeking for the solutions of the system $\phi(x_k)=x_{k+1}$;
2. Or applying Newton's method in a single (random) initial point $x$, and looking for the solution of $\phi^m(x)=x$ (thus having to compute the product of the Jacobians of $\phi$ on its first $m$ powers).

-
You could hitch your Newton iteration to a cycle finding algorithm... – J. M. Apr 8 '13 at 0:12

If you're specifically interested in attractive orbits then, of course, you can just iterate the function. A major theorem in one-dimensional dynamics states that attractive orbits always attract at least one critical point. This means that there's only a few initial seeds you need to run your test from. I'm not sure if there's an analog to this in higher dimensional dynamics, though. Furthermore, it wont't help for non-attractive orbits anyway.

In the more general case of (potentially) non-attractive orbits, the first thing to note is that it's probably not a good idea to try to find fixed points of $f^{m}$, as you seem to know. The reason is that Newton's method, like many numerical algorithms, is fairly sensitive to initial conditions. To get around this, suppose we are looking for a point of period $m$ and define a map $F:{\mathbb R}^{n m} \rightarrow {\mathbb R}^{n m}$ by

$$F(x_1,x_2,\ldots,x_m) = (f(x_n), f(x_1), \ldots, f(x_{n-1})).$$

Note that each $x_i$ here is a vector in $\mathbb R^n$ and that an orbit of period $m$ for $f$ forms a vector that is a fixed point of $F$. Shifting to a higher dimensional map like this lightens the problem of sensitivity.

Next, if you want all periodic points, there are algorithms that can simultaneously find all solutions to many types of equations and systems of equations. One simple technique, that works for a univariate polynomial, is to realize the polynomial as the characteristic polynomial of a matrix and then use techniques from numerical linear algebra to find the eigenvalues. A system of polynomials or even algebraic equations can be reduced using Groebner basis techniques, which is somewhat analogous to Gaussian elimination; the eigen-approach can then be used on this system.

Of course these techniques and others are implemented in many modern numerical tools. My personal favorite is Mathematica, which can also solve many systems of transcendental equations, if you're will to restrict the domain of search.

Example

As an example, let's explore the chaotic Henon map defined by

$$f(x,y) = (y+1-1.4x^2, 0.3x).$$

This is a fairly well known 2D chaotic map with a strange attractor that looks like so:

f[{x_, y_}] = {y + 1 - 1.4 x^2, 0.3 x};
hpic = ListPlot[NestList[f, {0, 0}, 1000]]


Note that I'll work this example out in Mathematica but will try to keep the code to a minimum. I think the code for the important steps is fairly clear. For example, we clearly see the definition of $f$ in the code above. And, again, there is certainly other numerical software out there that could do this.

Next, we define the function $F$:

F[xys:{{_, _} ..}] := f /@ RotateRight[xys];


The {{_,_}..} business defines $F$ for lists of pairs; given such a list, the list is rotated to the right and then $f$ is mapped onto the result. Now, suppose we want all orbits of period 7.

vars = Table[{x[k], y[k]}, {k, 1, 7}];
orbits = Select[Chop[
NSolve[F[vars] == vars, Flatten[vars]]],
FreeQ[#, _Complex] &];
firstOrbit = First[orbits]

(* Out: {x[1] -> -0.285392, y[1] -> -0.326186, x[2] -> 0.559786,
y[2] -> -0.0856176, x[3] -> 0.475678, y[3] -> 0.167936,
x[4] -> 0.851159, y[4] -> 0.142703, x[5] -> 0.128443,
y[5] -> 0.255348, x[6] -> 1.23225, y[6] -> 0.038533,
x[7] -> -1.08729, y[7] -> 0.369675}
*)


The key portion is simply the NSolve step, which find the orbits. Let check the one orbit displayed.

{x0, y0} = {-0.2853920648214598, -0.3261858137547133};
seventhIterate = Nest[f, {x0, y0}, 7]
seventhIterate - {x0, y0}

(* Out:
{-0.285392, -0.326186} *)
{1.04361*10^-14, 6.10623*10^-16}
*)


Not bad! Here's what the orbit looks like on the attractor itself:

There's a lot to do here, though. Notably, higher order iterates are going to be harder. Sometimes a homotopy approach can help there. I believe that you need a bit more knowledge and cleverness concerning the system, though. I sketched a homotopy approach to a different question on mathematica.stackexchange.com: http://mathematica.stackexchange.com/questions/19480/

-
Thanks for the answer! I'm not acquainted with Mathematica, though I'm starting to think I should be working with it. – paoloff Apr 8 '13 at 21:28
@Igorevitch No problem! Hope it helps. – Mark McClure Apr 8 '13 at 21:44

A software toolbox like GAIO could help.

-