A Fast Method to Find Periodic Orbits in Chaotic Attractors With Applications to the Rössler System

A systematic method to find short unstable periodic orbits embedded in chaotic attractors is proposed. The method is based on the construction of symbolic representations of trajectories, which are used to limit the number of symbol sequences considered to find periodic orbits with a given period and also to find candidates of periodic orbits. The Newton method is applied to find accurate positions of periodic orbits. Using the Rössler system as an example, it is shown that the proposed method outperforms existing methods in terms of the number of periodic orbits found and the computation time.


I. INTRODUCTION
U NDER certain assumptions a chaotic attractor is densely filled by unstable periodic orbits. Low period orbits give good approximations to the mean properties of chaotic trajectories [1], [2]. In this brief, we consider the problem of finding unstable periodic orbits embedded in chaotic attractors. For discrete-time dynamical systems, we assume that the system is defined by the map P. For continuous-time dynamical systems the problem under study is first reduced to the discrete time by selecting a return map P. In both cases, the goal is to find periodic orbits of P with periods p ≤ p max .
The standard method to find period-p orbits of a nonlinear map P (or fixed points of P p ) is based on the Newton method to find zeros of a nonlinear function. In this approach one constructs a function F whose zeros correspond to periodic orbits of P, selects an initial guess for the position of a periodic orbit and applies the Newton iteration in the hope that the method converges. The function F may be defined as F(u) k = w (k+1) mod p − P(w k ), k = 0, 1, . . . , p − 1. (1) where u = (w 0 , w 1 , . . . , w p−1 ). It is clear that u is a zero of F if and only if (w 0 , w 1 , . . . , w p−1 ) is a periodic orbit of P.
The Newton method is an iterative process in which successive approximations u (k) after the initial guess u (0) are obtained by applying the formula u (k+1) = u (k) − (F (u (k) )) −1 F(u (k) ). ( The convergence of the Newton method depends on the quality of the initial guess. If u (0) is close to the true position of a zero of F then the Newton method converges very fast. For a general nonlinear map P one can find periodic orbits by using different initial guesses. Common choices are to select them randomly or to locate them on a uniform grid.
To find periodic orbits embedded in a chaotic attractor usually a better choice is to use the trajectory monitoring approach [1], where one looks for pseudo periodic orbits in a long trajectory Pseudo periodic orbit are used as initial guesses for the Newton method applied to the map (1). This method was successfully applied to find short periodic orbit for the Chua's circuit with a cubic nonlinearity [3], [4]. The main drawback of this method is that a given orbit may be found several times, which slows down the search process. For specific classes of discrete and continuous dynamical systems more efficient methods have been proposed [5], [6], [7], [8].
In this brief, we propose a systematic method to find short unstable periodic orbits embedded in a chaotic attractor. The method is based on the construction of a symbolic representation of trajectories belonging to the attractor and can be applied to a wide class of dynamical systems. As an example we study the existence of periodic orbits for the Rössler system [9]. Studying dynamical phenomena including the existence of periodic orbits for the Rössler system is an active field of research. Low order periodic orbits for the Rössler system are extracted and encoding of periodic orbits by symbolic dynamics is studied in [10]. Rigorous analysis of the existence of short periodic orbits is carried out in [11]. Different routes to chaos and different kinds of chaotic attractors existing for this system are studied in [12]. Averaged properties of unstable periodic orbits existing in the Rössler attractor are studied in [2]. Unbounded dynamics and the homoclinic chaos in the Rössler model is observed in [13], [14]. The existence of chaos and hyperchaos in the 4D Rössler system is proved in [15]. The origin of homoclinic chaos in the classical Rössler model is studied in [14]. Conditions for the existence of a periodic solution bifurcating from the zero-Hopf equilibrium are formulated in [16]. The existence of infinitely many periodic orbits in the Rössler system for two sets of parameter values is proved in [17]. Practical applications of the Rössler system are considered in [18]. Nonlinear dynamics in weakly-synchronized This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Rössler system is investigated in [19]. Information transfer with respect to relative entropy in the Rössler system is studied in [20].
The remainder of this brief is organized as follows. In Section II a systematic method to find unstable periodic orbits is described in detail. The results obtained by applying this method to the analysis of the Rössler system are presented in Section III. The classical case with a symbolic dynamics defined on two symbols and a more complex case with a symbolic dynamics using three symbols are considered. For both cases a symbolic representation of trajectories is constructed and the proposed method is used to find low period orbits. The last section concludes the study.

II. FINDING UNSTABLE PERIODIC ORBITS
In this section a systematic procedure to find unstable periodic orbits embedded in a chaotic attractor is presented. The map P with a chaotic attractor is either a map defining a discrete time system or a return map for a continuous time system. The proposed procedure to find period-p orbits of P with p ∈ {1, 2, . . . , p max } consists of the following steps: 1) define a symbolic representation of trajectories, 2) find short forbidden symbol sequences, 3) generate non-forbidden symbol sequences with periods p ∈ {1, 2, . . . , p max }, 4) for each non-forbidden sequence s find a candidateũ for the position of the corresponding periodic orbit based on symbolic dynamic information, 5) for each candidateũ apply the Newton operator to find the positionū of the orbit, the convergence of the Newton method to the positionū with the symbol sequence s confirms the existence of a periodic orbit, 6) verify that all periodic orbits found are different. In the following part of this section, the steps listed above are described in detail.
The proposed method can be applied to systems for which a symbolic representation of trajectories exists. Defining a symbolic representation of trajectories is perhaps the most difficult part of the procedure. We seek for a splitting of the state space into m regions R 1 , R 2 , . . . , R m such that there is a one-to-one correspondence between trajectories and symbol sequences (see also [10], [21]). For a trajectory (w i ) the corresponding symbol sequence (σ i ) is defined by the condition w i ∈ R σ i . Since we are interested in finding short periodic orbits embedded in the attractor we can relax the condition on the correspondence between trajectories and symbol sequences and require that different symbol sequences are observed not for all trajectories but only for periodic orbits of interest.
Here, we propose a simple method to define symbolic dynamics for narrow attractors. The method is based on attractor parametrization. For simplicity, we assume that the attractor can be parameterized using one of the system variables, for example y. The construction of a symbolic dynamics starts by computing a long trajectory of P. If the attractor is narrow then the plot of y k+1 versus y k resembles a plot of a one-dimensional map. Positions of extreme values in this plot are used to divide the state space into regions corresponding to different symbols. Two examples of this procedure for the case of symbolic dynamics with 2 and 3 symbols are presented in Section III. Using the Rössler system as an example, we show that the proposed method works fine in case of narrow attractors. For systems with more complex attractors other methods to define symbolic dynamics on the attractor may be used.
Let us assume that each periodic orbit of P with the period p ≤ p max has a distinct symbolic representation with m symbols. The number of symbol sequences with the length p is m p . Symbol sequences s = (s 0 , s 1 , . . . , Points belonging to a given periodic orbit have symbolic representations which are not cyclically different. The period of a symbol sequence and corresponds to a period-k orbit. To find all period-p orbits of P it is sufficient to consider S all (p) Considering all cyclically different symbol sequences with the period p may be infeasible for large p due to the fast growth of S all (p) (it grows with p roughly as fast as p −1 m p ). Usually, it is possible to limit the number of symbol sequences which have to be considered. The idea is to identify short forbidden sequences and skip symbol sequences containing forbidden sequences. Forbidden sequences can be identified by monitoring a long trajectory (w i ) N−1 i=0 . First, we select the maximum length l max of forbidden sequences which we search for. Next, for each l = 2, 3, . . . , l max we find symbol sequences of the length l which are not present in the trajectory. Such sequences which do not contain forbidden sequences of a shorter length are added to the list of forbidden sequences. The idea of using forbidden sequences is based on the assumption that if a given short symbol sequence is not present in a long trajectory then it is perhaps not admissible by the dynamical system. The number of symbol sequences of the length l is m l . It follows that the length N of a trajectory which has to be considered to correctly detect all forbidden sequences with the length l ≤ l max grows exponentially with l max . One should also note that short forbidden sequences eliminate much more symbol sequences than the long ones. Therefore, it is usually sufficient to select a small value of l max , for example l max ≤ 10. Now, we describe how to find a periodic orbit u = (w 0 ,w 1 , . . . ,w p−1 ) with the symbol sequence s =  (s 0 , s 1 , . . . , s p−1 ). To obtain a good initial guess for the Newton method we use a trajectory (w i ) N−1 i=0 with the symbol sequence σ = (σ i ) N−1 i=0 . Let us select h ≥ 0. For k = 0, 1, . . . , p − 1 we find a position j k in the sequence σ such that s (k+i) mod p = σ j k +i for all i = −h, −h + 1, . . . , i max and that i max is as large as possible. The point u = (w 0 ,w 1 , . . . ,w p−1) = (w j 0 , w j 1 , . . . , w j p−1 ) is used as an initial guess for the Newton method to find the periodic orbit with the symbol sequence s = (s 0 , s 1 , . . . , s p−1 ).
To explain why this method works let us first assume that h = 0. The sequence (σ j k +i ) N−j k i=0 matches the sequence Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. (s (k+i) mod p ) ∞ i=0 for the first i max + 1 elements, which means that trajectories starting at w j k and at the periodic pointw k have the same future symbol sequences up to the position i max . The larger value of i max , the smaller is the distance between w j k andw k . To make sure that we obtain long matched sequences (large i max ) for all admissible periodic orbit we have to consider a sufficiently long trajectory (w i ) N−1 i=0 . Selecting h = 0 works fine when the attractor is very narrow and the plot of y k+1 versus y k is close to a plot of a one-dimensional map. If this is not the case and branches in the plot of y k+1 versus y k are far apart then one should select h > 0. This way we force that symbols of preimages of w j k andw k agree, which ensures that we select a correct branch of the attractor. In the examples considered in Section III it is sufficient to use h = 1. For h = 0 some of the existing periodic orbits are not found due to the lack of convergence of the Newton method.
In the final step, the Newton method with the initial guess u (0) =ũ is applied to find the position of the periodic orbit. Convergence of the Newton method to a periodic orbitū with the symbol sequence s confirms the existence of the periodic orbit in question.

III. SHORT PERIODIC ORBITS FOR THE RÖSSLER SYSTEM
As an example let consider the Rössler system [9] defined by the following set of ordinary differential equations: We study the system with two sets of parameter values (a, b) = (5.7, 0.2) and (a, b) = (5.7, 0.32). Example trajectories are shown in Fig. 1. We choose the return map P defined by the plane = {v = (x, y, z) ∈ R 3 : x = 0,ẋ > 0}. On the plane we use the local coordinate system w = (y, z) ∈ R 2 . The return map is defined as follows. Let us select w = (y, z) ∈ R 2 such that v = (0, y, z) ∈ . The image of w = (y, z) ∈ R 2 under the map P is P ((y, z) Fig. 2(a). Fig. 2(b) shows a plot of y n+1 versus y n .  The plot presented in Fig. 2(b) is visually indistinguishable from a plot of a one-dimensional map with a single minimum. It is therefore natural to construct a symbolic dynamics with two symbols using the coordinates of the minimum to define the splitting of the state space into regions corresponding to different symbols. From a trajectory (w k ) N−1 k=0 with N = 10 6 we select the point w m 0 with the minimum value of the y coordinate (y m 0 ≈ −10.57497). The value y m 0 −1 ≈ −6.73872 is selected to define the symbolic dynamics with two symbols. The symbol associated with the point w = (y, z) is The splitting of the state space into regions corresponding to different symbols is shown in Fig. 3(b) with a red vertical line.
In the second part of the procedure forbidden sequences are found. This is done by scanning a trajectory (w k ) N−1 k=0 with N = 10 7 and identifying symbol sequences which are not present in this trajectory. The search for forbidden sequences is limited to sequences of the length l ≤ 14. Seven forbidden sequences are found: (111), (1100), (110101), (11010001),  (1101000001), (1101000000011), (11010000000100).
Next, all non-forbidden cyclically different symbol sequences with periods p = 1, 2, . . . , p max = 30 are generated. In Table I, we report the numbers S all (p) of all symbol sequences and the numbers S nf (p) of non-forbidden symbol sequences with the period p. One can see that S all (p) grows much faster than S nf (p). It follows that removing forbidden sequences is an important part of the procedure since it permits elimination of many sequences and thus reduces the computation time needed to find all periodic orbits with a given period.
In the final step of the search procedure, for each nonforbidden sequence a candidate for the position of the periodic orbit is constructed based on the symbolic dynamics information in the way described in Section II and the Newton method is applied to find the position of the corresponding periodic orbit. The numbers P sd (p) of period-p orbits found are given in Table I. The total number of periodic orbits found is 386887. The total computation time needed to find all period-p orbits with p ≤ p max is 136 seconds for p max = 20, 16 minutes for p max = 25 and 3.5 hours for p max = 30. The most time consuming parts of the search procedure are steps 4 and 5 (see Section II), which for p max = 30 require 39.5% and 59.7% of the total computation time, respectively. Other steps are fast and require less than 1% of the total computation time.
The number of periodic orbits found with the period p ≤ 20 is 3655. This is in agreement with the results reported in [11], where all periodic orbits with periods p ≤ 20 were found using rigorous interval arithmetic based calculations. This consistence confirms that the proposed method is successful in finding all short periodic orbits.
For comparison, we also report results obtained using the trajectory monitoring approach where initial guesses for the Newton method are δ pseudo periodic orbits with δ = 0.05 found in a chaotic trajectory of the length N = 2·10 7 . The total computation time is 40 hours. The numbers P tr (p) of periodp orbits found are presented in Table I. The total number of periodic orbits found is P tr (≤ 30) = 102813, which is 26% of the number P sd (≤ 30) = 391018 of orbits found using the symbolic dynamics based method. This is in spite of the fact that the computation time is several times longer.
The numbers of periodic orbits found using both methods are the same for p ≤ 18. The length n tr (p) of a trajectory needed to find all period-p orbits using the trajectory monitoring approach grows exponentially with p and exceeds N = 2 · 10 7 for p > 18. The inferior performance of the trajectory monitoring approach is a consequence of multiple application of the Newton method to find a given periodic orbit, whereas in the symbolic dynamics based approach each symbol sequence is considered only once.
B. The Case a = 5.7, b = 0.32 Let us now consider the case a = 5.7, b = 0.32. A trajectory of the return map P is plotted in Fig. 3(a). The plot of y n+1 versus y n presented in Fig. 3 The division of the state space into regions corresponding to different symbols is shown in Fig. 3(b). In the region s = 1 one can see two branches of the attractor.
Scanning a trajectory with the length N = 10 7 , we find 11 forbidden sequences with the length l ≤ l max = 8:  Table II. Due to a much larger number of short periodic orbits the search for periodic orbits is limited to periodic orbits with periods p ≤ p max = 20. The numbers P sd (p) of period-p orbits found are presented in Table II. The total number of periodic orbits found is P sd (≤ 20) = 985572. The total time needed to find all periodic orbits with periods p ≤ p max is 9 minutes for p max = 15 and 8.2 hours for p max = 20.
For comparison, the trajectory monitoring approach with a trajectory of the length N = 4 · 10 7 is applied to find periodic orbits with periods p ≤ 20. The total computation time is 19.5 hours. The number P tr (p) of periodic orbits found using this method is the same as for the symbolic dynamics based method for p ≤ 10. The total number of periodic orbits found using the trajectory monitoring approach is P tr (≤ 20) = 246739, which is 25% of the number P sd (≤ 20) = 985572.

C. Topological Entropy
Based on the number of periodic orbits, one may estimate the value of the topological entropy of the return map  P using the formula h(P) = lim sup p→∞ p −1 log Q p , where Q p is the number of fixed points of P p [21]. The estimates h p = p −1 log Q p obtained from the results presented in Tables I and II are plotted in Fig. 4. The convergence of the estimates h p is another indication that the proposed method is successful in finding the majority of short periodic orbits embedded in a chaotic attractor.
From the results obtained, it follows that the topological entropy of the return map P is approximately h(P) ≈ 0.5098 for the classical case and h(P) ≈ 0.808 for the second case. We may conclude that the dynamics for the second case is topologically much more complex than for the classical case.

IV. CONCLUSION
A systematic method to find unstable periodic orbits embedded in chaotic attractors was proposed. The method is based on symbolic representation of trajectories belonging to the attractor. Symbolic representations are used to find non-forbidden sequences and to construct initial guesses for the Newton method. The Newton method is applied to find accurate positions of unstable periodic orbits. The method was applied to the analysis of the Rössler system. Two cases were considered with symbolic representations defined on two and three symbols, respectively. It was confirmed that the method is capable of finding the majority (perhaps all) of short periodic orbits in a relatively short time. It was shown that the proposed method outperforms the trajectory monitoring based approach both in terms of the number of periodic orbits found and the computation time. The proposed method is applicable to a wide class of discrete and continuous dynamical system for which a symbolic representation of trajectories can be defined.