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.83; 
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.

Lab 2 --- Bifurcations and Chaos

In this lab you will have a look at some bifurcation diagrams and then try out some symbolic dynamics both for dynamical systems on an interval and on an invariant Cantor set.  As before, the 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 Critical Orbit Diagrams.  The bracket will expand to reveal Mathematica cells containing explanatory text and commands.  Just work your way through the section reading the explanations and executing the commands (you do this by clicking in the cells and then hitting the Enter key or Shift-Return).  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 !

Critical Orbit Diagrams

In this section we will have a look at the amazing details of the orbit diagram for the logistic map.  First we define the logistic dynamical system.  The bifurcation plotting software is set up to vary a parameter called c so we will use that letter in the formula.

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

To construct an orbit diagram we need to start with some initial condition.  In general this could vary with c so we define an initial condition function ic[c] which gives the desired starting point.  In this case it is a constant function at the critical point x0 = 0.5.

ic[c_]  :=  0.5

Now the function bifurcation plot can be used to draw the orbit diagram. Here is some information about this function, whose actual definition is contained in the initialization section of this notebook.

?BifurcationPlot

BifurcationPlot[f, ic, {0, 4}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, 
      PointSize -> 0.001, 
PlotRange -> {0, 1}]

If you don't see anything in this plot you may have to increase the PointSize parameter by editing the previous line.  For example you could try PointSize->0.005 or PointSize->0.01.  Also you can use the mouse to expand the size of the picture.

	Of course the most interesting part of the orbit diagram is from c=2 to c=4.  You can zoom in on that part by changing the range for c.

BifurcationPlot[f, ic, {2, 4}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, PointSize -> 0.001, 
      PlotRange -> {0, 1}]

It is endlessly amusing to blow up small regions of the diagram and to ponder the incredible complexity contained in it.  For example here is a tiny window where a period 5 orbit is born in a saddle-node bifurcation and then proceeds to period double producing attractors of period 10, 20 , 40 , etc.

BifurcationPlot[f, ic, {3.73, 3.75}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, PointSize -> 0.001, 
      PlotRange -> {0, 1}]

You can blow up just one of the five period-doubling regions by changing the   PlotRange as well as the range for c.  Here a blow-up of the part of the last diagram near the center revealing further smaller windows.  For deep blow-ups it is sometimes good to increase the Delay to give more time to converge to the attractor and the number of OrbitPoints, since not all of them will appear in the chosen  Plot Range.

BifurcationPlot[f, ic, {3.742, 3.745}, ParameterPoints -> 300, 
      OrbitPoints -> 400, Delay -> 200, PointSize -> 0.001, 
      PlotRange -> {0.45, 0.55}]

I suggest that you go back to the previous cell and  explore other regions of the orbit diagram by making changes to the c range and the Plot Range.

Of course the logistic map is not the only interesting case.  Here is the critical orbit diagram for the quadratic maps f(x) = c-x^2.  For these, the critical point is at x=0 instead of at x=1/2.

g1[x_]  :=  c - x^2
g = DS[g1];
ic[c_]  :=  0.0;  
BifurcationPlot[g, ic, {-0.25, 2}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, 
      PointSize -> 
        0.001, 
PlotRange -> {-2, 2}]

And here a collection of multiples of the cosine function.  These maps have more than one critical point but we can just pick one, say x=0,  and make the plot.

h1[x_]  :=  c*Cos[x]
h = DS[h1];
ic[c_]  :=  0.0;  
BifurcationPlot[h, ic, {1, 3}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, 
      PointSize -> 
        0.001, 
PlotRange -> {-Pi, Pi}]

Symbolic Dynamics and Itineraries

In this section we will look at symbolic dynamics for the logistic map f(x) = 4x(1-x) and the corresponding tent map T(x).  First we define the maps.

f1[x_]  :=  4*x*(1 - x)
f =  DS[f1]
T1[x_]  :=  If[x < 1/2, 2*x, 2 - 2*x]
T  =  DS[T1]

Next we define a function to produce the itinerary of an orbit.  To do this we just test whether or not each point is less than or equal to 1/2.

IntervalTest[x_]  :=  
      If[x ≤ 1/2, 0, 1]
Itinerary[Orbit[pts_]] :=  Map[IntervalTest, pts]

Now lets test it on an orbit of the tent map, say the one with x0 = 5 / 31 .

tentorb  =  T^{0, 25}[5/31]
Itinerary[tentorb]
TimeSeriesPlot[tentorb];

Recall that the tent map and the logistic map are conjugate via the change of coordinates y = sin^2(π*x/2).  So if we compute the logistic orbit with initial condition y0 = sin^2(5π/62) we should get the same itinerary.  For this we will have to use approximate numbers so we set the initial condition to be N[sin^2(5π/62)]

y0  =  N[Sin[5*Pi/62]^2]
logorb  =  f^{0, 25}[y0]
Itinerary[logorb]
TimeSeriesPlot[logorb];

Due to the numerical approximation, this orbit does not repeat exactly.  Eventually, the errors accumulate and the itinerary becomes different. This is a reflection of the sensitive dependence on initial conditions, of course.

logorb  =  f^{0, 100}[y0]
Itinerary[logorb]
TimeSeriesPlot[logorb];

Now we will turn the problem around seek out orbits with a given itinerary (or at least a finite part of an itinerary).  A given finite sequence {a0,a1,...,an} is realized by a closed interval of initial points.  We can compute this interval by using the two branches of the inverse map.  One branch maps [0,1] into the interval I0 = [0,1/2] and the other into the interval I1 = [1/2,1].  By applying these to the initial points 0 and 1 in the backwards order {an,...,a1,a0} we can find the endpoints of our closed interval !  Since the tent map can be handled using exact arithmetic we do that first and then use the conjugacy to handle the logistic map.

TentInv[0][x_]  :=  x/2
TentInv[0][{x1_, x2_}] := {x1/2, x2/2}
TentInv[1][x_]  :=  1 - x/2
TentInv[1][{x1_, x2_}] := {1 - x1/2, 
        1 - x2/2}
TentInterval[itin_]  :=  
      Module[{int = {0, 1}}, 
For[i = Length[itin], 
          i ≥ 
            1, i--, 
int  =  TentInv[
                itin[[i]]][int]]; 
If[
          int[[1]] < int[[2]], Return[int], 
          Return[Reverse[int]]]]

For example, here are all the interval corresponding to itineraries of length 2.

TentInterval[{0, 0}]
TentInterval[{0, 1}]
TentInterval[{1, 0}]
TentInterval[{1, 1}]

To get the intervals for the logistic map we apply the conjugacy and numerical approximation.  We will set the precision of the numerical approximation so as to resolve the tiny interval well.

LogisticInterval[itin_]  :=  
      With[{prec = Floor[15 + Length[itin]/1.5]}, 
N[
          Sin[Pi*TentInterval[itin]/2]^2, prec ]]
LogisticInterval[{0, 0}]
LogisticInterval[{0, 1}]
LogisticInterval[{1, 0}]
LogisticInterval[{1, 1}]

Here is an interval for a randomly chosen itinerary 67 bits long.

itin  =  
      Table[Random[Integer], {i, 1, 67}]
LogisticInterval[itin]

Now you can have some fun by picking interesting itineraries and plotting the corresponding orbits of the logistic map. Start off with a period three orbit with itinerary 101101101....  To plot an orbit we will just use the midpoint of the interval as initial point.

itin  =  {1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 
          0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 
          1};
logint  =  LogisticInterval[itin]
y0  =  (logint[[1]] + logint[[2]])/
          2;
logorb  =  
        f^{0, Length[itin] - 1}[y0];
TimeSeriesPlot[logorb];

Now we will produce an orbit which starts out close to this period 3 behavior and later makes a transition to nearly fixed point behavior and finally returns to the period 3 behavior.  It can be found just by designing the right itinerary.  Note the long string of 1's in the middle to make it nearly fixed.

itin  =  {1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 
          0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 
          0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 
          1};
logint  =  LogisticInterval[itin]
y0  =  (logint[[1]] + logint[[2]])/
          2;
logorb  =  
        f^{0, Length[itin] - 1}[y0];
TimeSeriesPlot[logorb];

Try going back to the previous cells and typing in your own itinerary.  Try a "random" one, some periodic ones, anything you like, they all exist !

Cantor Sets

When the parameter r > 4, the logistic map f(x) = r x (1-x) ceases to have attractors.  Most of the orbits escape to minus infinity.  But there is a Cantor set of orbits left behind and on this invariant set the behavior is chaotic. A similar phenomenon occurs for the slope 3 tent map.  Since the latter can be handled with exact arithmetic, it will be used in this section. In this case the invariant set is The Cantor Set, that is , the classic middle-thirds Cantor set.

T1[x_]  :=  
      If[x ≤ 1/2, 3*x, 3 - 3*x]
T = DS[T1]

First we compute a few random orbits and note their tendency to escape.

orb1  =  T^{0, 20}[0.123]
orb2  =  T^{0, 20}[0.678]
TimeSeriesPlot[{orb1, orb2}, 
        LineColor -> {Blue, 
            Green}, 
PlotRange -> {{0, 20}, {-1, 
              2}}];
CobwebPlot[T, {orb1, orb2}, {-1, 2}, 
        LineColor -> {Blue, Green}];

To get orbits which remain bounded, it is necessary to choose the initial point in the Cantor set.  This can be done by specifying a base 3 expansion consisting of 0's and 2's.  We will create some functions to compute exact values of finite and periodic base three numbers given the list of their "tricimals".

MakeTrigit[trilist_]  :=  
      Module[{n = Length[trilist], x = 0, i}, 
For[i = n, 
          i ≥ 1, i--, 
x  =  
            x/3  +  trilist[[i]]]; 
Return[
          x/3];]
MakePeriodicTrigit[
        trilist_]  :=  
With[{q = 
            Length[trilist]}, 
MakeTrigit[
            trilist]*3^q/(3^q - 1)]

MakeTrigit[{0, 2, 2, 0, 2}]
MakePeriodicTrigit[{0, 2, 2, 0, 2}]

Now we can choose any base three expansion involving only 0's and 2's and plot the corresponding orbit.

x0  =  MakeTrigit[{0, 2, 2, 0, 2}];
orb1  =  T^{0, 20}[x0]
x0  =  
        MakePeriodicTrigit[{0, 2, 2, 0, 2}];
orb2  =  T^{0, 20}[x0]
TimeSeriesPlot[{orb1, orb2}, 
        LineColor -> {Blue, 
            Green}, 
PlotRange -> {{0, 20}, {0, 
              1}}];
CobwebPlot[T, {orb1, orb2}, {-0.5, 1.5}, 
        LineColor -> {Blue, Green}];

Here is an orbit which moves around a lot within the Cantor set.  The trigits are chosen at random.

tlist  =  
      Table[2*Random[Integer], {i, 1, 200}]
x0  =  MakeTrigit[tlist];
orb1 =  T^{0, 200}[x0];
TimeSeriesPlot[
        orb1, 
PlotRange -> {{0, 200}, {0, 
              1}}];
CobwebPlot[T, orb1, {-0.5, 1.5}, 
        LineColor -> Blue];

If you put even a single trigit 1, the orbit will escape.

tlist[[198]] = 1;
tlist
x0  =  MakeTrigit[tlist];
orb1 =  T^{0, 200}[x0];
TimeSeriesPlot[
        orb1, 
PlotRange -> {{0, 200}, {0, 
              1}}];
CobwebPlot[T, orb1, {-1, 1.5}, 
        LineColor -> Blue];

Because the Cantor set consists of a dust of points in a line, it is rather hard to visualize.  We will make a bar-code picture of it where each point in the Cantor set is replaced by a line segment.  First we have to compute all the endpoints in the n-th stage of the construction of the set.  We will do this by repeatedly applying the two branches of the inverse of T to the initial endpoints {0,1}.

Tinv1[x_]  :=  x/3;
Tinv2[x_]  :=  
        1 - x/3;

CantorStep[
        prev_]  :=  Map[Tinv1, prev]~Join~Reverse[
          Map[Tinv2, 
            prev]]

CantorSet[n_] :=  Nest[CantorStep, {0, 1}, n]

Here are the first five steps of the construction.

CantorSet[1]
CantorSet[2]
CantorSet[3]
CantorSet[4]
CantorSet[5]

To make an image we will define a function which turns each point of the set into a vertical line segment.

MakeSegment[x_]  :=  
      Line[{{x, 0}, {x, 1}}]
MakeCantorLines[n_]  :=  Map[MakeSegment, CantorSet[n]]

Here are the pictures for the first five steps of the construction.

Show[Graphics[MakeCantorLines[1]], 
        Axes -> True];
Show[Graphics[MakeCantorLines[2]], 
        Axes -> True];
Show[Graphics[MakeCantorLines[3]], 
        Axes -> True];
Show[Graphics[MakeCantorLines[4]], 
        Axes -> True];
Show[Graphics[MakeCantorLines[5]], Axes -> True];

Finally here is a later stage and several magnifications showing the "fractal" structure, that is, the blown up pictures are similar to the whole.

c9  =  
        Show[Graphics[MakeCantorLines[9]], 
          Axes -> True];
Show[c9, 
        PlotRange -> {{0.66, 0.78}, {0, 
              1}}];
Show[c9, PlotRange -> {{0.74, 0.755}, {0, 1}}];

Lab 3 --- Fractals and Iterated Function Systems

In this section we will use iterated function systems to construct some of the famous fractals.  First we have to enter a collection of similarity transformations and collect them together into an IFS object.  Here are the three maps used to make the Sierpinski triangle.

sf1[{x_, y_}]  :=  {x/2, y/2};
sf2[{x_, y_}]  :=  {0.5 + x/2, y/2};
sf3[{x_, y_}]  :=  {0.25 + x/2, 0.25*Sqrt[3] + y/2};
sierpinski  = IFS[ sf1, sf2, sf3]

The program know how to apply some iterate of the IFS to a given initial point.  The result is a list of points obtained by applying all the maps to the given point.  In this case the nth iterate would produce a list of 3^n points.

sierpinski[{0, 0}]
sierpinski^2[{0, 0}]
sierpinski^3[{0, 0}]

You can plot such an iterate to get an approximate picture of the attractor of the iterated function system.

ListPlot[sierpinski^7[{0, 0}], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];

It is interesting to use an initial shape other than a point.  Simple shapes can be represented as a Mathematica Line object.  Here are some examples

point  =  Line[{{0, 0}, {0, 0}}];
segment  =  Line[{{0, 0}, {1, 0}}];
triangle  =  
        Line[{{0, 0}, {1, 0}, {0.5, N[Sqrt[3]/2]}, {0, 0}}];
square  =  
        Line[{{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}}];
\(\(diamond\  = \ 
        Line[{{0.5, 0}, {0, 0.5}, {\(-0.5\), 0}, {0, \(-0.5\)}, {0.5, 
              0}}];\)
\)
Show[Graphics[{triangle, square, diamond}], AspectRatio -> 1, 
        Axes -> True];

Now we can apply the IFS to one of these shapes.

frac  =  sierpinski[triangle]
Show[Graphics[frac], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];

Here is a movie showing the convergence of the iterates to the attractor.

\(Do[frac\  = \ 
        sierpinski^i[triangle]; 
\(Show[Graphics[frac], AspectRatio -> 1, 
          PlotRange -> {{0, 1}, {0, 1}}];\), {i, 1, 7}]\)

To animate them you have to select them by clicking with the mouse on the second blue bracket on the right of the cell.  Once they are selected, go into the Cell menu and choose the item called Animate Selected Graphics.  You should see a movie of your plots showing how the square is transformed.  To stop the movie just click the mouse.  

If your movie goes too fast or too slow you can fix it as follows.  Go to the Preferences item in the Edit menu. Select GraphicsOptions and then  Animation.  You will find an item called AnimationDisplayTime.  Click on the old value, enter a new one and hit Enter (or Return).  Increasing it will slow the movie down.

Next we can demonstrate the fact that the attractor is independent of the initial shape by plotting a similar picture starting from a square or a diamond.

\(\(iterations\ \  = \ 3;\)
\)
frac  =  sierpinski^iterations[square];
\(\(Show[Graphics[frac], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];\)
\)
frac  =  sierpinski^iterations[diamond];
Show[Graphics[frac], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];

Go back and increase the number of iterations to see that, in the end, it all comes out the same.

Now we will have a look at a couple other fractals. Here are some Sierpinski carpets using 5 and 8 squares.  First we define 9 maps and then take subsets of them to construct the IFS.

s5f1[{x_, y_}]  :=  {x/3, y/3};
s5f2[{x_, y_}]  :=  {x/3 + 1/3, y/3};
s5f3[{x_, y_}]  :=  {x/3 + 2/3, y/3};
s5f4[{x_, y_}]  :=  {x/3, y/3 + 1/3};
s5f5[{x_, y_}]  :=  {x/3 + 1/3, y/3 + 1/3};
s5f6[{x_, y_}] :=  {x/3 + 2/3, y/3 + 1/3};
s5f7[{x_, y_}]  :=  {x/3, y/3 + 2/3};
s5f8[{x_, y_}]  :=  {x/3 + 1/3, y/3 + 2/3};
s5f9[{x_, y_}]  :=  {x/3 + 2/3, y/3 + 2/3};
square5  =  IFS[s5f1, s5f3, s5f5, s5f7, s5f9];
square8  =  
        IFS[s5f1, s5f2, s5f3, s5f4, s5f6, s5f7, s5f8, s5f9];

And now the fractals.  With an IFS like square8 which has 8 functions in it, you need to be careful about how many iterations you use.  Even n = 4 iterations gives you 8^4 = 4096 separate squares in your picture.

frac  = square5^4[square];
\(\(Show[Graphics[frac], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];\)
\)
frac  =  square8^4[square];
Show[Graphics[frac], AspectRatio -> 1, 
        PlotRange -> {{0, 1}, {0, 1}}];

Finally, we will plot the Koch curve.  This time the functions involve some rotation.

kf1[{x_, y_}]  :=  {x/3, y/3};
kf4[{x_, y_}]  :=  {x/3 + 2/3, y/3};
\(\(kf2[{x_, y_}]\  := \ {x/6 - Sqrt[3]*y/6 + 1/3, 
			Sqrt[3]*x/6 + 
            y/6};\)\)
\(\(kf3[{x_, y_}]\  := \ {x/6 + Sqrt[3]*y/6 + 
            1/2, 
			\(-Sqrt[3]\)*x/6 + y/6 + 
            1/\((2*Sqrt[3])\)};\)\)
koch  =  IFS[kf1, kf2, kf3, kf4];

And here's the movie:

\(Do[frac\  = 
        koch^i[segment]; 
\(Show[Graphics[frac], AspectRatio -> 1, 
          PlotRange -> {{0, 1}, {\(-1\)/2, 1/2}}];\), {i, 1, 5}]\)

You can animate as before.

Other

In this lab we will look at some dynamical phenomena peculiar to two dimensions and also plot bifurcations for one and two dimensional systems.

Orbits of Linear Maps

We will plot some orbits of simple linear maps of the plane. In Mathematica, vectors are represented as lists of numbers enclosed in set brackets.  For example, {1,2} is a vector in the plane.  A matrix is given as a list of vectors which constitute the rows of the matrix.  The most general possible 2x2 matrix would be {{a,b},{c,d}}.  The built-in Mathematica function MatrixForm just prints the matrix in a nice way.

\(A\  := \ {{a, b}, {c, d}}; 
MatrixForm[A]\)

There are many Mathematica functions to analyze matrices.  After setting specific values for a,b,c,d, a few useful commands are demonstated below.  The chosen values represent the Fibonacci matrix.

\(a\  = \ 1; \ b = 1; 
c = 1; \ d = 0; 
"A =" MatrixForm[A]\)"inverse =" MatrixForm[Inverse[A]]Eigenvalues[A]Eigenvectors[A]

To view a matrix as a dynamical system, we define a linear function which takes as input a vector {x,y} and produces as output, the vector obtained by multiplying by the matrix A.  Then we define a dynamical system object using this function.  As before, this allows the program to recognize the commands for iteration and orbit computation.

l1[{x_, y_}]  :=  A . {x, y}\(L\  = \ DS[l1]
\)L[{x, y}]L[{1, 5}]

Just as for one-dimensional maps, you can use the standard exponent notation to find an iterate of this map applied to an initial vector x0.   To compute an orbit just put a range of values in the exponent. The result is then an Orbit object containing a list of vectors which are the individual points of the orbit.  In this case, the coordinates of these points are the Fibonacci numbers.

\(x0\  = \ {1, 0}; 
L^100[x0]\)orb  =  L^{0, 50}[x0]

Next we plot the orbit using a vector plotting function called PhasePlot.  Since the origin is a saddle point, the orbit just heads to infinity in the direction of the unstable subspace.

PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True];

To make an interesting plot of a saddle point, you need to start close to the stable subspace.  The eigenvectors were computed above and the stable eigenvector is {(1-Sqrt[5])/2,1}.  The point x0 below is approximately, but not exactly, equal to this.  The orbit hops back and forth near the stable subspace (it's a flip saddle) and then takes off along the unstable direction.

\(x0\  = \ {\(-6.18\), 10}; 
orb\  = \ L^{0, 25}[x0]\)PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True];

It is more interesting to look at a case where the origin is an attractor. Here is a spiral sink, i.e., the eigenvalues are complex with norm less than 1.  By changing the numbers  r and theta below you can produce some really nice pictures.

\(
r\  = \ 0.95; \ theta\  = 2*Pi/7; 
a\  = \ r*Cos[theta]; 
    b\  = \ \(-r\)*Sin[theta]; 
c\  = \ r*Sin[theta]; d\  = \ r*Cos[theta]; 
    
Eigenvalues[A]\)Abs[Eigenvalues[A]]\(x0\  = \ {1, 0}; 
orb\  = \ L^{0, 50}[x0]; 

    PhasePlot[orb, PlotRange -> All, PointSize -> 0.02, Lines -> True]; \)

Here is a sink with real eigenvalues.   Several different orbits are shown to  illustrate how they all approach the origin tangent to the weaker eigenspace (in this case, the x-axis).  Try switching the values of a and d so that the other eigenvector becomes the weaker one.

\(a = 0.8; b = 0; 
c = 0; d = 0.7; 
Eigenvalues[A]\)\(x1\  = \ {1, \(-1\)}; orb1 = L^{0, 20}[x1]; 
x2\  = \ {1, \(-0.1\)}; 
    orb2 = L^{0, 20}[x2]; 
x3\  = \ {1, 1}; orb3 = L^{0, 20}[x3]; 

    x4\  = \ {\(-1\), 1}; orb4 = L^{0, 20}[x4]; 
x5\  = \ {\(-1\), 0.1}; 
    orb5 = L^{0, 20}[x5]; 
x6\  = \ {\(-1\), \(-1\)}; orb6 = L^{0, 20}[x6]; 
    

PhasePlot[{orb1, orb2, orb3, orb4, orb5, orb6}, 
      PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}, PointSize -> 0.02, 
      Lines -> True]; 
\)

Orbits the Henon Map

We begin by entering the formulas for the Henon map and declaring it to be a dynamical system object.    A mapping of the plane such as the Henon map takes a vector as input and produces another vector as output.  It can be entered as shown below.

h1[{x_, y_}]  := {a - x^2 + b*y, x}henon  =  DS[h1]

Next we set the values of the parameters, define  an initial point, compute an orbit and plot it. Since we will be plotting long orbits we will plot them with small dots not connected by lines. In this case we get the Henon attractor.

\(a = 1.4; \ b\  = \ 0.3; 
x0\  = \ {1, 1}; 

    orb\  = \ henon^{0, 5000}[x0]; 

    PhasePlot[orb, PlotRange -> {{\(-1.6\), 2}, {\(-1.6\), 2}}]; \)

To help understand this attractor we will  plot a sequence of pictures where the value of the parameter a is gradually increased.  At each stage, we plot the omega limit set of the orbit with initial conditions {1,1}.  This is done by skipping the first 100 points along the orbit to allow time to converge to whatever attractor may be present.
	To make the simpler attractors visible, we can use larger dots in the plot.  Starting from the value a = 0.3 for which there is an attracting fixed point, we progress through a sequence of period-doubling bifurcations arriving at a strange attractor with several parts.  There are interludes where a periodic attractor returns eventually leading to the  attractor found above.  Try stepping through the a values suggested below.

\(a = 0.3; 
    \ \  (*\ try\ a\  = \ 0.4, \ 0.9, \ 0.95, \ 1.05, \ 1.055, \ 1.07, \ 1.1, 
      \ 1.25, \ 1.27, \ 1.3, \ 1.4*) \ \ 
b\  = \ 0.3; 
x0\  = \ {1, 1}; 

    orb\  = \ henon^{100, 2000}[x0]; 

    PhasePlot[orb, PlotRange -> {{\(-1.6\), 2}, {\(-1.6\), 2}}, 
      PointSize -> 0.02]; \)

For a close-up look at the attractor, you can plot a very long orbit and plot blow-ups of it by changing the PlotRange option.

\(a = 1.4; \ b\  = \ 0.3; 
x0\  = \ {1, 1}; 

    orb\  = \ henon^{0, 20000}[x0]; 

    PhasePlot[orb, PlotRange -> {{0.3, 0.9}, {0.6, 1.2}}]; 

    PhasePlot[orb, PlotRange -> {{0.4, 0.6}, {0.9, 1.1}}]; \)

Bifurcations and Orbit Diagrams

In this section we will have a look at the amazing details of the orbit diagram for the logistic map and also construct a similar diagram for the two-dimensional Henon map.  First we define the logistic dynamical system.  The bifurcation plotting software is set up to vary a parameter called c so we will use that letter in the formula.

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

To construct an orbit diagram we need to start with some initial condition.  In general this could vary with c so we define an initial condition function ic[c] which gives the desired starting point.  In this case it is a constant function at the critical point x0 = 0.5.

ic[c_]  :=  0.5

Now the function bifurcation plot can be used to draw the orbit diagram. Here is some information about this function, whose actual definition is contained in the initialization section of this notebook.

?BifurcationPlot

BifurcationPlot[f, ic, {0, 4}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]

Of course the most interesting part of the diagram is from c=2 to c=4.  You can zoom in on that part by changing the range for c.

BifurcationPlot[f, ic, {2, 4}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]

It is endlessly amusing to blow up small regions of the diagram and to ponder the incredible complexity contained in it.  For example here is a tiny window where a period 5 orbit is born in a saddle-node bifurcation and then proceeds to period double producing attractors of period 10, 20 , 40 , etc.

BifurcationPlot[f, ic, {3.73, 3.75}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, PlotRange -> {0, 1}]

You can blow up just one of the five period-doubling regions by changing the   PlotRange as well as the range for c.  Here a blow-up of the part of the last diagram near the center revealing further smaller windows.  For deep blow-ups it is sometimes good to increase the Delay to give more time to converge to the attractor and the number of OrbitPoints, since not all of them will appear in the chosen  Plot Range.

BifurcationPlot[f, ic, {3.742, 3.745}, ParameterPoints -> 300, 
      OrbitPoints -> 400, Delay -> 200, PlotRange -> {0.45, 0.55}]

I suggest that you go back to the previous cell and  explore other regions of the orbit diagram by making changes to the c range and the Plot Range.

Next we will make a similar plot for a one-parameter family of Henon maps.  Here is the family which was studied using simple forward iteration in the last section.  Once again, we have called the parameter c.

h1[{x_, y_}]  :=  {c + 0.2*y  - x^2, x}h  =  DS[h1]

Now the orbits for the Henon map consist of vectors.  To get a 2D plot we have to plot c versus one or the other of the x, y coordinates (or some combination of them).  Here is a  function which projects a 2D vector onto its x coordinate and another which takes a vector orbit and produces its list of x coordinates

xproj[{x_, y_}]  :=  xXProj[orb_Orbit]  :=  Map[xproj, orb[]]

This can be supplied as an option to the BifurcationPlot function to construct an orbit diagram for the Henon family.  We also need to supply an initial conction somewhere in the basin of the attractor. The point {1,1} will do.

ic[c_]  :=  {1.0, 1.0}

BifurcationPlot[h, ic, {0, 1.6}, ParameterPoints -> 300, 
      OrbitPoints -> 200, Delay -> 100, Transformation -> XProj]

By blowing-up parts of this picture, you can find subattractors of the Henon attractor with various periods.

Stable and Unstable Manifolds

Next we will plot some stable and unstable manifolds for  2D diffeomorphisms.   The first step is to find a saddle fixed or periodic point for the map.  Of course we begin with the Henon map.

h1[{x_, y_}]  := {a - x^2 + b*y, x}henon  =  DS[h1]

After setting some familiar values for the parameters we search for a fixed point using the function SeekPeriodicOrbit.  This function looks for a solution of the fixed point equations near a given initial point.  If it finds one it returns an Orbit object containing the point and its first iterate (so you can see it really is fixed).

\(a = 0.8; 
b\  = \ 0.2; 

    orb\  = \ SeekPeriodicOrbit[henon, 1, {1, 1}]\)

You can test to see if it is a saddle point by computing its multipliers.

Multiplier[henon, orb]

Since one of them has absolute value bigger than one and one has absolute value less than one, it is indeed a saddle point.  The first entry of the orbit is the actual fixed point which we will call X0.

X0  =  orb[0]

Now we can use the magic function UnstableManifold to compute a piece of the unstable manifold near the fixed point.

?UnstableManifold

Wu  =  UnstableManifold[henon, X0, 1]

The result is a list of 40 points appoximately on  the unstable manifold of the saddle,  20 on each side of the fixed point.  We can plot this to see what is happening.

PhasePlot[Wu, Lines -> True, Points -> False, LineColor -> Blue, 
      PlotRange -> {{-1, 2}, {-1, 2}}];

This is only a small piece of the unstable manifold.  To extend it just change the third argument of the UnstableManifold function to increase the number of iterations it uses to extend the manifolds.

\(Wu\  = \ UnstableManifold[henon, X0, 3]; 

    PhasePlot[Wu, Lines -> True, Points -> False, LineColor -> Blue, 

      PlotRange -> {{\(-1\), 2}, {\(-1\), 2}}]; \)

Go back and gradually increase the number of iterations used up to about 10 to see the manifold grow.

To get a stable manifold, we just compute the unstable manifold for the inverse of the map.  So first we have to input the formulas for the inverse of the Henon map.

hinv1[{x_, y_}]  :=  {y, (x + y^2 - a)/b}henoninv  =  DS[hinv1]

The next command will compute pieces of both the stable and unstable manifolds and plot them on the same plot.  The unstable manifold is blue and the stable manifold is green.

\(a\  = 0.8; 
b = \ 0.2; 

    orb\  = \ SeekPeriodicOrbit[henon, 1, {1, 1}]; 
X0\  = \ orb[0]; 

    Wu\  = \ UnstableManifold[henon, X0, 6, Subdivisions -> 200]; 

    Ws\  = \ UnstableManifold[henoninv, X0, 4]; 

    PhasePlot[{Ws, Wu}, Lines -> True, Points -> False, 
      LineColor -> {Green, Blue}, 
PlotRange -> {{\(-2\), 2}, {\(-2\), 2}}]; 
    \)

As you can see, the stable manifold grows a lot faster (under backward iteration) than the unstable manifold does.  Now go back and gradually increase the value of the parameter a over the range [0.8,1.6] to watch how the manifolds change.  If you get ragged picture try increasing the subdivisions used to compute the unstable manifold or decreasing the number of iterations used to extend it.

Now the classic Henon attractor occurs for the parameters given below.  Here is a picture showing an orbit on the attractor as well as the stable and unstable manifolds of the saddle fixed point inside the attractor.

\(a\  = \ 1.4; 
b = \ 0.3; 

    orb\  = \ SeekPeriodicOrbit[henon, 1, {1, 1}]; 
X0\  = \ orb[0]; 

    Wu\  = \ UnstableManifold[henon, X0, 6, Subdivisions -> 100]; 

    Ws\  = \ UnstableManifold[henoninv, X0, 4]; 

    orb\  = \ henon^{100, 1000}[{1.0, 1.0}]; 

    PhasePlot[{Ws, Wu, orb}, Lines -> {True, True, False}, 
      Points -> {False, False, True}, LineColor -> {Green, Blue, Black}, 
      PointSize -> 0.01, 
PlotRange -> {{\(-2\), 2}, {\(-2\), 2}}]; \)

As you can see, the attractor is virtually indistinguishable from the unstable manifold of the saddle point.

Another interesting 2D diffeomorphism is the so-called standard mapping of the annulus.  It is a kind of discrete dynamical system analogue of the forced, undamped pendulum differential equation.  The x-variable is an angle (analogous to the position of the pendulum) so the state space is really a cylinder or annulus instead of a plane.  First we enter the formulas for the map and its inverse.  The Mod command in the first component forces the angular coordinate x to always lie in the interval from 0 to 2*Pi

st1[{x_, y_}]  :=  {Mod[x + y + a*Sin[x], 2*Pi], y + a*Sin[x]}stmap  =  DS[st1]stinv1[{x_, y_}]  :=  {Mod[x - y, 2*Pi], y - a*Sin[x - y]}stinv  =  DS[stinv1]

First we just plot some orbits.

\(a\  = \ 0.1; 
orb1\  = \ stmap^{0, 1000}[{0.1, 0}]; 

    orb2\  = \ stmap^{0, 1000}[{0, 0.1}]; 

    orb3\  = \ stmap^{0, 1000}[{1, 0}]; 
orb4\  = \ stmap^{0, 1000}[{2, 0}]; 
    
orb5\  = \ stmap^{0, 1000}[{0, 1}]; 

    orb6\  = \ stmap^{0, 1000}[{0, 2}]; 

    orb7\  = \ stmap^{0, 1000}[{0, \(-1\)}]; 

    orb8\  = \ stmap^{0, 1000}[{0, \(-2\)}]; 

    PhasePlot[{orb1, orb2, orb3, orb4, orb5, orb6, orb7, orb8}, 
      PointSize -> 0.002]; \)

As you can see the plot resembles the familiar  picture for the pendulum.  Now go back to the last cell and try some other values for  the parameter a.  For example, try  a = 0.5. You will nottice a small chaotic band associated with the orbits which started out near the origin.  Later we will see that this is an example of homoclinic chaos produced by a transversal homoclinic point to the saddle fixed point at the origin.  Try some larger values of a such as a = 0.8,  1.0, 1.2, 1.5.  You will find that the chaotic zone gradually takes over thestate space.

To understand the origin of the chaos we will plot the stable and unstable manifolds of the fixed point at the origin.  It turns out that the best pictures are obtained by thinking of the points
(0,0), (2*Pi,0), (4*Pi,0) etc as distinct.  For this reason we will first redefine the dynamical system so that it does not force the x-coordinate to lie in the interval [0, 2*Pi].

st2[{x_, y_}]  :=  {x + y + a*Sin[x], y + a*Sin[x]}stmap2  =  DS[st2]sti2[{x_, y_}]  :=  {x - y, y - a*Sin[x - y]}stinv2  =  DS[sti2]

First we check that the origin is a saddle point.

\(a = 0.8; 
orb\  = \ SeekPeriodicOrbit[stmap2, 1, {0, 0}]\)Multiplier[stmap2, orb]

Now we will plot the stable and unstable manifolds of the  origin and also those of its translate at (2*Pi,0).  You can easily see the transverse intersections of the stable and unstable manifolds and also how these manifolds accumulate on themselves near the fixed point.

\(a = 0.8; 
iters\  = \ 3; 

    wu\  = \ UnstableManifold[stmap2, {0, 0}, 7, Radius -> 0.5]; 

    ws\  = \ UnstableManifold[stinv2, {0, 0}, 7, Radius -> 0.5]; 

    wu1\  = \ UnstableManifold[stmap2, {2*Pi, 0}, 7, Radius -> 0.5]; 

    ws1\  = \ UnstableManifold[stinv2, {2*Pi, 0}, 7, Radius -> 0.5]; 


    PhasePlot[{wu, ws, wu1, ws1}, Lines -> True, \ Points -> False, 

      LineColor -> {Blue, Green, Blue, Green}, 
      PlotRange -> {{\(-2\), 8}, {\(-3\), 3}}]; 


    PhasePlot[{wu, ws, wu1, ws1}, Lines -> True, \ Points -> False, 

      LineColor -> {Blue, Green, Blue, Green}, 
      PlotRange -> {{\(-1\), 1}, {\(-1\), 1}}]; \)

Go back and try increasing a to a = 1.0 and a = 1.2.

Delay Embedding and Attractor Reconstruction

In this section we will compute some orbits of vector dynamical systems, and then try to reconstruct the dynamics using a scalar "measurement function".  For example here is a simple rotation map of the plane producing an orbit which is dense in a circle.

f1[{x_, y_}]  := 
       {Cos[2*Pi*omega]*x - Sin[2*Pi*omega]*y, 
        Sin[2*Pi*omega]*x + Cos[2*Pi*omega]*y}f  =  DS[f1];

Here is a rotation by one radian, an irrational fraction of a whole rotation.

\(omega\  = \ 1/\((2*Pi)\); 
orb\  = \ f^{0, 1000}[{1.0, 0.0}]; 

    PhasePlot[orb]; \)

Next we construct a scalar time series by applying some scalar measurement function to each point of this orbit.  We will try the measurement function h = x+y^2

\(h[{x_, y_}]\  := \ x + y^2; 
timeseries\  = \ Orbit[Map[h, orb[]]]; 
    \)

Now that we have a scalar time series we replot it as a  series of delay vectors and compare it with the original plot.

DelayPlot[timeseries, 2, PointSize -> 0.004];

This is not homeomorphic to the original circle.  Moving up to three dimensions solves the problem however.

DelayPlot[timeseries, 3, PointSize -> 0.004];

Now for a more interesting orbit.  Here is an orbit from the Henon attractor.

h1[{x_, y_}]  := {a - x^2 + b*y, x}\(henon\  = \ DS[h1]; 
a = 1.4; b = 0.3; 

    orb\  = \ henon^{100, 10000}[{1, 1}]; 

    PhasePlot[orb, PointSize -> 0.002]; \)

As a measurement function we will use  x+y

h[{x_, y_}]  :=  x + y\(timeseries\  = \ Orbit[Map[h, orb[]]]; 

    DelayPlot[timeseries, 2, PointSize -> 0.002]; \)

Again going to 3D eliminates the self-intersections producing an object homeomorphic to the original attractor.

DelayPlot[timeseries, 3, PointSize -> 0.002];

Finally we consider a higher dimensional example.  You can construct a 4D system analogous to the rotation of the circle, which involves two separate rotations.  The resulting orbit will generally be dense in the surface of a torus (the product space of two circles).

\(f1[{x1_, y1_, x2_, y2_}]\  := 
      \ {Cos[2*Pi*omega1]*x1 - Sin[2*Pi*omega1]*y1, 
        Sin[2*Pi*omega1]*x1 + Cos[2*Pi*omega1]*y1, 
		
        Cos[2*Pi*omega2]*x2 - Sin[2*Pi*omega2]*y2, 
        Sin[2*Pi*omega2]*x2 + Cos[2*Pi*omega2]*y2}\)f  =  DS[f1];

\(omega1\  = \ 1/\((2*Pi)\); 
omega2\  = \ 1.618; 

    orb\  = \ f^{0, 5000}[{1.0, 0.0, 1.0, 0.0}]; 
orb[{0, 20}]\)

The orbit is a list of 4D points so we can't plot it directly at all.  We will use a measurement function to reduce it to a scalar time series.

h[{x1_, y1_, x2_, y2_}]  :=  x1 + x2timeseries  =  Orbit[Map[h, orb[]]];

Here are the delay plots in 2 and in 3 dimensions.  You can almost see the torus but it seems that even in 3D there are still self-intersections.  According to the theory we should really be using an embedding dimension of 5 !

\(DelayPlot[timeseries, 2, PointSize -> 0.002]; 

    DelayPlot[timeseries, 3, PointSize -> 0.002]; \)

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]