COMPUTATION OF LONG PERIODIC ORBITS

IN CHAOTIC DYNAMICAL SYSTEMS

Brian A. Coomes, Hüseyin Ko\c cak, and Kenneth J. Palmer

A new method for computing long unstable periodic orbits of chaotic systems is developed. First, a global Newton's method, coupled with back and forth shooting, is used to compute suitable approximate periodic orbits with local errors near the machine precision. Then a periodic shadowing theorem is invoked to establish the existence of true periodic orbits near these approximate periodic orbits.

Periodic orbits are ubiquitous in chaotic dynamical systems. Although they are unstable and hence not visible in numerical simulations, they carry much vital information about the underlying dynamics. While this realization dates back to Poincaré, there have been some recent important efforts to compute quantitative features of specific chaotic systems from the knowledge of unstable periodic orbits~[1,4,5]. For example, for a chaotic map f, if P_N is the number of periodic points of (not necessarily least) period N, then, under suitable conditions, the topological entropy h_T of f is given by the formula h_T = \limsup_{N \to +\infty} (\ln P_N / N). In application of this formula one needs to locate periodic orbits of higher and higher periods to compute better and better approximations to the topological entropy. Not only is an efficient algorithm needed to compute these orbits but also a method to verify that a numerically computed apparent periodic orbit is near a true periodic orbit.

Using shadowing techniques, we established in [2] (see [7,8,9] for related results) a periodic shadowing theorem which guarantees the existence of a true periodic orbit near an approximate periodic orbit with sufficiently small local errors-a pseudo periodic orbit. In Section~2, after recalling the notions of shadowing in the context of periodic orbits, we restate this periodic shadowing theorem. The main contribution in this article is contained in Section~3 where we develop a global Newton's method, coupled with multiple back and forth shooting, for computing pseudo periodic orbits with local errors near the machine precision. In Section~4, we illustrate the combined effectiveness of these numerical and theoretical tools by rigorously establishing the existence of true periodic orbits of several specific systems, including a period 100,003 orbit of the famous chaotic map of Hénon.

In this brief section, we first recall the notions of a pseudo periodic orbit and shadowing in discrete dynamical systems. Then we state a periodic shadowing theorem which guarantees the existence of a true periodic orbit near a pseudo periodic orbit of a C^2 map f:{\Bbb R}^n\to {\Bbb R}^n.

Definition. A sequence $\\bold y _k\}^N_{k=0} is said to be a {\it $\delta$ pseudo periodic orbit/} of~f if} \parallel {\bold y}_{k+1} - f({\bold y}_k) \parallel \le \delta\quad \textor} k=0,\ldots,N-1 and \parallel {\bold y}_0 - f({\bold y}_N) \parallel \le \delta.

Definition. For a given $\varepsilon$, a $\delta$ pseudo periodic orbit $\\bold y _k\}^N_{k=0} is said to be {\it $\varepsilon$-shadowed} by a true periodic orbit if there are points \\bold x}_k\}^N_{k=0} with {\bold x}_{k+1} = f({\bold x}_k) for k= 0, \ldots, N-1 and {\bold x}_0 = f({\bold x}_N) such that} \parallel {\bold x}_k - {\bold y}_k \parallel \le \varepsilon \text{for} k=0,\ldots,N.

If \\bold y}_k\}^N_{k=0} is a \delta pseudo orbit, we define the linear operator L\colon({\Bbb R}^n)^{N+1}\to({\Bbb R}^n)^{N+1} by
\align (L{\bold u})_k &= {\bold u}_{k+1} - Df({\bold y}_{k}){\bold u}_k \quad \text{for} k=0,\ldots,N-1 \\ (L{\bold u})_N &= {\bold u}_{0} - Df({\bold y}_{N}){\bold u}_N,align where {\bold u} = \\bold u}_k\}_{k=0}^N \in ({\Bbb R}^n)^{N+1}.

Periodic Shadowing Theorem. Suppose~$f\colon {\Bbb R}^n\to {\Bbb R}^n$ is a~$C^2$ map. Let $\\bold y _k\}^N_{k=0} be a \delta pseudo periodic orbit for f and suppose that the operator L is invertible. Set} \varepsilon = 2 \parallel L^{-1}\parallel \delta, where the norm of $L^{-1}$ is the operator norm with respect to the supremum norm on $({\Bbb R}^n)^{N+1}$. Now, let M=\sup\bigl\{\parallel D^2f({\bold x})\parallel : {\bold x}\in {\Bbb R}^n, \parallel {\bold x} - {\bold y}_{k}\parallel \le\varepsilon \text{for some} k=0,\ldots ,N\bigr\}.

Then if
2 M \parallel L^{-1}\parallel^2\delta<1,
the pseudo periodic orbit $\\bold x _k\}^N_{k=0} is \varepsilon-shadowed by a unique true periodic orbit \\bold y}k\}^N_{k=0}.}

A proof of this theorem is given in [2]. A major consideration in applying this theorem to a specific pseudo periodic orbit is the verification of the invertibility of the operator L and the estimation of \parallel L^{-1} \parallel. A method for accomplishing this task effectively is expounded also in [2]. However, the important numerical problem still remains: How does one obtain a good pseudo periodic orbit in the first place? This is the question to which we now turn.

We first describe the global Newton's method for a one-dimensional differentiable map f : {\Bbb R} \to {\Bbb R}. We consider the function {\Cal G} : {\Bbb R}^{N+1} \to {\Bbb R}^{N+1} acting on sequences of length N+1 defined by
\aligned [{\Cal G}({\bold x})]_k &={\bold x}_{k+1}-f({\bold x}_k) \quad \text{for}\quad k=0,\ldots ,N-1\\ [{\Cal G}({\bold x})]_N &={\bold x}_0-f({\bold x}_N), aligned \tag1
where {\bold x} = \\bold x}_k\}^N_{k=0} \in {\Bbb R}^{N+1}. Notice that a sequence of N+1 points \\bold x}_k\}^N_{k=0} is a true periodic orbit of (not necessarily least) period N+1 for f if
{\Cal G}({\bold x})={\bold 0}. \tag2

Now, in order to solve this equation numerically, starting with a crude pseudo periodic orbit \\bold x}_k\}^N_{k=0} of period N+1, we apply Newton's method to the function {\Cal G}. After applying Newton's method once, the next approximation to the solution of Eq.(2) will be the sequence {\bold y} = \{{\bold y}_k\}_{k=0}^N given by
{\bold y}={\bold x}-D{\Cal G}({\bold x})^{-1}{\Cal G}({\bold x}).

For brevity, we let
{\bold h}=-D{\Cal G}({\bold x})^{-1}{\Cal G}({\bold x}).

So {\bold h} is the solution of the equation
D{\Cal G}({\bold x}){\bold h}=-{\Cal G}({\bold x}). \tag3

If we put
\aligned {\bold g}_k&=f({\bold x}_k)-{\bold x}_{k+1} \quad \text{for} k=0,\ldots ,N-1\\ {\bold g}_N&=f({\bold x}_N)-{\bold x}_0 aligned

and {\bold h}=\{{\bold h}_k\}^N_{k=0}, we see that Eq.(3) can be written as
\aligned {\bold h}_{k+1}&=Df({\bold x}_k){\bold h}_k+{\bold g}_k \quad \text{for} k=0,\ldots ,N-1\\ {\bold h}_0&=Df({\bold x}_N){\bold h}_N+{\bold g}_N. aligned \tag4

For simplicity of notation, if we set
a_k=Df({\bold x}_k),

then Eq.(4) can be solved using backward shooting as
{\bold h}_0=-(1-a_0^{-1}\ldots a_N^{-1})^{-1} \sum^N_{m=0}a_0^{-1}\ldots a_m^{-1}{\bold g}_m

and
\align {\bold h}_N&=a_N^{-1}({\bold h}_0-{\bold g}_N)\\ {\bold h}_k&=a_k^{-1}({\bold h}_{k+1}-{\bold g}_k) \quad \text{for} k=N-1,\ldots ,1. align

If f is a chaotic map, then |a_k| > 1 for most k and the Newton corrections {\bold h}_k would be computable from the formulas above. Now, the new improved pseudo periodic orbit {\bold y} = \\bold y}_k\}^N_{k=0} is simply given by {\bold y}_k = {\bold x}_k + {\bold h}_k.

The global Newton's method can be generalized for computation of pseudo periodic orbits of higher dimensional chaotic maps. The basic setup is identical to the one just presented in dimension one. To compute a pseudo periodic orbit of period N+1 of a chaotic map f : {\Bbb R}^n \to {\Bbb R}^n, consider the map {\Cal G} : ({\Bbb R}^n)^{N+1} \to ({\Bbb R}^n)^{N+1} defined exactly as before in Eq.(1) and solve Eq.(2). However, care must be taken in the computation of {\bold h} = \\bold h}_k \}_{k=0}^{N} from Eq.(4). The dynamics of higher dimensional chaotic maps usually consist of expanding and contracting components. So we use back and forth shooting to solve Eq.(4). To illustrate the necessary procedures, we now describe the global Newton's method in dimension two.

The contracting and expanding components of a two dimensional map along an approximate periodic orbit {\bold x} = \\bold x}_k\}^N_{k=0} can be identified by triangularizing the matrices Df({\bold x}_k). We begin with an arbitrarily chosen 2\times 2 orthogonal matrix S_0. Then we obtain a sequence \{S_k\}^{N+1}_{k=0} of orthogonal matrices by applying the Gram-Schmidt process to get
Df({\bold x}_k)S_k=S_{k+1}A_k \quad \text{for} k=0,\ldots ,N

where
A_k=\left[\aligned a_k& \quad b_k\ 0& \quad c_k aligned \right] \tag5

is upper triangular. We repeat this procedure with S_{N+1} as the new S_0. After several repetitions we expect the columns of S_{N+1} to be close to those of S_0 except for the signs. Now if we set
{\bold h}_k=S_k \bar {\bold h}_k \quad \text{for} k=0,\ldots ,N
in Eq.(4), we obtain the equations
\aligned \bar{\bold h}_{k+1}&=A_k \bar{\bold h}_k+S_{k+1}^*{\bold g}_k \quad \text{for} k=0,\ldots ,N-\\ \bar{\bold h}_0 &=A_N \bar{\bold h}_N+S_0^*{\bold g}_N, aligned\tag6

where, for k=0,\ldots ,N-1, the matrix A_k=S_{k+1}^*Df({\bold y}_k)S_k is as in Eq.(5) but A_N is taken as the upper triangular part of S_{0}^*Df({\bold y}_N)S_N. Then if we write
\aligned \bar {\bold g}_k&=S_{k+1}^*{\bold g}_k \quad \text{for} k=0,\ldots ,N-1\\ \bar {\bold g}_N&=S_0^*{\bold g}_N, aligned

Eq.(6) can be written in components as
\aligned \bar{\bold h}_{k+1}^{(1)}&=a_k \bar{\bold h}_k^{(1)}+b_k \bar{\bold h}_k^{(2)} +\bar {\bold g}_k^{(1)} \quad \text{for} k=0,\ldots ,N-1 \\ \bar{\bold h}_0^{(1)}&=a_N \bar{\bold h}_N^{(1)}+b_N \bar{\bold h}_N^{(2)}+ \bar {\bold g}_N^{(1)} aligned \tag7

and
\aligned \bar{\bold h}_{k+1}^{(2)}&=c_k \bar{\bold h}_k^{(2)}+\bar {\bold g}_k^{(2)} \quad \text{for} k=0,\ldots ,N-1 \\ \bar{\bold h}_0^{(2)}&=c_N \bar{\bold h}_N^{(2)}+\bar {\bold g}_N^{(2)}.aligned\tag8

Now, for a chaotic map we expect |a_k| > 1 and |c_k| < 1 for most k. We solve Eq.(8) forwards to get
\bar{\bold h}_0^{(2)}=(1-c_N\ldots c_0)^{-1} \sum^N_{m=0}c_N\ldots c_{m+1}\bar {\bold g}_m ^{(2)}

and for k=0,\ldots ,N-1 obtain \bar {\bold h}^{(2)}_k by recursion using the first formula in Eq.(8). We substitute the result in Eq.(7) and solve the resulting equations by backward shooting, as in the scalar case.

To apply the Periodic Shadowing Theorem to the pseudo periodic orbit \\bold y}_k \}_{k=0}^{40}, we must determine rigorous bounds for M and \Vert L^{-1} \Vert. Clearly, here M can be taken as 8. Using the algorithm given in [2], we compute that \Vert L^{-1} \Vert \le 30.478601129435479001. Thus
2 M \Vert L^{-1} \Vert^2 \delta \le 1.969852272634231400 \times 10^{-11}.
Now the theorem implies the existence of a true periodic orbit of period 41 of the Logistic map within
\epsilon = 2 \Vert L^{-1} \Vert \delta \le 8.078833179829062600 \times 10^{-14}
of the pseudo periodic orbit.

Example 2: A periodic orbit of period 100,003 of the Hénon map. Consider the Hénon map f(x_1, x_2) = (1+ x_2 - \alpha x_1^2, \beta x_2) with the parameter values \alpha = 1.4 and \beta = 0.3. Starting with an initial point near the attractor we computed a crude orbit by simply iterating the map 100,002 times. Then we applied the global Newton's method eight times to this crude orbit to obtain a refined pseudo periodic orbit of period 100,003. The refined pseudo periodic orbit had the uniform rigorous local error bound of
\delta=6.5115307479103834 \times 10^{-15}.
In lieu of tabulating the coordinates of the refined pseudo periodic orbit we plot them in Fig.1. It is interesting to compare this picture with the famous picture of the Hénon attractor~[6].

Using the algorithm for estimating \Vert L^{-1}\Vert described in~[2], we find that
2M\Vert L^{-1}\Vert^2 \delta \le 0.000474776016367796,

\midspace{11cm}

\midspace{11cm}

and
\varepsilon=2\Vert L^{-1}\Vert \delta \le 1.486010629876734\times 10^{-9}.
Now the Periodic Shadowing Theorem implies the existence of a true periodic orbit of period 100,003, which is prime, of the Hénon map within a uniform distance of 1.486010629876734\times 10^{-9} of the refined pseudo periodic orbit.

This work was supported in part by the National Science Foundation and was partly written during a visit to CADSEM at Deakin University.

References