Convergence Time towards Periodic Orbits in Discrete Dynamical Systems

We investigate the convergence towards periodic orbits in discrete dynamical systems. We examine the probability that a randomly chosen point converges to a particular neighborhood of a periodic orbit in a fixed number of iterations, and we use linearized equations to examine the evolution near that neighborhood. The underlying idea is that points of stable periodic orbit are associated with intervals. We state and prove a theorem that details what regions of phase space are mapped into these intervals (once they are known) and how many iterations are required to get there. We also construct algorithms that allow our theoretical results to be implemented successfully in practice.


Introduction
Periodic orbits are the most basic oscillations of nonlinear systems, and they also underlie extraordinarily complicated recurrent dynamics such as chaos [1][2][3][4][5]. Moreover, they occur ubiquitously in applications throughout the sciences and engineering. It is thus important to develop a deep understanding of periodic dynamics.
It is important and common to question how long it takes a point in phase space to reach a stable periodic orbit from an arbitrary initial condition. When studying synchronization and other forms of collective behavior, it is crucial to examine not only the existence of stable periodic orbits but also the time that it takes to converge to such dynamics in both natural and human-designed systems [6][7][8]. For example, it is desirable to know how long it will take an engineered system that starts from an arbitrary initial condition to achieve the regular motion at which it is designed to work [9,10]. A system can also be perturbed from regular motion by accident, and it is important to estimate how long it will take to return to regular dynamics. Similar questions arise in physics [11,12], biology [6,13,14], and many other areas. It is also important to consider the time to synchronize networks [15][16][17] and to examine the convergence properties of algorithms for finding periodic orbits [2,18].
To study the problem of convergence time to periodic orbits, let's first consider the Hartman-Grobman Theorem [19,20], which states that the flow of a dynamical system (i.e., a vector field) near a hyperbolic equilibrium point is topologically equivalent to the flow of its linearization near this equilibrium point. If all of the eigenvalues of the Jacobian matrix evaluated at an equilibrium have negative real parts, then this equilibrium point is reached exponentially fast when one is in a small neighborhood of it. To determine convergence time to a hyperbolic equilibrium, we thus need to calculate how long it takes to reach a neighborhood of the equilibrium from an arbitrary initial condition. After reaching the neighborhood, the temporal evolution is then governed by a linear dynamical system (which can be solved in closed form). An analogous result holds for hyperbolic periodic orbits in vector fields [21]. To turn periodic orbits in vector fields into fixed points in maps, one can use Poincaré return maps, which faithfully capture properties of periodic orbits. A Poincaré map can be interpreted as a discrete dynamical system, so the problem of determining how long it takes to reach a hyperbolic stable periodic orbit from arbitrary initial conditions in a vector field is reduced to the problem of determining how long it takes to reach the neighborhood of a hyperbolic fixed point in a discrete dynamical system.
Our work considers how long it takes to reach a periodic orbit of a differential equation-starting from an arbitrary point in phase space-by using a Poincaré return map of its associated vector field. For simplicity, suppose that a return map (which is built from a Poincaré section) is unimodal. If we approximate the unimodal Poincaré map by using a unimodal function f (x), then we can use f (x) in our algorithm to estimate the convergence time to the periodic orbit. Periodic motion is ubiquitous in models (and in nature), and it is important to explore how long it takes to converge to such behavior.
In this paper, we prove a theorem for the rate of convergence to stable periodic orbits in discrete dynamical systems. Our basic strategy is as follows. We define the neighborhood I p of a hyperbolic fixed point, and we calculate what fraction v of the entire phase space I is mapped into I p after q iterations. Using m(w) and m(I), respectively, to denote the measures of w and I, a point that is selected uniformly at random from I has a probability of m(w)=m(I) to reach I p in q iterations. To illustrate our ideas, we will work with a one-dimensional (1D) discrete dynamical system x n z1~f (x n ; r) that is governed by a unimodal function f and is parametrized by a real number r. We focus on unimodal functions for two primary reasons: (i) many important results in dynamical systems are based on such functions; and (ii) it is simpler to illustrate the salient ideas using them than with more complicated functions.
To determine the set that is mapped into I p , we take advantage of the fact that points in periodic orbits are repeated periodically, so their corresponding neighborhoods must also repeat periodically. In theory, an alternative procedure would be to iterate backwards from I p , but this does not work because one cannot control successive iterations of f {1 . The function f is unimodal, so it is not bijective and in general one obtains multiple sets for each backward iteration of a single set. The number of sets grows geometrically, and one cannot in general locate them because an analytical expression for f {1 is not usually available.
To explain the main ideas of this paper and for the sake of simplicity, consider a stable periodic orbit O p of period p that is born in p saddle-node bifurcations of f p . Every point and I i all repeat periodically. Roughly speaking, we will build the interval I p from the interval I i .
Consider a plot in which points along the horizontal axis are mapped via f to points along the vertical axis (as is usual for 1D maps). The orbit O p is periodic with period p, so x j [O p implies that f(x j ,f q (x j )), q~0,1, . . .g yields p periodic points with a horizontal axis location of x j . We say that these points are located in the ''column" x j . Because f q (x j )~x i [O p for some q, we obtain p points located in the same column x j . These points are given by f(x j ,f q (x j ))~(x j ,x i ), i~1, . . . ,pg. As we have indicated above, each point (x j ,x i ) is associated with an interval I i . No matter how many iterations we do, the fact that the orbit is periodic guarantees that there are exactly p intervals in the same place (where the points (x j ,x i ) are located). We thereby know the exact number and locations of all intervals.
To complete the picture, we must also take into account that if there exists an interval W qij such that f q (W qij )~I i , then any point of W qij will reach a point of I i in at most q iterations. The geometric construction above yields the interval W qij , as one can see by drawing a pair of parallel line segments that intersect both f q and the endpoints of the interval I i . We will approximate f q by a set of such line segments so that we can easily calculate the intersection points.
The remainder of this paper is organized as follows. First, we give definitions and their motivation. We then prove theorems that indicate how long it takes to reach the interval I i from an arbitrary initial condition. We then construct algorithms to implement the results of the theorems. Finally, we discuss a numerical example and then conclude.

Definitions
Consider the discrete dynamical system x n z1~f (x n ; r), f : I?I , where f (x; r) is a one-parameter family of unimodal functions with negative Schwarzian derivative and a critical point at x~C. Without loss of generality, we suppose that there is a (both local and global) maximum at C. At a critical point of a map f , either f 0~0 (as in the logistic map) or f 0 does not exist (as in the tent map). Some of the results of this paper related with critical points only require continuous functions, which is a much weaker condition than the requirement of a negative Schwarzian derivative.
Remark 1 Because f has a negative Schwarzian derivative, f q does as well (because it is a composition of functions with negative Schwarzian derivatives). By using the chain rule, we obtain f q 0~0 only at extrema. Therefore, f q 0 =0 between consecutive extrema. The Minimum Principle [22] for a function with negative Schwarzian derivative then guarantees that there is only one point of inflection between two consecutive extrema of f q . If there were more than one point of inflection, then f q 00~0 at least two points. One of them would be a maximum of f q 0 , and the other one would be a minimum. This contradicts the Minimum Principle. Consequently, the graph of f q between two consecutive extrema has a sigmoidal shape (i.e., it looks like or ), which becomes increasingly steep as q becomes larger. This fact makes it possible to approximate f q between two consecutive extrema by a line segment near the only point of inflection that is located between two consecutive extrema.
Because the Schwarzian derivative of f is negative, Singer's Theorem [23] ensures that the system (1) has no more than one stable orbit for every fixed value of the parameter r. Additionally, the system (1) exhibits the well-known Feigenbaum cascade [24][25][26], which we show in the bifurcation diagram in Fig. 1.
For a particular value of the parameter r, the map f p has p simultaneous saddle-node (SN) bifurcations, which result in an SN p-periodic orbit. As r is varied, the SN orbit bifurcates into a stable orbit fS i g p i~1 and an unstable orbit fU i g p i~1 . The points S i and U i are, respectively, the node and the saddle generated in an SN bifurcation, so U i is the nearest unstable point to S i (see Fig. 2). In other words, the points in the stable orbits (called "node orbits") are node points, whereas the points in the unstable orbit (called "saddle orbits") are saddle points. From Remark 1, we know that the neighborhoods of these points are concave or convex.
The derivative of f p is 1 at the fixed point where the SN bifurcation takes place. As one varies r, the derivative evaluated at that bifurcation point changes continuously from 1 to {1. When the derivative is {1, the stable orbit (i.e., the node orbit) undergoes a period-doubling bifurcation. As a result, the stable orbit becomes unstable (yielding the orbit fU i g p i~1 ) and two new stable orbits (fS i1 g p i1~1 and fS i2 g p i2~1 ) appear. The points S i1 and S i2 are nodes, and the point U i is a saddle. From our geometric approach, the intervals (S i1 ,U i ) and (S i2 ,U i ) that are generated via the period-doubling bifurcation behave in the same way as the interval (S i ,U i ) that was generated in the SN bifurcation. Therefore, we can drop the indices "1" and "2" and write (S i ,U i ) for the orbits that arise from both the SN bifurcation and the period-doubling bifurcation. i that satisfy Because we can control the evolution inside I Pi , we can examine how long it takes to reach I Pi starting from an arbitrary point x[I. As we will see below, to obtain this result, we need to discern which subintervals of I are mapped by f q into I Pi for arbitrary q. The first step in this goal is to split the interval I in which f q is defined into subintervals in which f q is monotonic.
be the set of points at which f q has extrema. Let B~fq 1~a ,q k~b g, and we recall that we are considering the interval I~½a,b. We will call P mon{r~A |B~fq 1 ,q 2 ,q 3 , . . . ,q k{1 , q k jq 1 vq 2 v . . . vq k{1 vq k g the partition of monotonicity of f q . We will call I qj~½ q j ,q j z1 (where j~1, . . . ,k{1) the jth interval of monotonicity of f q . By construction, I~½a,b~[ j I qj , and f q is monotonic in I qj . As we will explain below, one can calculate intervals of monotonicity I qj easily by using Lemmas 1 and 2.    Once we know the intervals in which f q is monotonic, it is easy to obtain subintervals of I that are mapped by f q into I Pi .
We proceed geometrically (see Figs. 2, 4, and 5): N draw parallel lines through the points x 0 i and x i (i.e., through the endpoints of I Pi ); N obtain the points at which the lines intersect f q ; N calculate which points are mapped by f q into the intersection points of (ii), for which one uses the fact that f q is monotonic in I qj~½ q j ,q j z1; N determine, using the points obtained in (iii), the interval that is mapped by f q into I Pi .
Using this geometric perspective, we make the following definitions.

Definition 3 Let
where the points x i,qj ,x Remark: If we did not take point (ii) into account, then f q would not be monotonic in W qij .
By construction, all points x[W qij reach I Pi in at most q iterations (see Figs. 2, 4, and 5). That is, f l (W qij )~I Pi for lƒq.

Theorems
Once we know W R , we can calculate the probability that a point picked uniformly at random from phase space is located in W R . We can then calculate the probability that that point reaches a capture interval of O p in at most q iterations. We let m(W R ) denote the measure of W R , and we have the following theorem.
Theorem 1 Let O p~f S i g p i~1 be a stable p-periodic orbit of the system (1). Given an arbitrary point x[I, the probability to reach a capture interval of O p after at most q iterations is i , respectively. We thereby construct the subinterval W qij . We depict the mapping of the subinterval W qij using a filled green arrow. doi:10.1371/journal.pone.0092652.g005 Proof 1 From the definition (4) of W R , all x[W R satisfy f l (x)[I p for lƒq. There always exist values of lvq such that f l (x)[I p because extrema of f l (x) that satisfy lvq are necessarily also extrema of f q , and points belonging to the latter set of extrema reach a capture interval of O p after at most l iterations (see Lemma 1 below). Consequently, one reaches I p from x[W R after at most q iterations (and we note that it need not be exactly q iterations). Thus, the probability to reach I p from an arbitrary point x[I after at most q iterations (i.e., the probability that x[W R ) is Corollary 1 With the hypotheses of Theorem 1, the probability to reach a capture interval of O p in exactly q iterations is This answers the question of how long it takes to reach a capture interval of a p-periodic orbit from an arbitrary point. However, we also need to calculate m(W R ). To do this, we need to understand the structure of W qij . As the following lemma indicates, some of these subintervals are located where f q is monotonic and others contain extrema of f q .
Lemma 1 If f : I?I is an unimodal C 0 function with a critical point at C, then f q (x) has extrema (i) at points for which f q{1 (x)~C; (ii) at the same points at which f q{1 (x) has extrema. Proof 2 (i)For all x[I such that f q{1 (x)~C, we know that f q (x)~f (f q{1 (x))~f (C). Therefore, f q has an extremum because f has an extremum.
(ii)Write I~J L |fCg|J R , where J L~( a,C) and J R~( C,b), so f is a monotonic function on the intervals J L and J R .
(ii.a) If x[J L or x[J R and the function f q{1 (x) has an extremum, then we know that f q{1 (x) is a monotonically increasing function on one side of x and a monotonically decreasing function on the other. Consequently, f q (x)~f (f q{1 (x)) is the composition of two monotonic functions (f and f q{1 ), both of which are increasing (or decreasing) on one side of x. On the other side of x, one of them is increasing and the other is decreasing. Therefore, there is an extremum at x.
(ii.b) Otherwise, if f q{1 (x)~C, then we see straightforwardly that f q has an extremum.
We have just seen how to determine the locations of extrema of f q . We also need to know the values that f q takes at these extrema.
As we will see below, if the system (1) has a stable p-periodic orbit and qwp, then the values that f q takes at its extrema are the same as those that f p takes at its extrema. This makes it possible to calculate the subintervals W qij that are associated with extrema of f q by using I PiC and the derivative of f . Lemma 2 Let O p be a stable p-periodic orbit of the system (1). The coordinates of the extrema of f q (where qwp) are (x iC ,f q{i j p (C)), where x iC denotes the points x[I such that f i (x)~C, the index i takes values of i~0,1, . . . ,q{1 (where we note that f 0 :Id is the identity map), and q{ij p~( q{i)mod p.
Proof 3 According to Lemma 1, the extrema of f q are (i) x[I such that f q{1 (x)~C; (ii) x[I such that f q{1 (x) is an extremum. It thus follows that the extrema of f q{1 are (ia) x[I such that f q{2 (x)~C; (iib) x[I such that f q{2 (x) is an extremum. Repeating the process, we obtain that extrema of f q are located at x iC , where i~0,1, . . . ,q{1. The value of f q at x iC is Because O p is a stable p-periodic orbit, there exists one point of O p near C that is repated periodically after p iterations. Consequently, fC,f (C),f 2 (C), . . . ,f p{1 (C)g is a periodic sequence and f q{i (C)~f q{i j p (C).

Algorithms
As we discussed above, Lemmas 1 and 2 determine intervals of monotonicity (see definition 2), and they also make it possible to construct algorithms for calculating W qij .
For these algorithms, we approximate f q by line segments in the subintervals in which f q is monotonic. This approximation is very good unless one is extremely close to an extremum (see Fig. 6), and this is already the case even for relatively small q (as we will demonstrate below). Additionally, recall that W qij is determined by the intersection points of f q with line segments. Therefore, once we have approximated f q by a set of line segments, it is straightforward to calculate those intersection points.
Algorithm 1 (Calculating coordinates for extrema of f q ) Suppose that we know the coordinates of the extrema of f q{1 . According to Lemma 1, the extrema of f q are located at the points (i) x[I such that f q{1 (x)~C and f q{1 (x) is not an extremum; (ii) x[I such that f q{1 (x) is an extremum. We know the extrema in (ii) by hypothesis. To find the extrema in (i), we need to calculate the points x[I that satisfy f q{1 (x)~C. Because we know the coordinates of extrema of f q{1 , we construct the lines that connect two consecutive extrema (see Fig. 7). Let x nz1~a x n zb be the equation for such a line. We solve ax n zb~C to obtain a seed that we can use in any of the many numerous numerical methods for obtaining roots of nonlinear algebraic equations. Observe that f q is monotonic in the interval in which the line x nz1~a x n zb is defined. This circumvents any problem that there might otherwise be in obtaining a good seed to ensure convergence of the root solver. Moreover, we have as many seeds as there are points x[I that satisfy f q{1 (x)~C. Note that we need to construct both the line that connects (a,f (a)) with the first extremum of f q{1 and the line that connects (b,f (b)) with the last extremum of f q{1 .
To calculate the points x[I for which f q{1 (x) is an extremum, we apply this algorithm recursively, and we note that we know by hypothesis that f has an extrememum at C. We first build the line segments that connect (a,f (a)) with (C,f (C)) and (C,f (C)) with (b,f (b)). These two line segments give seeds from which to determine the points x[I that satisfy f (x)~C. We thereby obtain the coordinates for the extrema of f 2 . We then use the same procedure to obtain the coordinates for extrema of f 3 , f 4 , . . . , f q .
We will see below that if the system (1) has a stable p-periodic orbit and q&p, then the points x[I with f q{1 (x)~C are given to a very good approximation by the intersection points of two lines. Moreover, as one can see in Fig. 6, the value q~6 is already large enough to approximate f q very successfully by a set of line segments when f is the logistic map.
Algorithm 2 (Calculation of W qij in the system (1)) Suppose that we know the coordinates of the extrema of f q (e.g., by computing them using Algorithm 1). We want to obtain W qij from the definition (3), where and the ith capture interval is To determine the points x i,qj and x 0 i,qj , we first approximated them by replacing f q by line segments that connect consecutive extrema of f q (i.e., by the same procedure that we use in Algorithm 1 to obtain approximations of points). Using the approximations of x i,qj and x 0 i,qj , we construct the interval I app~( x i,qj ,x 0 i,qj ) and then check if there is an extremum of f q in I app . (This is trivial because we know the coordinates of the extrema of f q .) We need to consider two cases.
(i)The map f q has no extrema in I app . This is equivalent to case (i) of Algorithm 1. We use the approximations of x i,qj and x 0 i,qj as seeds in a numerical root-finding method.
(ii)The map f q has extrema in I app . This is equivalent to case (ii) of Algorithm 1.
If there is an extremum of f q in I app , then that extremum is necessarily one of the extrema given by Lemma 2: (x iC ,f q{i j p (C)). Because f i (x iC )~C and f is a continuous function, there must exist an interval I iC such that x iC [I iC and f i (I iC )5I Pi,C Taking into account that x iC is known, we construct the sequence Let L i,k be the linear map whose graph is the line of slope f 0 (x i,k ) that intersects the point x i,k . If the period p of the orbit is sufficiently large, then we can approximate f near x i,k (where k~0,1, . . . ,i{1) by the linear map L i,k . Thus, instead of iterating I iC with the map f to obtain I Pi,C , we iterate I iC with the linear map L i,k that approximates f . That is, Because each L i,k is a linear map, it is straightforward to compute L {1 i,k and hence to compute At the end of this section, we will discuss the error that is introduced by this approximation. Figure 6. Outside of the intervals W qi,j , we approximate the map f q using line segments. A line segments connects the upper endpoint of the interval W qi,j to the lower endpoint of W qi,jz1 . The map f 6 is very well approximated using line segments as long as one is not too close to an extremum. We again use the logistic map to illustrate our procedure. The blue curve is a period-6 supercycle and r&3:99758311825456726610. See The interval I iC that we have just constructed is the interval that we seek.
In Algorithm 1, we constructed line segments that connect two consecutive extrema of f q . They are located in the intervals ½q j ,q jz1 ) and ½q jz1 ,q jz2 ), respectively. We now have intervals W qi,j 5½q j ,q jz1 ) and W qi,j z1 5½q jz1 ,q jz2 ) that contain these two consecutive extrema of f q , so we construct the line segment that connects the upper endpoint of W qi,j to the lower endpoint of W qi,jz1 . (Note that we do not connect the two extrema directly via a line segment.) For q&1, this line segment approximates f q outside of the intervals W qi,j and W qi,jz1 . See Fig. 6, which illustrates (for the case when f is the logistic map) that we can approximate f 6 by a set of line segments for q~6. We can then use these line segments in Algorithm 1, and we do not need numerical computations to find the intersection points.
As we discussed previously, we can replace f q by linear expressions to approximate the intersection points when determining W R in Algorithms 1 and 2. Replacing f q by a linear approximation simplifies operations and reduces the amount of calculation. To determine the desired intersection points, we have thereby replaced a numerical method for obtaining roots of nonlinear algebraic equations by an analytical calculation that uses a system of two linear equations. We now estimate the error of replacing f q by lines segments. The line segments that replace the function f q intersect f q very close to the unique point of inflection between a pair of consecutive extrema of f q (see Remark 1 and Consequently, the error of approximating f q by the line f q (x inf )zf q 0 (x inf )(x{x inf ) is where we have taken into account that there are more than 2 q local extrema of f q in the interval ½a,b. The exponential growth of 2 q enforces a fast decay in the error. Consequently, using line segments to approximate f q is an effective procedure with only a small error.

Numerical Example
Algorithms 1 and 2 are based on the same procedure: approximate f q (x) by a line y(x)~axzb and solve y(x)~C to obtain an approximation of the f q (x)~C (instead of solving f q (x)~C directly). In this section, we consider an example application of Algorithm 1.
To obtain the critical points of f q z1, we need to calculate the points that satisfy f q~C . Suppose that q~6 (and again see Fig. 6 for an illustration of the line-segment approximation with q~6 for the logistic map). The biggest distance between consecutive Figure 7. Graphs of f 6 (blue) and f 10 (red) for the same value of the parameter r (when f is the logistic map) as in Fig. 6. The dark pink line joins two consecutive extrema of f 6 , and the black line is the tangent line that crosses through the inflection point. Both lines are approximations to f q . As expected, the approximation is better for the larger value of q. doi:10.1371/journal.pone.0092652.g007 extrema occurs near the critical point C, so we approximate f 6 by a line segment in this region to obtain an upper bound for the error. The extrema are located at (4:525|10 {1 ,2:414|10 {3 ) and (4:787|10 {1 ,9:994|10 {1 ), and they are connected by the line y&38:053 x{16, from which we obtain the approximation x app &0:453 for the solution of f 6 (x)~C. From direct computation, the value of x that satisfies f 6 (x)~C is x&0:465. The relative error is E rel &2:58%, and this is the largest error in this example from all of the approximating lines segments. As we showed in equation (10), the error decreases exponentially. Hence, when we approximate f 6zm using line segments, the relative error will be bounded above by E rel &2:58=2 m %. One can observe this decrease in error in Fig. 7, in which we plot both f 6 and f 10 for the logistic map and the same parameter value r. Observe that several extrema of f 10 lie between onsecutive extrema of f 6 , so using the line-segment approximation in f 10 induces a much smaller error than using it in f 6 .

Conclusions and Discussion
When studying dynamical systems, it is important to consider not only whether one converges to periodic orbits but also how long it takes to do so. We show how to do this explicitly in onedimensional discrete dynamical systems governed by unimodal functions. We obtain theoretical results on this convergence and develop practical algorithms to exploit them. These algorithms are both fast and simple, as they are linear procedures. One can also apply our results to multimodal one-dimensional maps by separately examining regions of parameter space near each local extremum.
Although we have focused on periodic dynamics, the ideas that we have illustrated in this paper can also be helpful for trying to understand the dynamics of chaotic systems. Two important properties of a chaotic attractor are that (i) its skeleton can be constructed (via a ''cycle expansion") by considering a set of infinitely many unstable periodic orbits; and (ii) small neighborhoods of the unstable orbits that constitute the skeleton are visited ergodically by dynamics that traverse the attractor [4]. In Refs. [21,22], Schmelcher and Diakonos developed a method to detect unstable periodic orbits of chaotic dynamical systems. They transformed the unstable periodic orbits into stable ones by using a universal set of linear transformations. One could use the results of the present paper after applying such transformations. Moreover, the smallest-period unstable periodic orbits tend to be the most important orbits for an attractor's skeleton [4], so our results should provide a practical tool that can be used to help gain insights on chaotic dynamics.
Once unstable orbits has been transformed into stable ones we can use results of this paper to answer the above question.