Dynamical Systems and Chaos

Lab 1 --- Dynamics of One-Dimensional Maps

In this lab you will use Mathematica to compute and plot orbits of one-dimensional dynamical systems.  The goal will be to illustrate some of the features described in class, such as attracting fixed points, periodic points and chaos.  To make things easier, you will use a package of Mathematica programs that I have written which are contained in the Initialization section of this notebook.  To get Mathematica started and to enter all these programs into memory, go to the Kernel menu and inside the Evaluation submenu you should select the item called Evaluate Initialization.  It is not necessary to look in detail at these programs now.  You will see how to use them as we go along.

This lab is organized into several sections whose titles are visible below.  To begin the first section just double-click on the bracket to the far right of the title Dynamical systems and orbits.  The bracket will expand to reveal lots of Mathematica cells containing explanatory text and commands.  Just work your way through the section reading the explanations and executing the commands.  When you finish a section you can close it up again by double clicking on the expanded bracket at the extreme right.  Then go on to the next section.  Enjoy !

Dynamical systems and orbits

To use the dynamical systems package, you need to enter a dynamical system.  The cell below defines the logistic function f1[x] and then uses it to define a "dynamical system object" called f.  This extra step lets the system know that the package I wrote should be applied to compute iterates, orbits, etc.  To execute a Mathematica command, click the mouse somewhere in the "cell" containing the command and do a Shift-Return, that is, hold down the Shift key and then hit the Return key.

f1[x_]  :=  a*x*(1 - x)f  =  DS[f1]

To compute, say, the 12-th iterate of f applied to the initial state x0 = 0.1, we can use the standard exponential notation below.  First we will set the value of a (the semicolon prevents the a value from being printed).

\(a\  = \ 2.9; 
f^12[0.1]\)

You can edit previous commands and reexecute them.  Go back to the last cell and compute some larger iterate of the same point, say the 1200-th one or the 12000-th one.

To compute an entire orbit, instead of just one point on it, just enclose the exponent in list brackets as shown below.

f^{100}[0.1]

The return value is an "orbit object", which is just a list of points contained inside an Orbit[ ] symbol.  The points represent the orbit x0,x1,.... Once you have computed an orbit object you can use a simple bracket notation to print out entries of it.  But you need to give the orbit a name first.

orb  =  f^{100}[0.1];

Now you can use, for example, orb[27] to read off x27.  You can also put a range of indices inside to get a sublist or use an empty bracket to print the orbit as a plain list.

orb[27]orb[{0, 3}]orb[]

It seems that this orbit is attracted to a fixed point of f located at  p=0.65517. This is even clearer when we plot the orbit.  There are many built-in plotting functions in Mathematica and I have included several more in the package.  For example, the cell below produces a time-series plot.

TimeSeriesPlot[orb, PlotRange -> All];

You can use the mouse to resize graphic images for clearer viewing.

Another interesting plot is shown below.

CobwebPlot[f, orb, {0, 1}];

Periodic Attractors

It is easy to use simple iteration and plotting to search for attracting periodic orbits.  In the cell below, we plot the time series for an orbit of the logistic map with a=2.5.

\(a\  = \ 2.5; 
orb\  = \ f^{100}[0.1]; 

    TimeSeriesPlot[orb, PlotRange -> All]; \)

The commands in the previous cell can be used to explore the changing dynamics of the logistic map as the parameter  a is changed.  Go back and change the value of  a and see what happens.  Try the values a=2.9,3.1,3.4,3.5,3.55,3.83,3.845.In each case you should be able to identify the period of the attractor.

Cobweb plots are also interesting.  Try several   a values in the command below.

\(a\  = \ 2.5; 
orb\  = \ f^{100}[0.1]; 
CobwebPlot[f, orb, {0, 1}]; 
    \)

For high period attractors it is sometimes hard to see the attracting periodic orbit because the "transient" part of the orbit messes up the picture.  A nice way to avoid this is to ignore the first part of the orbit to give it time to be attracted to the periodic orbit.  Another variation on the iteration commands makes this easy.  The command below plots iterates   x1000,..x1050 along the orbit, ignoring the first 1000 iterates.

\(a\  = \ 3.1; 
orb\  = \ f^{1000, 1050}[0.1]\)CobwebPlot[f, orb, {0, 1}];

Have a look at some of the other  a values using this technique.

Periodic Orbits via the equation f^k[x] = x

Of course iteration will only find attracting periodic orbits.  To find repellers you have to solve the equation f^k[x]=x.  You can find a formula for f^k[x] by just using the variable x as your initial condition. The first line below unsets   a so it is a variable again.

a =. f[x]f^2[x]f^3[x]

If f[x] is a polynomial, as it is for the logistic map, you can try to solve   f^k[x]=xby factorization , at least for small values of  k.

Factor[f^2[x] - x]Factor[f^3[x] - x]Factor[f^4[x] - x]

Unfortunately, the polynomials are still pretty big.  For specific values of   a you can use Mathematica's NSolve function to look for solutions.

\(a\  = \ 3.5; 
NSolve[f^2[x] - x == 0, x]\)NSolve[f^3[x] - x == 0, x]NSolve[f^4[x] - x == 0, x]

These results show where the period 1, 2 and 4 orbits are for this value of  a.  Apparently there are no period 3 orbits since the roots of  f^3[x]-x other than the fixed points are all complex numbers.  Try a=3.83,  which seemed to have a period 3 attractor in the plots from the previous section. You will actually find two  period-3 orbits -- one is a repeller which was not previously detected. Try   a=4 as well.

Another way to study periodic points  is to plot the graphs of the iterates and look for places where they cross the diagonal.  For this the built-in plotting function   Plot can be used.  You can plot several iterates on the same picture.  Plotting the 0-th iterate is one way to draw the "diagonal".

a = 3.5;
Plot[{f^0[x], f[x], f^2[x]}, {x, 0, 1}, 
        AspectRatio -> 1];
Plot[{f^0[x], f[x], f^3[x]}, {x, 0, 1}, 
        AspectRatio -> 1];

The graph of f^2[x] intersects the diagonal four times, at the fixed points and at two other points.  The latter form an orbit of period 2 (a repeller).  On the other hand, the graph of   f^3[x] crosses the diagonal only at the fixed points.  Now compare  a=3.83.

The abundance of periodic points for  a=4 is visible from the graphs of its iterates.  For example, the fifth iterate has a graph consisting of 32 monotonic pieces, each crossing the diagonal.

\(a = 4.0; 
Plot[{f^0[x], f[x], f^5[x]}, {x, 0, 1}, AspectRatio -> 1]; 
    \)

Try looking at some higher iterates of this map.

Finding Periodic Orbits and Multipliers Numerically

Another way to look for periodic points is to use a  special function included with this package called SeekPeriodicOrbit.  Instead of trying to find all possible solutions to the f^k[x]=x, it looks for one solution close to a given  point which you supply.  This is sometimes feasible when a more general solution is not, say when   f[x] is not a polynomial or when k is a bit larger.

We will try to find both of the period three orbits for the logistic map with    a=3.83.  Here is the graph of the third iterate and a blowup of that graph near the center

\(a = 3.83; 
Plot[{f^0[x], f^3[x]}, {x, 0, 1}, AspectRatio -> 1]; 

    Plot[{f^0[x], f^3[x]}, {x, 0.4, 0.6}, AspectRatio -> 1]; \)

From these graphs you can see some good guesses for where to search for the period 3 points.  In the blown up picture, the point on the left is on the attracting period 3 orbit and the point on the right is on the repeller.  We will use 0.5 as a starting point for a search.

orb1 =  SeekPeriodicOrbit[f, 3, 0.5]

This is apparently the attractor.  To find the repeller we search a bit to the right starting at 0.55.

orb2 =  SeekPeriodicOrbit[f, 3, 0.55]

Thus, using two different initial guesses leads to two different period 3 orbits which must be the ones we are after.  You can test them for stability using a package function called  Multiplier.

Multiplier[f, orb1]Multiplier[f, orb2]

This confirms that orb1 is an attractor and orb2 is a repeller.  Finally, we can display them both in a cobweb plot using red for the attractor and  blue for the repeller.

CobwebPlot[f, {orb1, orb2}, {0, 1}, LineColor -> {Red, Blue}];

This demonstrates another feature of the plotting functions in this package.  You can plot several orbits just by supplying a list of orbits as the argument instead of a single orbit.  Also there are various options which you can manipulate.  To get more information about a command you can type:

?CobwebPlot

Sensitive Dependence on Initial Conditions

In this section we will have a look at the dynamics of some dynamical systems which exhibit chaos.  A first sign of chaos is just the absence of a periodic attractor.  Here is a cobweb plot of the logistic maps with  a=4.0.

\(a = 4.0; 
orb\  = \ f^{1000, 1200}[0.1]; 
CobwebPlot[f, orb, {0, 1}]; 
    \)

In spite of the fact that we have skipped the first 1000 points on the orbit to give it a chance to converge to an attractor, there is no evidence of convergence -- the orbit seems to wander all over the interval [0,1].  Go back and try a similar plot for a=3.6785 and a=3.927737.  You will find similar behavior but on a slightly smaller interval.  Finally, have a look at a=3.57345 which seems to hop back and forth between two intervals of chaos.  Just from the pictures, it is not clear that these maps are really chaotic. Perhaps there is an attracting periodic orbit with a very large period so the resulting cobweb plot looks complicated.  For special values like a=4 we can prove that no attractors exist.  Also there is a theorem which says that there are many other such a values (a so-called "positive measure set" of them) without specifying exactly which ones.

It is easy to demonstrate sensitive dependence on initial conditions by simply comparing orbits with nearby initial conditions.  We will compute two such orbits and plot them together on the same graph.

a = 4.0;
orb1  =  f^{50}[0.1];
orb2 =  f^{50}[0.10001];
TimeSeriesPlot[{orb1, orb2}, 
        LineColor -> {Red, Green}];

You can see that after only about 12 or 13 iterations, the two orbits separate.  After that they are more or less indepedent, sometimes close and sometimes far apart.  Here is a cobweb plot demonstrating the same thing.  Perhaps we just need to take the initial conditions closer.  Go back and try changing the initial condition of orb2 to 0.1000001.  Keep adding zeroes and see what happens (machine precision is 15 places or so).

CobwebPlot[f, {orb1, orb2}, {0, 1}, 
        LineColor -> {Red, Green}];

If you try the same experiment for a parameter value with an attractor, the results are quite different.

a = 3.845;
orb1 =  f^{70}[0.1];
orb2 =  f^{70}[0.2];
TimeSeriesPlot[{orb1, orb2}, 
        LineColor -> {Red, Green}];

After a while the orbits get attracted to the period 3 attractor.  However, at any given iteration they are far apart because they are "out of phase".  In other words, one orbit has every third iterate x0,x3,x6,... converging to one point, say p0, on the period-3 orbit while the other orbit has x0,x3,x6,... converging to  p1.  However, this can be remedied by taking the two initial points closer.  Go  back up to the previous cell and try the initial condion 0.100001.  Now the long term behaviors are the same.

We can demonstrate exponential divergence and the idea of Lyapunov exponents by looking quantitatively at the differences bewteen two nearby orbits.    The plot below looks at the logaritbms of the absolute values of the differences  bewteen two nearby orbits for the chaotic logistic map  with a=4.

a = 4.0;
orb1 =  f^{50}[0.1];
orb2 =  f^{50}[0.10001];
diff  =  Abs[orb1[] - orb2[]]
ListPlot[Log[diff]];

Note that the first part of this plot looks roughly like a straight line.  We can estimate its slope as follows:

slope  =  (Log[diff[[11]]] - Log[diff[[1]]])/10

Compare this with the Lyapunov exponent.  This is approximated by the logarithm of the absolute value of the product of the derivatives along the orbit divided by the length of the orbit.  The product of the derivatives along the orbit is just what the function Multiplier is designed to compute.  It can also be applied to nonperiodic orbits.

mu  =  Multiplier[f, orb1]lyap  =  Log[ Abs[mu]]/50

Note that the lyapunov exponent is almost equal to the slope of the graph we found above. This shows how the Lyapunov exponent along the orbit orb1 measures the exponential rate at which nearby orbits (like orb2) diverge from it.

Some More Things to Try

If you have more time you could try some of the following experiments.  For these you will have to type the commands yourself.

1. Define a new dynamical system g = DS[Cos] based on the cosine function.  Plot graphs some of its iterates and explain what you find.

2. Using SeekPeriodicOrbit, find several periodic points of the map 4x(1-x).  Compute their multipliers using Multiplier.  The results are surprisingly simple.

3. Using iteration and plotting, convince yourself that the parameter value where the logistic map ax(1-x) first has an attracting period 2 orbit is a=3. Then by increasing a, find a slightly larger value where a period 4 attractor appears.  Continue to periods 8 and 16, if you can.  After finding these values ( call them a2, a4, a8, a16 ) have a look at the ratios  (a4-a2)/(a8-a4) and (a8-a4)/(a16-a8).  These approximate Feigenbaum's number.  The pattern continues.

Initialization

This section contains the definitions of a package of Mathematica functions for computing and plotting orbits of dynamical systems.  These functions will be used throughout the lab.  To evaluate them go up to the kernel menu and in the Evaluation submenu choose the command Evaluate Initialization.  It is not necessary to look at the detailed definitions in this section but each one is preceded by some words of explanation.

Dynamical Systems, Orbits and Iteration

One of the interesting aspects of Mathematica is that the various mathematical objects which arise (numbers, functions, etc.) are stored with a kind of label in front.  For example the number 5 is really Integer[5] internally.  The label Integer is what makes it possible to handle operations differently for integers than, say, for real numbers.  To make the same kind of distinctions possible in this package of functions, we are going to introduce new kinds of objects appropriate to dynamical systems theory.  Specifically we are going to use dynamical system objects, of the form DS[function here] and orbit objects, of the form Orbit[list here].  This extra layer of notation will usually not be visible while using this package.

A discrete dynamical system is just a mapping from some set to itself.  We will represent it in the form f = DS[f1] where f1 is the mapping.  For example we could have a mapping of the real line like f1[x_] := 4*x*(1-x) or a mapping of the plane to itself, such as f1[{x_,y_}]:=f[{x_,y_}]:={1.28-x^2+0.3*x,y}.  Then the command f = DS[f1] turns the mapping into a "dynamical system object".  You can think of DS as an operator which does nothing, except label the mapping as a dynamical system.

We can use this labelling to teach Mathematica how to iterate.  If f=DS[f1] is a dynamical system object, we want to be able to use the notation f[x] to evaluate the mapping and the notation f^n[x] to evaluate the n-th iterate of the mapping.  The following commands make this possible.  The mapping f1 is recovered from f as f[[1]], that is, the first element inside f. The notation DS/: in front of the definition of f^n assures that the definition is associated with the "operator" DS rather than with the exponentiation operator.  Iteration is basically the built-in function Nest.

Off[General::spell] 
Off[
      General::spell1]  (*first turn off annoying spelling messages*)

\(f_DS[x_]\  := \ \(f[\([1]\)]\)[x]\ \ 
\)
DS /:  f_DS^n_Integer[x_]  :=  Nest[f[[1]], x, n]

An orbit will be a list of states enclosed in an Orbit[] label.  We will use the notation F^{n}[x0] to represent the orbit Orbit[{x0,x1,...,xn}] of some initial point and F^{m,n} to represent Orbit[{xm,x1,...,xn}].  Orbits are computed using the built-in function NestList.

\(DS /: \ 
      f_DS^{n_Integer}[x_]\  := \ \(\(NestList[f[\([1]\)], x, 
          n]\)\(//\)\(Orbit\)\(
\)\)\)
DS /:  f_DS^{m_Integer, n_Integer}[x_]  :=  
      NestList[f[[1]], Nest[f[[1]], x, m], n - m] // 
        Orbit

Suppose that orb is an orbit object Orbit[{x0,x1,...,xn}].  We will introduce the notations orb[i] for xi,  orb[] for the list {x0,x1,...,xn} without the enclosing Orbit[ ] bracket,  orb[{m}] for the partial list {x0,x1,...,xm}, and orb[{m,n}] for {xm,x1,...,xn}.  Since orbits are indexed starting with x0 and Mathematica lists are indexed from index 1, an n+1 is needed.

\(l_Orbit[]\  := \ l[\([1]\)]
		\)
\(l_Orbit[n_Integer]\  := \ l[\([1, n + 1]\)]
\)
\(l_Orbit[{n_Integer}]\  := \ Take[l[\([1]\)], n + 1]
\)
l_Orbit[{m_Integer, n_Integer}]  :=  
      Take[l[[1]], {m + 1, n + 1}]

Sometimes we want to find the iterates or orbits of a list of initial states.  The following definitions extend the built-in Mathematica function Map to handle this situation.

Periodic Orbits

This group of functions is concerned with finding periodic orbits and computing their multipliers.  For the purposes of these functions, an orbit of period k will be represented by an Orbit of length k+1 of the form Orbit[{p0,p1,p2,....,pk}] with pk=p0.

The first function just tests an orbit to see if it seems to be periodic.  The test is just to see when it comes back with a certain small distance of itself.  For vectors the distance is measured using the maximum absolute value of the differences of the coordinates instead of with the Euclidean distance.

\(\(MinimalPeriod[orb_Orbit, tol_:  10^\((\(-6\))\)]\  := 
	
        With[{diff\  = \ 
              Map[Max[Abs[#]] &, 
                Drop[orb[] - Table[orb[0], {Length[orb[]]}], 1]]}, 
		Min[
            Position[
              diff, _?\((Function[d, d <= tol])\)]]]\)\(
\)\(		\ \)\)
MinimalPeriod::usage  =  "MinimalPeriod[orb,tol] \
tries to find the  minimal period of orb = Orbit[{x0,x1,...}] by searching \
for the an item xk in the sequence which agrees with x0 to the given \
tolerance.  The xi may be vectors."

If you have a periodic orbit you can compute the multipliers by multiplying the derivatives at all points of the orbit.  Mathematica can compute derivatives if you supply the formula and the following functions make use of this.  The first function just computes the derivative of a given function at a given point.  In the vector case, the derivative is the Jacobian matrix.  There are two versions of the function for that reason.

\(
ScalarQ[x_]\  := \ \(\(\(! \((Head[x] === List)\)\) // 
        LogicalExpand\)\(
\)\(	\)\)\)
\(DerAtPoint[f_DS, x0_?ScalarQ]\  := 
      Module[{x}, D[f[x], x] /. x -> x0]
\)
\(DerAtPoint[f_DS, x0_List]\  := 
      Module[{x, xx}, x\  = \ Array[xx, Length[x0]]; 
        Transpose[
          Map[\ D[f[x], #] &, \ x] /. 
            Table[xx[i] -> x0[\([i]\)], {i, 1, 
                Length[x0]}]]
		]\)

Next we multiply the derivatives at {p0,p1,...,p(k-1)} for a given periodic orbit {p0,p1,...,pk}.  For the vector case we need to use matrix multiplication, called Dot.

\(ScalarOrbitQ[
        orb_]\  := \ \(\(\((Head[orb] === Orbit)\) && \((Length[orb[]] >= 
              1)\) && ScalarQ[
            orb[1]]\)\(//\)\(LogicalExpand\)\(
\)\)\)
\(VectorOrbitQ[
        orb_]\  := \ \(\(\((Head[orb] === Orbit)\) && \((Length[orb[]] >= 
              1)\) && VectorQ[
            orb[1]]\)\(//\)\(LogicalExpand\)\(
\)\)\)
\(DerAlongOrbit[f_DS, orb_?ScalarOrbitQ] := \ 
	
      If[Length[orb[]] === 1, 2, 
        Apply[Times, 
          Map[DerAtPoint[f, #] &, 
            Drop[Reverse[orb[]], 1]]]
		]
\)
\(DerAlongOrbit[f_DS, orb_?VectorOrbitQ]\  := 
	
      If[Length[orb[]] === 1, 
        IdentityMatrix[Length[orb[1]]], \ 
		Apply[Dot, 
          Map[DerAtPoint[f, #] &, 
            Drop[Reverse[orb[]], 1]]]
		]\)

Finally, the function below computes the multiplier of a given periodic orbit.  In the vector case, it returns the eigenvalues of the product of the derivative matrices along the orbit.

\(Multiplier[f_DS, orb_?ScalarOrbitQ]\  := \ 
      DerAlongOrbit[f, orb]
\)
\(Multiplier[f_DS, 
        orb_?VectorOrbitQ]\  := \ 	
		With[{A = DerAlongOrbit[f, orb]}, 
        Eigenvalues[A]]
\)
Multiplier::usage  =  "Multiplier[f,orb] computes the multiplier or \
multipliers of a periodic orbit, orb = Orbit[{x0,...,xq}] of a dynamical \
system f by multiplying the derivatives of f at x0,...,x(q-1). For maps in \
dimension greater than 1 the derivatives will be matrices and Multiplier \
returns the list of eigenvalues."

Next we have a fairly complicated function called SeekPeriodicPoint which searches for periodic points.  It tries to solve the equation f^k(x) = x using Mathematica's built-in FindRoot command which is based on Newton's method.  The search is conducted starting at a given point, but if that fails nearby points are chosen at random.  So we need to define a RandomPoint function first, which produced a number or vector chosen at random from a given range.  In the scalar case, the range should be an interval given as a vector {xmin,xmax}.  In the vector case the range will be specified by a list of two vectors -- a vector of minima and a vector of maxima.  For example {{xmin,ymin,zmin},{xmax,ymax,zmax}} represents a cube in three dimensions.

\(\(RandomPoint[r_?VectorQ]\  := Random[Real, r];\)\(
\)\)
\(\(RandomPoint[r : {{___}, {___}}]\  := 
        Map[Random[Real, #] &, Transpose[r]]\)\(
\)\)
RandomPoint::usage  =  "RandomPoint[range] returns \
a random vector in the given range.  For an n-dimensional point, the range \
should be a pair of the form {{x1min,x2min,...},{x1max,x2max,...}}."

There are several options for SeekPeriodicPoint which constrain the search procedure.

\(\(Options[SeekPeriodicOrbit]\  = \ {
	
              Radius -> Automatic, 
	MaxTries -> 10, 
	
              Tolerance -> 0.000001, 
	MaxIterations -> 100
	};\)

        \(spooptpattern\  = \ \((Radius -> _)\) | \((MaxTries -> _)\) | \((Tolerance -> _)\);\)



        SeekPeriodicOrbit[f_DS, k_Integer, xguess_, opts___]\  := \ 
          Module[{fk, guess, sol, x, xx, x0, po, ic, mt, mi, tol, ra, xrange, 
              fropts}, 
		mt\  = \ \(MaxTries /. {opts}\) /. 
                Options[
                  SeekPeriodicOrbit]; 
		tol\  = \ \(Tolerance /. {opts}\) /. Options[SeekPeriodicOrbit]; 
		ra\  = \ \(Radius /. {opts}\) /. 
                Options[SeekPeriodicOrbit]; 
		fropts\  = \ 
              DeleteCases[{opts}, spooptpattern]; 
		 (*\ 
              Try\ the\ given\ initial\ point\ first\ *) 
		guess\  = \ 
              xguess; 
		If[
              VectorQ[guess], 
			x = Array[xx, Length[guess]]; 
              ic = Transpose[{x, guess}], 
			ic = {{x, 
                    guess}}]; 
		fk\  = \ 
              f^k[x] - x // 
                Simplify; 
		sol\  = \ \(FindRoot[fk, ##] &\)\  @@ \ 
                Flatten[{ic, fropts}, 1]; 
		x0\  = \ 
              x /. sol; 
		po\  = \ f^{k}[x0]; 
		If[
              MinimalPeriod[po] === k, Return[po]]; 
		 (*\ 
              Failing\ this, \ try\ random\ guesses\ nearby\ *) 
		If[
              ra === Automatic, 
			ra = 
                2*Max[Abs[xguess - f^k[xguess]]]; 
			 (*\(ra\  = \ {0.1, 
                      0.1};\)*) 
			xrange\  = \ {xguess - ra, 
                  xguess + ra}]; 
		While[
              mt > 1, 
			guess\  = \ RandomPoint[xrange]; 
			If[
                VectorQ[guess], 
				ic = 
                  Transpose[{x, guess}], 
				ic = {{x, 
                      guess}}
				]; 
			x0\  = \ 
                x /. \((\(FindRoot[fk, ##] &\)\  @@ 
                      Flatten[{ic, fropts}, 1])\); 
			po\  = \ 
                f^{k}[x0]; 
              Print[po]; 
			If[MinimalPeriod[po] === k, 
                Return[po]]; 
			\(mt--\)]; 
		 (*\ Failing\ this, \ 
              return\ the\ empty\ list\ *) 
		Return[Orbit[{}]];
		]\)

\(SeekPeriodicOrbit::usage\  = \ "SeekPeriodicOrbit[f,k,xguess,opts] tries to find a periodic orbit of period k for the dynamical system f near the point xguess using the Mathematica function FindRoot. Works well only for small k and low dimensional systems.  If it succeeds,it returns Orbit[{x0,x1,...,xk}] of length k+1.  If it fails,it returns Orbit[{}]."\)

Unstable Manifolds

Here is a function for computing unstable manifolds of 2D maps.  Supposing that we have already located a period q point x0, we first test it to make sure it is a saddle point.  Then we locate the unstable eigenvector, construct a short line segment in that direction, subdivide it into many smaller segments and then extend these by iteration of the dynamical system.

\(IsSaddle[f_DS, x0_, q_]\  := \ 
      With[
	{mu\  = \ Multiplier[f, f^{0, q}[x0]]}, 
	
        If[Im[mu] != {0, 0}, Return[False], ]; 
	
        If[	Abs[mu[\([1]\)]] > 1 && 
              Abs[mu[\([2]\)]] < 1 || 
	\ \ 	Abs[mu[\([2]\)]] > 1 && 
              Abs[mu[\([1]\)]] < 1, True, False]
	]
	\)
\(Subdivide[l_, i_, 
        n_:  2]\  := \(\(Take[l, i - 1]\)\(~\)\(Join\)\(~\)\(Table[
          N[l[\([i]\)]*\((1 - t/n)\) + l[\([i + 1]\)]*t/n], {t, 0, 
            n}]\)\(~\)\(Join\)\(~\)\(Drop[l, i + 1]\)\(
\)\)\)
\(ListOrbit[f_, l_List, n_]\  := \ 
      Flatten[Transpose[
			Map[\ NestList[f, #, n] &, l]], 
        1]
\)
\(\(Options[UnstableManifold]\  = \ {
	Period -> 1, 
	
          Subdivisions -> 20, 
	Radius -> 0.01};\)
\)
\(UnstableManifold[f_DS, x0_, n_, opts___]\  := \ 
	
      Module[{q, su, ra, A, mu, ev, vu, l, l1, l2, g}, 
	
        q\  = \ \(Period /. {opts}\) /. Options[UnstableManifold]; 
	
        su\  = \ \(Subdivisions /. {opts}\) /. Options[UnstableManifold]; 
	
        ra\  = \ \(Radius /. {opts}\) /. Options[UnstableManifold]; 
	
        If[\(! IsSaddle[f, x0, q]\), Print["Not a Saddle"]; 
          Return[{}]]; 
	A\  = \ DerAlongOrbit[f, f^{0, q}[x0]]; 
	
        mu\  = \ Eigenvalues[A]; 
	ev\  = \ Eigenvectors[A]; 
	
        vu\  = \ If[Abs[mu[\([1]\)]] > 1, ev[\([1]\)], 
            mu[\([1]\)] = mu[\([2]\)]; ev[\([2]\)]]; 
	
        g\  = \ If[mu[\([1]\)] > 0, f^q[#] &, f^\((2*q)\)[#] &]; 
	
        l1\  = \ ListOrbit[g, Subdivide[{x0 + ra*vu, g[x0 + ra*vu]}, 1, su], 
            n]; 
	l2\  = \ 
          ListOrbit[g, Subdivide[{x0 - ra*vu, g[x0 - ra*vu]}, 1, su], n]; 
	
        Orbit[Reverse[l2]~Join~l1]
	]
	\)
UnstableManifold::usage  =  "UnstableManifold[f,x0,n,opts] returns \
a list of points on the unstable manifold of the periodic point x0 for the \
map f. By default, x0 is a fixed point; period q can be handled by setting \
the option Period->q. A number of points are chosen along the unstable \
eigenvector close to x0. The number of points is determined by the option \
Subdivisions which defaults to 20. The distance from x0 is controlled by the \
option Radius which defaults to 0.01. The argument n is the number of \
iterations used to extend the manifold outside this radius."

Plotting

This section contains some functions for plotting the orbits of dynamical systems in interesting ways.  There are many different types of plots, all of which take similar kinds of options, such as the ranges for the variables and the sizes and colors of the dots and lines.  So all the plotting functions below are based on a GenericPlot function which handles all of these details.  Then the only difference between the various plots is some preliminary transformations which they perform on the orbits.  And you can even make up new types of plots by defining new transformations.

Generic Plotting Functions

First there are some functions for turning a list of 2D or 3D vectors into a Mathematica graphics object -- either a list of Point 's, a single Line or a Polygon.  There are optional arguments to control the color and size.  There are two versions of each function, one for 2D and one for 3D.

ScalarListQ[
          l_]  :=  (Head[l] === List) && (Length[l] >= 1) && 
            ScalarQ[l[[1]]] // LogicalExpand

TwoDListQ[
          l_]  :=  (Head[l] === List) && (Length[l] >= 
                1) && (Head[l[[1]]] === 
                List) && (Length[l[[1]]] === 2) // LogicalExpand

ThreeDListQ[
          l_]  :=  (Head[l] === List) && (Length[l] >= 
                1) && (Head[l[[1]]] === 
                List) && (Length[l[[1]]] === 3) // LogicalExpand


\(\(MakePoints[l_?TwoDListQ, color_: {0, 0, 0}, 
          size_:  0.01]\)\(\ \)\(:=\)\(
\)\(	\)\(Graphics[{Apply[RGBColor, 
              color], PointSize[size], Map[Point, l]}]\)\(
\)\)
\(\(MakeLine[l_?TwoDListQ, color_: {0, 0, 0}, 
          thickness_:  
            0.005]\)\(\ \)\(:=\)\(
\)\(	\)\(Graphics[{Apply[RGBColor, 
              color], Thickness[thickness], Line[l]}]\)\(
\)\)
\(\(MakePolygon[l_?TwoDListQ, color_: {0, 0, 0}, 
          size_:  0.01]\  := 
	
        Graphics[{Apply[RGBColor, color], PointSize[size], 
            Polygon[l]}]\)\(
\)
      \(	\)\(
\)\( (*\ 3  D\ versions\ *) \)\(									\)\)
\(\(MakePoints[l_?ThreeDListQ, color_: {0, 0, 0}, 
          size_:  0.01]\)\(\ \)\(:=\)\(
\)\(	\)\(Graphics3D[{Apply[
              RGBColor, color], PointSize[size], Map[Point, l]}]\)\(
\)\)
\(\(MakeLine[l_?ThreeDListQ, color_: {0, 0, 0}, 
          thickness_:  
            0.005]\)\(\ \)\(:=\)\(
\)\(	\)\(Graphics3D[{Apply[RGBColor, 
              color], Thickness[thickness], Line[l]}]\)\(
\)\)
\(MakePolygon[l_?ThreeDListQ, color_: {0, 0, 0}, 
          size_:  0.01]\  := 
	
        Graphics3D[{Apply[RGBColor, color], PointSize[size], Polygon[l]}]\)

Here are a few definitions of popular RGB colors so you don't have to supply explicit R,G,and B values.

\(\(ColorTable\  = \ {Black -> {0, 0, 0}, White -> {1, 1, 1}, 
          Red -> {1, 0, 0}, 
	Green -> {0, 1, 0}, Blue -> {0, 0, 1}, 
          Yellow -> {1, 1, 0}, Purple -> {0.5, 0, 1}, 
	
          BlueGreen -> {0, 1, 1}, LightGray -> {0.75, 0.75, 0.75}, 
          Gray -> {0.5, 0.5, 0.5}, 
	
          DarkGray -> {0.25, 0.25, 0.25}};\)\)

Some of the plots included decorations like a copy of a graph of a function, a diagonal line, or a circle.  These also need to be converted to Mathematica graphics objects.

\(MakeDiagonal[x_, color_: {0, 0, 0}, 
          thickness_:  0.005]\  := \ 
	
        Graphics[{Apply[RGBColor, color], 
            Thickness[
              thickness], 
										Line[{{x[\([1]\)], 
                  x[\([1]\)]}, {x[\([2]\)], x[\([2]\)]}}]}]\)

\(MakeCircle[rad_:  1.0, pts_:  100, color_: {0, 0, 0}, 
            thickness_:  0.005]\  := 
	
          With[{d\  = \ N[2  Pi/pts]}, 
            Graphics[{Apply[RGBColor, color], Thickness[thickness], 
                Line[Table[{Cos[d*i], Sin[d*i]}, {i, 0, pts}]]}]]\)
										
\(MakeFunctionGraph[f_, x_, 
        color_: {0, 0, 0}, thickness_:  0.005]\  := \ 
	
      With[{r\  = \ 
            Range[x[\([1]\)], x[\([2]\)], 
              0.005*\((x[\([2]\)] - x[\([1]\)])\)]}, 
        Graphics[{Apply[RGBColor, color], Thickness[thickness], 
            Line[Transpose[{r, Map[f, r]}]]}]]\)

Here is the generic plotting function.  The first version is for plotting a single orbit.  There are numerous options

\(\(Options[
          GenericPlot]\  = \ {
		Points -> True, 
		Lines -> 
            True, 
		PointColor -> Black, 
		LineColor -> 
            Black, 
		PointSize -> 0.01, 
		LineSize -> 
            0.002, 
		Transformation -> 
            IdentityTransformation, 
		MoreGraphics -> {}
		};\)\)
\(\(genoptpattern\  = \ \((Points -> _)\) | \((Lines -> _)\) | \((PointColor -> _)\) | 
		\((LineColor -> _)\) | \((PointSize -> _)\) | \((LineSize -> _)\) | 
		\((Transformation -> _)\) | \((MoreGraphics -> _)\);\)
\)
\(GenericPlot[orb_, 
        opts___]\  := \ 
		Module[{showopts, pf, lf, pc, lc, ps, ls, pr, 
          tr, mg, 
				xmin, xmax, trorb, 
          gr}, 
			pf\  = \ \(Points /. {opts}\) /. 
            Options[GenericPlot]; 
			pc\  = \ \(\(PointColor /. {opts}\) /. Options[GenericPlot]\) /. 
            ColorTable; 
			ps\  = \ \(PointSize //. {opts}\) /. 
            Options[GenericPlot]; 
			lf\  = \ \(Lines /. {opts}\) /. 
            Options[GenericPlot]; 
			lc\  = \ \(\(LineColor /. {opts}\) /. 
              Options[GenericPlot]\) /. 
            ColorTable; 
			ls\  = \ \(LineSize /. {opts}\) /. 
            Options[GenericPlot]; 
			pr\  = \ 
          PlotRange /. {opts}; 
			tr\  = \ \(Transformation /. {opts}\) /. 
            Options[GenericPlot]; 
			mg\  = \ \(MoreGraphics /. {opts}\) /. Options[GenericPlot]; 
			showopts\  = \ 
          DeleteCases[{opts}, genoptpattern]~
            Join~
						{Axes -> True, 
              AspectRatio -> 1}; 
			xmin\  = \ Min[orb[]]; 
        xmax\  = \ Max[orb[]]; 
			xmin\  -= \ 0.1*\((xmax - xmin)\); 
        xmax\  += \ 0.1*\((xmax - xmin)\); 
			If[Head[pr] == List, 
          xmin = pr[\([1, 1]\)]; 
          xmax = pr[\([1, 2]\)], ]; 
			trorb\  = \ 
          tr[orb]; 
			gr\  = \ {
				If[lf, 
                MakeLine[trorb, lc, ls], {}], 
				If[pf, 
                MakePoints[trorb, pc, ps], {}]
				}~Join~
            mg; 
			Show[gr, showopts]
			]
\)

The default identity transformation is defined here.

IdentityTransformation[orb_]  :=  orb[]

Here is the version for plotting a list of orbits.

OrbitListQ[
          l_]  := (Head[l] === List) && (Length[l] >= 
                1) && (Head[l[[1]]] === Orbit) // LogicalExpand

\(\(GenericPlot[orbs_?OrbitListQ, 
          opts___]\)\(\ \)\(:=\)\(\ \)\(
\)\(		\)\(Module[{num, showopts, 
            pf, lf, pc, lc, ps, ls, pr, tr, mg, 
				xmin, xmax, trorbs, 
            gr}, 
			pf\  = \ \(Points /. {opts}\) /. 
              Options[GenericPlot]; 
			pc\  = \ \(\(PointColor /. {opts}\) /. Options[GenericPlot]\) /. 
              ColorTable; 
			ps\  = \ \(PointSize //. {opts}\) /. 
              Options[GenericPlot]; 
			lf\  = \ \(Lines /. {opts}\) /. 
              Options[GenericPlot]; 
			lc\  = \ \(\(LineColor /. {opts}\) /. Options[GenericPlot]\) /. 
              ColorTable; 
			ls\  = \ \(LineSize /. {opts}\) /. 
              Options[GenericPlot]; 
			pr\  = \ 
            PlotRange /. {opts}; 
			tr\  = \ \(Transformation /. {opts}\) /. Options[
                GenericPlot]; 
			mg\  = \ \(MoreGraphics /. {opts}\) /. 
              Options[GenericPlot]; 
				 (*\ 
            make\ the\ options\ into\ lists\ if\ necessary\ *) 
			num\  = \ Length[orbs]; 
			If[Head[pf] =!= List, 
            pf = Table[pf, {num}]]; 
			If[VectorQ[pc], 
            pc\  = \ Table[pc, {num}]]; 
			If[Head[ps] =!= List, 
            ps = Table[ps, {num}]]; 
			If[Head[lf] =!= List, 
            lf = Table[lf, {num}]]; 
			If[VectorQ[lc], 
            lc = Table[lc, {num}]]; 
			If[Head[ls] =!= List, 
            ls = Table[ls, {num}]]; 
			showopts\  = \ 
            DeleteCases[{opts}, genoptpattern]~
              Join~
						{Axes -> True, 
                AspectRatio -> 1}; 
			xmin\  = \ Min[Union[orbs]]; 
          xmax\  = \ Max[Union[orbs]]; 
			xmin\  -= \ 
            0.2*\((xmax - xmin)\); 
          xmax\  += \ 0.2*\((xmax - xmin)\); 
			If[Head[pr] == List, 
            xmin = pr[\([1, 1]\)]; 
            xmax = pr[\([1, 2]\)], ]; 
			trorbs\  = \ 
            Map[tr, orbs]; 
			gr\  = \ 
            Flatten[{
				MapThread[
                    If[#2, MakeLine[#1, #3, #4], {}] &, 
								{trorbs, lf, lc, ls}], 
				MapThread[
                    If[#2, 
                        MakePoints[#1, #3, #4], {}] &, 
								{trorbs, pf, pc, ps}]
				}]~Join~mg; 
			Show[gr, 
            showopts]
			]\)\(
\)\)

Plots for One-dimensional Orbits

Each plot is mainly determined by the preliminary transformation it does to its orbits.  The CobwebPlot transforms a scalar orbit  Orbit[{x0,x1,...}] to a list of pairs {{x0,x0},{x0,x1},{x1,x1},{x1,x2},...} suitable for plotting in the plane.  It also contains a copy of the diagonal and a copy of the graph of the dynamical system, added with the MoreGraphics option of GenericPlot.

\(CobwebTransformation[orb_?ScalarOrbitQ]\  := \ 
        With[{dup = 
              Flatten[Transpose[{orb[], 
                    orb[]}]]}, 
					Transpose[{Drop[dup, \(-1\)], 
              Drop[dup, 1]}]]\)

\(\(CobwebPlot[f_DS, orb_, xrange_, opts___]\  := \ 
	
        GenericPlot[orb, opts, Transformation -> CobwebTransformation, 
          Points -> False, 
		PlotRange -> {xrange, 
              xrange}, 
		PlotLabel -> 
            StringForm["f = ``", 
              Evaluate[f["x"]]], 
		MoreGraphics -> {MakeDiagonal[
                xrange], 
						MakeFunctionGraph[f, 
                xrange]}
				]\)\(
\)\(				\)\)
CobwebPlot::usage  =  "CobwebPlot[f,orb,xrange,opts] \
produces a  graphical analysis or cobweb plot showing the graph of the \
function f, the diagonal and the orbit orb.  xrange is a list giving the left \
and right boundaries of the plot.  opts is any  number of options.  A list of \
orbits can be given instead of a  single orbit.  For example if f is a \
dynamical system and orb1 and orb2 are orbits which have already been \
computed, then  CobwebPlot[f,{orb1,orb2},{0.0,1.0},LineColor->{Black,Red}] \
will make a plot showing the two orbits in the given colors."

Another useful plot for scalar dynamical systems is the time series plot.  This time the scalar orbit Orbit[{x0,x1,...}] is transformed into the list of pairs {{0,x0},{1,x1},{2,x2},...}

TimeSeriesTransformation[orb_?ScalarOrbitQ]  :=  
        Map[Flatten, Transpose[{Range[0, Length[orb[]] - 1], orb[]}]]


\(\(TimeSeriesPlot[orb_, opts___]\  := \ 
	
        GenericPlot[orb, opts, Transformation -> TimeSeriesTransformation, 
          PointSize -> 0.02, AspectRatio -> 0.5, 
          AxesLabel -> {"n", "x"}]\)\(
\)\(				\)\)
TimeSeriesPlot::usage  =  "TimeSeriesPlot[orb,opts] \
produces a time series plot showing the graph of the orbit orb  as a function \
of time index n. opts is any number of options. A list of orbits can be given \
instead of a single orbit.   For example if orb1 and orb2 have already been \
computed  using Orbit[...], then  \
TimeSeriesPlot[{orb1,orb2},LineColor->{Black,Red}] will make a plot showing \
the two orbits in the given colors."

The final plots for scalar data are two and three dimensional delay plots.  In the two-dimensional case, the scalar orbit Orbit[{x0,x1,...}] is transformed into the list of pairs {{x0,xd},{x1,x(d+1)},{x2,x(d+2)},...} where d is a delay which defaults to d=1.

\(DelayTransformation[orb_?ScalarOrbitQ, dim_:  2, 
          delay_:  1]\  := \ 
	
        With[{ln = Length[orb[]] - \((dim - 1)\)*delay}, 
          Map[Flatten, 
            Transpose[
              Table[Take[orb[], {1 + i*delay, ln + i*delay}], {i, 0, 
                  dim - 1}]]]
	]\)
		
\(\(DelayPlot[orb_, dim_Integer:  2, delay_Integer:  1, 
          opts___]\  := 
	
        Switch[dim, \ 
		2, 
          GenericPlot[orb, opts, 
            Lines -> 
              False, 
			Transformation -> \((DelayTransformation[#, dim, 
                    delay] &)\), 
			AxesLabel -> {"x[n]", "x[n+1]"}], 
		3, 
          GenericPlot[orb, opts, 
            Lines -> 
              False, 
			Transformation -> \((DelayTransformation[#, dim, 
                    delay] &)\), 
			AxesLabel -> {"x[n]", "x[n+1]", "x[n+2]"}]
			]\)\(
\)\(			\)\)
DelayPlot::usage  =  "DelayPlot[orb,dim,delay,opts] \
produces a  delay plot showing the orbit orb embedded in dim dimensions, \
using a given delay. opts is any number of options. dim defaults to 2 and \
delay to 1 so, for example, the command DelayPlot[orb] will make a plot \
showing {orb_n,orb_n+1} in the plane. A list of orbits can be given instead \
of a single orbit."

\(CircleOrbit[orb_]\  := \ 
      With[{tp\  = \ N[2  Pi]}, 
						Map[{Cos[tp*#], Sin[tp*#]} &, 
          orb]]\)

Plots for Vector Orbits

To plot a vector orbit, just project out two of the coordinates (the first two by default) and plot some dots.

\(ProjectionTransformation[orb_?VectorOrbitQ, 
          indices_: {1, 2}]\  := \ 
Map[Part[#, indices] &, orb[]]\)

\(PhasePlot[orb_, opts___]\  := 
	
        GenericPlot[orb, opts, 
          Transformation -> ProjectionTransformation, 
				PointSize -> 
            0.007, AspectRatio -> 1, 
				Lines -> False]\)

Bifurcation Plot

\( (*\ A\ function\ for\ plotting\ attractor\ bifurcation\ 
	
          diagrams\ for\ one - 
        parameter\ families\ of\ maps\ *) 
	
\(Options[
          BifurcationPlot]\  = \ {
	ParameterPoints -> 100, 
	
          OrbitPoints -> 100, 
	Delay -> 100, 
	ClearOut -> True, 
	
          RoundDigits -> 8, 
	
          PlotRange -> All, 
\ \ PointSize -> 
            0.001, 
	
          Transformation -> IdentityTransformation};\)
	\)
\(\(bpoptpattern\  = \ \((ParameterPoints -> _)\) | \((OrbitPoints -> _)\) | 
	\((Delay -> _)\) | \((ClearOut -> _)\) | \((RoundDigits -> _)\) | \((PointSize -> _)\) | \((Transformation -> _)\);\)
	\)
\(BifurcationPlot[f_DS, ic_, pr_, opts___]\  := 
	
      Module[{showopts, pp, op, de, co, rd, plr, ps, tr, fc, orb, Rnd, Clp, 
          ptlist = {}}, 
		pp\  = \ \(ParameterPoints /. {opts}\) /. 
            Options[BifurcationPlot]; 
		op\  = \ \(OrbitPoints /. {opts}\) /. Options[BifurcationPlot]; 
		de\  = \ \(Delay /. {opts}\) /. 
            Options[BifurcationPlot]; 
		co\  = \ \(ClearOut /. {opts}\) /. 
            Options[BifurcationPlot]; 
		rd\  = \ \((\(RoundDigits /. {opts}\) /. Options[
                  BifurcationPlot])\)*1.0; 
		plr\  = \ \(PlotRange /. {opts}\) /. Options[
              BifurcationPlot]; 
	\ \ \ ps\  = \ \(PointSize /. {opts}\) /. 
            Options[BifurcationPlot]; 
		tr\  = \ \(Transformation /. {opts}\) /. Options[BifurcationPlot]; 
		showopts\  = \ 
          DeleteCases[{opts}, bpoptpattern]~
            Join~
						{Axes -> True, AspectRatio -> 1}; 
		If[co, 
          Unprotect[Out]; Clear[Out]; 
          Protect[Out]]; 
		 (*\(fc\  = \ 
              Compile[{x, c}, Evaluate[f[x, c]]];\)*) 
		Rnd\  = \ 
          Function[x, Round[10^rd*x]*10^\((\(-rd\))\)]; 
		Clp\  = \ 
          If[ListQ[plr], 
			Function[l, 
              Select[l, \((plr[\([1]\)] <= #\  && \ # <= 
                        plr[\([2]\)])\) &]], 
			Function[l, 
              l]]; 
		For[c = pr[\([1]\)], c <= pr[\([2]\)], 
          c = c + \((pr[\([2]\)] - pr[\([1]\)])\)/pp, 
			orb\  = \ 
            tr[f^{de, de + op}[ic[c]]]; 
			orb\  = \ 
            Clp[Union[Rnd[orb]]]; 
			ptlist\  = \ 
            ptlist~Join~
				Transpose[{Table[c, {Length[orb]}], 
                  orb}];
			]; 
		 (*\(Print[
              Length[ptlist]];\)*) 
		Show[
          Graphics[{PointSize[ps], 
				Map[Point, ptlist]}], 
          showopts];
		]
			\)
\(BifurcationPlot::usage\  = \ "BifurcationPlot[f,ic,pr,opts] produces an orbit diagram for a one parameter family of maps depending on a parameter, c. Here f should be a dynamical system object whose underlying formula contains the variable c. ic[c] is a function of one variable specifying the initial point to use in iterating. pr is the range of c values to use.  There are several options.  For example, if f has been defined as the logistic map f= DS[f1] where f1[x_] := c*x*(1-x) and ic[c_] := 0.5 the command: BifurcationPlot[f,ic,{0,4},ParameterPoints->300,OrbitPoints->200,Delay->100, PlotRange->{0,1}] will produce an orbit diagram with the parameter c ranging over the interval {0,4} in 300 steps.  For each step an orbit with initialcondition ic[c] is computed. The first 100 points are not plotted and the next 200 are. Only points lying in the range 0<=x<=1 are plotted. This function is memory intensive."
		\)

Iterated Function Systems

Here we define an IFS object representing an iterated function system.  To apply and IFS to a point we simply apply each function to the point and return the list of answers.  To iterate this we need to define what it means to apply and IFS to a list of vectors.  Finally ,we define how to use an IFS to transform a line or a list of lines

\(SwitchHead[f_, h_]\  := \ 
      f /. \(\(f[\([0]\)]\)\(->\)\(h\)\(
\)\)\)
f_IFS[x : {__Line}]  := Flatten[Map[f, x], 1]
f_IFS[x : {__List}]  :=  Flatten[Map[f, x], 1]
f_IFS[x_List /; VectorQ[x]]  :=  
      Map[#[x] &, SwitchHead[f, List]]
\(f_IFS[Line[l_]]\  := \ Map[Line, Transpose[Map[f, l]]]
\)
\(IFS /: \ f_IFS^n_Integer[x_]\  := \ Nest[f, x, n]
\)
IFS /:  f_IFS^n_Integer[x : {__Line}]  :=  Nest[f, x, n]