T. M. Jonassen

## A Mathematica Implementation of Nonlinear Dynamical Systems Theory via the Spider Algorithm and Finding Critical Zeros of High-Degree Polynomials NB   CDF   PDF

Important properties pertaining to families of discrete dynamical systems are furnished here by studying the kneading theory developed by Milnor and Thurston, and subsequently implementing the spider algorithm, developed by Hubbard and Schleicher. The focus is on identifying crucial combinatorial and numerical properties of periodic critical orbits in one-dimensional discrete dynamical systems, which are generated by iterating real quadratic polynomial maps that constitute an important class of unimodal systems.

### Introduction

The mathematical concepts and notations, which facilitate the introduction of the kneading theory developed by Milnor and Thurston [1] and the spider algorithm developed by Hubbard and Schleicher [2], are summarized first.

A mapping of an interval given by gives rise to a dynamical system on by iteration. An orbit of a point is a sequence of points given by We use the standard notation , where indicates the -times composition of with itself. An orbit is called periodic with a primitive period if for some , but for . A critical periodic orbit is a periodic orbit that contains the critical point of . A critical periodic orbit is also called a superstable periodic orbit. A point is called critical if the derivative of the map is zero at the point.

Symbolic dynamics are techniques used in the study of dynamical systems. The simplest form converts orbits into sequences of symbols from an alphabet . In our unimodal case, the alphabet consists of three symbols, . The iteration process (the dynamical system) translates into a simple operator on the symbol space consisting of infinite words from . In this article we look at some elementary symbolic dynamics for simple one-dimensional systems , where is a smooth function of the interval such that has only one critical point in . We may write such that is decreasing on and increasing on . Let and . The address of a point , denoted by , is given by if , if , and if . Since by assumption, the orbit of a point is contained in , and we may assign an infinite sequence of symbols to the orbit according to the rule The sequence is also called the itinerary of the point There is a natural mapping on the space of symbol sequences compatible with the dynamics on , the shift map , defined by Let ; then clearly . Hence the action corresponds to shifting the sequence of symbols one place to the left and forgetting the first symbol in the sequence . The critical point plays a special role in the dynamics of The symbol sequence , which gives the dynamics of the critical orbit, is called the kneading sequence of and is denoted by The kneading sequence of for a superstable periodic orbit is periodic. A periodic kneading sequence is written with an overbar, as in the sequence of period 7. There is a periodic orbit for this sequence, while there is no orbit for the dynamical system corresponding to the sequence . One of the main problems in this field is to decide if a given periodic symbol sequence corresponds to an orbit in the dynamical system A symbol sequence with a corresponding orbit in the dynamical system is called admissible. We apply the theory developed in [1] and [3] to obtain an algorithm that decides if a given symbol sequence is admissible. This algorithm is based on the fact that the kneading sequence is minimal with respect to the lexicographic order denoted by on the symbol space . In particular, if is the shift map on the symbol space, then , , for any admissible kneading sequence of length

The spider algorithm [2, 4] was designed to study certain properties of the Mandelbrot set for families of dynamical systems . We apply the spider algorithm to real unimodal systems generated by quadratic polynomials using admissible kneading sequences. The construction of a periodic kneading sequence of length corresponds to finding a certain real solution of a polynomial equation of degree . This leads to a numerical method for studying the parameter space of unimodal systems generated by quadratic polynomials. The Sharkovsky theorem [5, 6] is illustrated in [7] as a special case, “period 3 chaos”.

It is easy to see that it is sufficient to study a single representative among the nondegenerate quadratic polynomials: let , where . Then for any such map there is a homeomorphism (in fact, of the simple form ) such that , where the quantities , , and depend on , , and . Hence and are topologically conjugate and have the same dynamics. In the following, we use as the representative for the quadratic polynomials.

### Some Simple Properties of the Family

A fixed point is called hyperbolic if . We may replace by any iterate of and hence an -periodic orbit is called hyperbolic if A hyperbolic invariant set is called hyperbolic (in the one-dimensional case) if the derivative of is greater (smaller) than 1 in absolute value for all points in the invariant set. In higher dimensions, with , where is a compact smooth manifold, hyperbolicity means that the tangent bundle of the invariant set has a -invariant splitting into a contracting subbundle and an expanding subbundle. It might of course happen that one of these is empty (the invariant set is an expeller or attractor). For a more precise definition of hyperbolicity, see [8], and for an application showing how hyperbolicity ensures stability of the topological type of the system under small perturbations, see [9].

We now state some simple properties of the dynamical system All of the properties listed are easily proved using elementary calculus, so we omit the calculations. The statements on symbolic dynamics and hyperbolicity can easily be proved using the techniques in [8] or a modification of the arguments in the next section. In the following, let , , and whenever these quantities are real (i.e., when ). Clearly, is the only critical point of .

1. If , then for all .
2. If , then has two fixed points given by and .
3. If , then the interval is invariant under ; that is, .
4. If , then every periodic orbit of is repelling; has periodic orbits of every primitive order, but none of them contains the critical point.
5. If , then has an invariant Cantor set such that the restriction is topologically conjugate to a one-sided shift on two symbols. Furthermore, the set is hyperbolic. For any point , we have .
6. The mapping is topologically conjugate to the mapping via the homeomorphism , , with , , and , where
.

### Real Spiders and the Spider Map

Consider the -periodic orbit containing the critical point under the map for a suitable choice
of :

Since , we have , and the correct point to choose in the fiber is given by the kneading symbol at that location in the periodic orbit. In our case, the fiber is empty if , contains exactly one point if , and contains two points if . The coding, , of an orbit under is done according to the rule

The kneading sequence of is the symbolic orbit of the critical point, . In our setting the kneading sequence is periodic. We use the notation to denote that the finite symbol sequence under the bar is repeated an infinite number of times. It is easily seen that not all symbol sequences are compatible with the underlying dynamical system. In fact, it can be shown that there is at most one order of points that is compatible with a kneading sequence.

A real spider is a very special case of the spiders defined on the Riemann sphere for complex systems. On the Riemann sphere, a spider is an equivalence class of curve systems connected in ∞, the “body” of the spider, and the curves going out from this point may be thought of as the “legs”. The legs are used to impose an ordering of the points in . However, in there is a natural ordering, so the space of real spiders associated with the dynamical system takes the form of -tuples of real numbers subject to a set of inequalities , where .

#### The Spider Space

Let and be an index vector, where with if Let be the subset The space equipped with the natural inherited topology from is called the real spider space associated with . A mapping is called a spider mapping. We will later index the space by a periodic admissible kneading sequence, writing

Example 1. Consider the real spider space and let The map is clearly well defined on as ; the first component in the image is negative and the second component is positive. Suppose has a fixed point Then , , and . We find , , and By rearranging these equations, we have , , and This corresponds exactly to the orbit , and hence any such fixed point corresponds to a superstable period-3 orbit under .

Example 2. Consider the real spider space and let be the map defined by

We check that this map is well defined; that is, we show that . First, and means and . Next, because and , we have that ; hence, . Therefore the map is well defined. Assume as in Example 1 that has a fixed point , that is, . The fixed point equation gives us that , , , , and . Substituting, we rewrite these equations as , , , , and . This is exactly a critical period-5 orbit for our polynomial family , where
Hence, a fixed point for this spider map corresponds to a superstable period-5 orbit with kneading sequence .

Example 3. In Examples 1 and 2 we used spider maps with fixed points that correspond to periodic critical orbits obeying certain combinatorial properties. In the two first cases we used kneading sequences that are compatible with the dynamics of the system . We now choose a map that corresponds to an incompatible kneading sequence (later to be called an inadmissible kneading sequence). Consider the sequence . This gives us the real spider space and suggests that we define the spider map as

We show that is not well defined. Note that, since and , we have so . Hence, . However, so , that is, , implying that is not closed under . As a consequence, some of the roots of the component functions may become complex after a few iterations of . Indeed, this is exactly what happens as we show later on.

### Finding Admissible Kneading Sequences Using the Minimality of the Kneading Sequence with Respect to the Lexicographic Order

We apply the theory developed in [1] to obtain an algorithm that decides if a kneading sequence is admissible. This algorithm is based on the fact that the kneading sequence is minimal with respect to the lexicographic order denoted by as defined in the next subsection. In particular, if is the shift map on the symbol space, then , , for any admissible kneading sequence of length .

#### The Lexicographic Order

Let be a three-letter alphabet with the ordering , and let be the set of infinite words from subject to the restriction: suppose and are two words in containing the letter , say and , where and do not contain the letter ; then .

Let , where we write . Assume that for . We define , where if and if . In other words, as is the sequence of addresses coming from the dynamical orbit , the quantity determines the orientation properties for at the point . Note that is decreasing (orientation reversing) for (corresponding to the symbol ) and increasing (orientation preserving) for (corresponding to the symbol ).

We can now define a signed lexicographic ordering, denoted by < (less) and (less or equal), for two elements . Assume that for , then if either and , or and . We write if or .

The following lemma is used to construct the general algorithm for finding admissible kneading sequences.

Lemma 1. Let be a kneading sequence of a unimodal map let be the itinerary of a point , and let denote the shift map on the symbol space. Then for all and In particular, .
Proof. See [3] or [1].

We may use the special case of Lemma 1 to decide if a given candidate for a periodic kneading sequence of length is admissible. We simply need to test if for .

##### Mathematica Programs

We first generate possible candidates for admissible periodic kneading sequences and then reduce the number of candidates by excluding some sequences that cannot be candidates. Clearly all sequences must start with if , and it is easy to prove that sequences of the form and , where W is a word from including the empty word, cannot be admissible for . Here are the programs that generate the candidates.

These programs find the signed lexicographic order for the candidate strings.

These programs test the candidates and return the ones that are admissible sequences.

Here is what we wanted to obtain. Due to the large number of elements in the output, we only count the number of words in each list.

### The Ordering of Points in an Admissible Periodic Kneading Sequence (Or, Forming the Real Spiders)

Let be an admissible kneading sequence with and a word from the alphabet with . Let the corresponding dynamical orbit be , where in our case. The problem is to order the points in the orbit on the real line , in other words, to find a bijective function based on the kneading information. Let and denote the number of and in the word , respectively. Clearly , and if We have some trivial information about the function . Clearly , , and . We compute the image of using the same trick used for computing the admissible kneading sequences, that is, we sort points according to their lexicographic order. The following lemma from [1] relates the order of points in the dynamic space with the order of addresses in the symbol space.

Lemma 2. Let be a unimodal map where the critical point is a minimum, let with , and let and denote their itineraries. Then with respect to the signed lexicographical order.
Proof. See [1], Section 3.

We apply Lemma 2 to find the ordering of points in the dynamical space of an admissible kneading sequence given by , where is a word from the alphabet as follows. The symbols in the periodic word are assigned to symbolic points in the dynamic space, that is, to indices , and these are split into three groups according to their symbol in the kneading sequence. For example, the sequence is mapped to . The problem is then reduced to sorting the first and third groups according to their relative positions in the dynamic space. Now we just compare two versions of the symbol sequence by rotating left the correct number of times according to the symbol position in the string, so this symbol becomes the first symbol, given by the indices we have already found. We then apply Lemma 2 to determine their relative positions in the dynamic space.

##### Mathematica Programs

The following program gives a version of the map operating on words that correspond to dynamical orbits.

Here is an example using the admissible sequence .

We use the function JSortMap to produce a function for generating a suitable element of the spider space associated with a given admissible kneading sequence . The element is then used as an initial point for the spider algorithm given in the next section to generate the dynamical orbit for the system . We choose this spider because it has equally spaced points in each of the intervals and according to the numbers of and in the word . The programming here is straightforward; we only need to find an “inverse” to the map described by JSortMap.

Here is an example with the admissible sequence .

### Mathematica Implementation of a Spider Map

The simple Examples 1, 2, and 3 suggest how we should define the spider map associated with a periodic kneading sequence. Consider the periodic dynamical sequence

where we have for with . Hence, we have , so , where if the corresponding kneading symbol is , , or .

Implementing a real spider map that chooses correct roots according to a given kneading sequence is easily done in Mathematica. Here, we do not perform any error or sanity checks, so our map simply takes the form

where SpiderRoot returns the correct root according to the symbol in the kneading sequence. In our case, this map can be defined by

### The Spider Algorithm

We now briefly describe the spider algorithm for the system .

We have seen in Examples 1 and 2 that a critical periodic orbit is a fixed point for the spider map constructed by choosing the correct roots according to the combinatorics of the dynamical orbit. Example 1 shows that this fixed point, in the special case of a period-3 orbit, is stable. We might hope that this is true in general for any critical periodic orbit, and hence suggest the following algorithm.

Problem: Find a parameter value such that the real dynamical system has a periodic kneading sequence .

The Real Spider Algorithm:
A)
Choose a finite string of length of symbols, where the first two symbols are and , the next symbols are chosen from the two-letter alphabet , and the last symbol is .
B) Form the map , where the component function is and if the symbol in the string is , , or , respectively.
C) Choose a vector , where the points are ordered according to the dynamics of the periodic orbit for the dynamical system .
D) Form the sequence and stop the iteration process when the sequence of vectors converges to some point (or in ).
E) The parameter with the desired critical orbit is given by .

In Mathematica this is implemented by the following iteration process.

`First[    FixedPoint[        RealSpiderMap[Characters[kneading_sequence]],spider]     ]`

Here kneading_sequence is a string of symbols and spider is an ordered list of real numbers. If the kneading sequence is not compatible with the dynamics, then the returned number is nonreal, that is, a number in .

We define three functions associated with the numerical computation of the spider algorithm: , , and . In all cases is a kneading sequence and is an optional integer passed to FixedPoint or FixedPointList to control the maximum number of iterations. This is necessary in some cases because there is a “bit-flip” on the least significant bit at the fixed point, causing a nonstopping condition in FixedPoint. Note that we check if is an admissible periodic kneading sequence. The function SpiderIterationList returns a list of all steps in the iteration process taken to find the fixed point of the spider map. The function SpiderFixedPoint returns the fixed point (the orbit) associated with the kneading sequence . The function CriticalParameter returns the first component of the fixed point, that is, the parameter for corresponding to the periodic kneading sequence.

### Examples

We now give some examples using the programs. Some variations of the functions are also used; open this cell to see the programs.

#### Dynamics of a Period-15 Orbit

In this example we display the dynamics of the period-15 orbit in two different ways. We first check that this string represents an admissible periodic kneading sequence.

We can represent the superstable orbit graphically in the usual way by drawing the graph of for the value of corresponding to this sequence in red, drawing the diagonal in blue, and then following the orbit through the critical point with green lines. Doing it this way makes it hard to keep track of the orbit properties.

Another way is to simply display how points on the critical orbit move around under the map without taking their exact location into account. We just show their relative locations, which point is mapped to which, and in what direction. In the following diagram the points in the orbit of are shown in black along the horizontal axis and the critical point is marked with a C. The curved lines represent how the points on the orbit are mapped with lines above the horizon meaning a movement from left to right, and those below a movement from right to left.

#### Admissible Sequences of Length 9

In this example we compute all admissible sequences of length 9 and the corresponding critical parameter. Note that this corresponds to finding certain real solutions of a polynomial of degree . There are 28 admissible sequences of length 9. (The function SortedAdmissibleSequences is described in the next subsection.)

We use the function CriticalParameter to find the corresponding values for the dynamical system .

Note that these solutions are all the (real) zeros of the polynomial of degree 256 corresponding to periodic orbits through the zeros of the dynamical system. The polynomial can be computed with Nest.

The following table shows the periodic kneading sequence and the corresponding parameter .

#### Sorted Lexicographical Ordering

We define a function that returns the admissible periodic kneading sequences of length n in sorted lexicographical order. In addition, we have associated the symbols < and with the lexicographical order.

Here is an example with SortedAdmissibleSequences for sequences of length 15. There are 1091 different admissible sequences of this length. The function sorts these with respect to the lexicographical order. The notation means n strings are omitted in the output.

We may use the symbols < (Precedes) and (PrecedesSlantEqual) to test the lexicographic order of two strings. These relations work on any nonempty string from the alphabet . The strings do not need to be of equal length.

#### The Sharkovsky Ordering of

The bifurcation diagram in the following figure shows the attracting set of the critical orbit. We cannot see the repelling periodic orbits in this diagram. The Sharkovsky theorem states the relationship between coexisting periodic orbits without considering stability properties.

The natural numbers are ordered as follows by . With , the Sharkovsky ordering is given by

##### The Sharkovsky Theorem

Let be a continuous map of some interval If has a periodic orbit of primitive period then has periodic orbits of primitive period for all with in the Sharkovsky ordering. In particular, if has a periodic orbit of primitive period three, then has periodic orbits of all periods.

See [5], [6], or [8] for a proof. Because [5] was written in Russian the result was unknown for a long time in the West. A proof of a special case of the Sharkovsky theorem, the theorem named “period-3 implies chaos”, was given in [7] because the authors were unaware of the result in [5]. However, the proof of this special case is much easier than the general proof of the Sharkovsky theorem.

Consider the dynamical system . The Sharkovsky ordering of has the odd numbers (1 excluded) as its greatest numbers in reverse order, , . The Sharkovsky theorem implies that the first period- orbit must come into existence before or at the same time (with respect to the parameter) as a period- orbit when the parameter is varied from to . As shown in our previous examples, there is more than one admissible period- kneading sequence if . Let denote the last occurrence (with respect to the usual order in ) of a superstable period- orbit, an odd number, for . We will find the sequence

using kneading sequences and the lexicographical order. Note that this sequence does not give the bifurcation points in the parameter space, but they are quite good approximations, as the width of the windows containing the attracting periodic windows becomes very narrow on the left side of the bifurcation diagram shown earlier.

Let be an admissible sequence of length , with . It is easily shown that the maximal strings of this length are of the form This means that these critical periodic orbits are the first to appear when moving in the parameter space from the right to the left (see the bifurcation diagram). This fact saves a lot of computation as we do not have to use the function SortedAdmissibleSequences. We may simply generate each sequence of length that is needed.

Equipped with the fact that the maximal odd kneading sequences are of the form , we generate these for a consecutive sequence of odd numbers and compute the corresponding critical parameter using the spider algorithm.

Note that these calculations correspond to finding particular nontrivial real solutions of polynomials of degrees in the range .

We now consider what happens if the spider algorithm is applied to an inadmissible sequence. The sequence of length 5 is not admissible.

In the next computation we apply the spider algorithm with this configuration and we easily see that we should use an initial spider of the form .

We obtain an orbit in . Even if we use a different initial spider, we obtain the same orbit.

This orbit is a critical orbit for the system viewed as a complex dynamical system.

### Conclusion

The main issue in this work was not the implementation of the spider algorithm, which was trivial, but rather the implementation of algorithms that decide if a given string of symbols is compatible with the dynamics of . We used symbolic techniques to show how a dynamical orbit is ordered and how to use symbolic dynamics to obtain numerical results. Mathematica provides excellent tools for this purpose.

### References

 [1] J. Milnor and W. Thurston, “On Iterated Maps of the Interval,” in Dynamical Systems, Proceedings of the Special Year, College Park, MD, Lecture Notes in Mathematics, 1342, Berlin/Heidelberg: Springer, 1988 pp. 465-563. doi .10.1007/BFb0082819. [2] J. H. Hubbard and D. Schleicher, “The Spider Algorithm,” in Complex Dynamics: The Mathematics behind the Mandelbrot and Julia Sets, Proceedings of Symposia in Applied Mathematics, 49, New York: AMS, 1994 pp. 155-180. [3] H-M. Xie, Grammatical Complexity and One-Dimensional Dynamical Systems (Directions in Chaos), Vol. 6, Singapore: World Scientific Publishing Company, Ptc. Ltd., 1996. [4] D. A. Brown, “Using Spider Theory to Explore Parameter Spaces,” Ph.D. thesis, Ithaca, NY: Cornell University, 2001. [5] A. N. Sharkovsky, “Coexistence of Cycles of a Continuous Transformation of a Line into Itself,” Ukrainshii Mathematicherhii Zhurna, 16(1), 1964 pp. 61-71 (in Russian). [6] A. N. Sharkovsky, Y. L. Maistrenko, and E. Yu. Romanenko, Difference Equations and Their Applications (Mathematics and Its Applications, Vol. 250), New York: Kluwer Academic Press, 1993. [7] T. Y. Li and J. A. Yorke, “Period Three Implies Chaos,” The American Mathematical Monthly, 82(10), 1975 pp. 985-992 www.its.caltech.edu/~matilde/LiYorke.pdf. [8] R. Devaney, An Introduction to Chaotic Dynamical Systems, Menlo Park, CA: Benjamin/Cummings Publishing Company, 1986. [9] T. M. Jonassen, “Lifts of One-Dimensional Systems: I. Hyperbolic Behaviour,” Journal of Physics A: Mathematical and General, 30(3), 1997 pp. 937-948. T. M. Jonassen, “A Mathematica Implementation of Nonlinear Dynamical Systems Theory via the Spider Algorithm and Finding Critical Zeros of High-Degree Polynomials,” The Mathematica Journal, 2010. dx.doi.org/doi:10.3888/tmj.11.3-2.