Dynamical Systems and Chaos

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  Evaluation submenu choose the command Evaluate Initialization Cells.  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::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}]

Tests for different kinds of orbits

ScalarOrbitQ[orb_] := (Head[orb]===Orbit)&&(Length[orb[]]>=1)&&ScalarQ[orb[0]]//LogicalExpand

VectorOrbitQ[orb_] := (Head[orb]===Orbit)&&(Length[orb[]]>=1)&&VectorQ[orb[0]]//LogicalExpand

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[0]]//LogicalExpand

VectorOrbitQ[orb_] := (Head[orb]===Orbit)&&(Length[orb[]]>=1)&&VectorQ[orb[0]]//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[{}]."

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

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

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

MakePoints[l_Orbit,color_:{0,0,0},size_:0.01] :=MakePoints[l[[1]],color,size]
MakeLine[l_Orbit,color_:{0,0,0},size_:0.01] :=MakeLine[l[[1]],color,size]
MakeArrows[l_Orbit,color_:{0,0,0},size_:0.01] :=MakeArrows[l[[1]],color,size]




MakePoints[l_?TwoDListQ,color_:{0,0,0},size_:0.01] :=
	Graphics[{Apply[RGBColor,color],PointSize[size],Map[Point,l]}]

MakePoints[l_?ComplexListQ,color_:{0,0,0},size_:0.01] :=
	Graphics[{Apply[RGBColor,color],PointSize[size],Map[Point[{Re[#],Im[#]}]&,l]}]

MakeLine[l_?TwoDListQ,color_:{0,0,0},thickness_:0.005] :=
	Graphics[{Apply[RGBColor,color],Thickness[thickness],Line[l]}]

MakeArrows[l_?TwoDListQ,color_:{0,0,0},thickness_:0.005] :=
	With[{ll={Drop[l,-1],Drop[l,1]}//Transpose},Graphics[{Apply[RGBColor,color],Thickness[thickness],Map[Arrow,ll]}]]

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

MakeLine3D[l_?ThreeDListQ,color_:{0,0,0},thickness_:0.005] :=
	Graphics3D[{Apply[RGBColor,color],Thickness[thickness],Line[l]}]

MakePolygon3D[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[2Pi/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

Clear[GenericPlot]

IdentityTransformation[orb_] := orb

Options[GenericPlot] = {
		Points->True,
		Lines->True,
	  Arrows->False,
		PointColor->Black,
		LineColor->Black,
		PointSize->0.01,
		LineSize->0.002,
		Transformation->IdentityTransformation,
		MoreGraphics->{}
		};
genoptpattern = (Points->_)|(Lines->_)|(Arrows->_)|(PointColor->_)|
		(LineColor->_)|(PointSize->_)|(LineSize->_)|
		(Transformation->_)|(MoreGraphics->_);

GenericPlot[orb_,opts___] := 
		Module[{showopts,pf,lf,af,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];
		    af = Arrows/.{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,ImagePadding->Automatic};
			trorb = tr[orb];
			gr = {
				If[lf,MakeLine[trorb,lc,ls],{}],
			If[af,MakeArrows[trorb,lc,ls],{}],
				If[pf,MakePoints[trorb,pc,ps],{}]
				}~Join~mg;
			Show[gr,showopts]
			]

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,af,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];
		     af = Arrows/.{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[Head[af]=!=List,af=Table[af,{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,ImagePadding->Automatic};
			trorbs = Map[tr,orbs];
			gr = Flatten[{
				MapThread[If[#2,MakeLine[#1,#3,#4],{}]&,
								{trorbs,lf,lc,ls}],
		     MapThread[If[#2,MakeArrows[#1,#3,#4],{}]&,
								{trorbs,af,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,GenericPlot3D[orb,opts,Lines->False,
			Transformation->(EmbedOrbit[#,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[2Pi]},
						Map[{Cos[tp*#],Sin[tp*#]}&,orb]]

Plots for Vector and Complex 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[]]

ComplexTransformation[orb_] := Map[{Re[#],Im[#]}&,orb[]]

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

ComplexPlot[orb_,opts___] :=
	GenericPlot[orb,opts,Transformation->ComplexTransformation,
				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->False,
	RoundDigits->6,
	PlotRange->All,
	Transformation->IdentityTransformation};
	
bpoptpattern = (ParameterPoints->_)|(OrbitPoints->_)|
	(Delay->_)|(ClearOut->_)|(RoundDigits->_)|(Transformation->_);
	
(*BifurcationPlot[f_DS,ic_,pr_,opts___] :=
	Module[{showopts,pp,op,de,co,rd,plr,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];
		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]]];
		Print[orb,Map[Rnd,orb]];
			orb = Clp[Union[Rnd[orb[[1]]]]];
			ptlist = ptlist~Join~
				Transpose[{Table[c,{Length[orb]}],orb}];
			];
		(*Print[Length[ptlist]];*)
		Show[Graphics[{PointSize[0.001],
				Map[Point,ptlist]}],showopts]
		]*)

BifurcationPlot[f_DS,ic_,pr_,opts___] :=
	Module[{showopts,pp,op,de,co,rd,plr,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];
		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]]];
		(*Print[Union[Map[Rnd,orb]]];*)
			orb = Clp[Union[Map[Rnd,orb[[1]]]]];
			ptlist = ptlist~Join~
				Transpose[{Table[c,{Length[orb]}],orb}];
			];
		(*Print[Length[ptlist]];*)
		Show[Graphics[{PointSize[0.001],
				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:{__Complex}] := Flatten[Map[f,x],1]
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]]]

f_IFS[x_Complex] := Map[#[x]&,SwitchHead[f,List]]
f_IFS[x_Real] := Map[#[x]&,SwitchHead[f,List]]
f_IFS[x_?NumberQ] := Map[#[x]&,SwitchHead[f,List]]


IFS/: f_IFS^n_Integer[x_] := Nest[f,x,n]

IFS/: f_IFS^n_Integer[x:{__Line}] := Nest[f,x,n]

Julia Sets and the Mandelbrot Set

(* Compute escape times for finding Julia sets *)
Clear[EscapeTime]
EscapeTime[f_DS,z0_,er_,nmax_] :=
	Module[{n=0,z=z0},
		While[Abs[z]<=er && n<nmax,
			z = f[z];
			n++;
			];
		n
		]
EscapeTime::usage = "EscapeTime[f,z0,er,nmax] finds the number of 
iterations of the complex function f required for the orbit of z0
to leave the ball of escape radius, er.  It returns nmax if the orbit
does not escape after that many iterations." 

(* This constructs a compiled function for finding escape times
for the given function f[z] and for a given escape radius and 
maximum number of iterations *)

MakeEscapeTime[f_,er_,nmax_] := Module[{g,compilecommand,z,n},
	compilecommand = Hold[
		Compile[{{z0,_Complex}},
			Block[{n=0,z=z0},
			While[Abs[z]<=er && n<nmax,
				z=g;
				n++;
				];
			n
			],{{z,_Complex},{n,_Integer}}]]/.g->f[z];
	ReleaseHold[compilecommand]
	]

(* And one for finding Mandelbrot sets.  Here we have a family, f[z,c]
and a function returning the location of the critical point, cp[c].*)

CriticalEscapeTime[f_,cp_,c0_,er_,nmax_] :=
	Module[{n=0,z=cp[c0]},
		While[Abs[z]<=er && n<nmax,
			z = f[z,c0];
			n++;
			];
		n
		]
CriticalEscapeTime::usage = "CriticalEscapeTime[f,cp,c0,er,nmax]
finds the number of interations of the family of complex functions
f[z,c0] required for a critical orbit to leave the ball of escape
radius, er.  Here cp[c] is a function giving the location of the 
critical points. It returns nmax if the critical orbit has not
escaped after that many iterations."

(* And here is a compiled function for Mandelbrot escape times *)

MakeCriticalEscapeTime[f_,cp_,er_,nmax_] := 
Module[{g,h,compilecommand,z,n},
	compilecommand = Hold[
		Compile[{{c0,_Complex}},
			Block[{n=0,z=h},
			While[Abs[z]<=er && n<nmax,
				z=g;
				n++;
				];
			n
			],{{z,_Complex},{n,_Integer}}]]/.{g->f[z,c0],h->cp[c0]};
	ReleaseHold[compilecommand]
	]

(* Some standard colors as a table of substitution rules *)

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 coloring functions for filled Julia sets and the Mandelbrot set.
	These are used for the bands around the sets.  The color of the sets
	themselves is determined by the option MaxIterationColor.*)

(* Sixteen shades of gray 0=white, 15=dark gray *)
GrayArray = Table[RGBColor[1-0.06*i,1-0.06*i,1-0.06*i],{i,0,15}];

GrayArray = Table[List[1-0.06*i,1-0.06*i,1-0.06*i],{i,0,15}];
Grays[et_] := GrayArray[[1+Floor[Mod[et,16]] ]]

(* 10 blues *)
Blues[et_] :=
{RGBColor[0.8,1,1],RGBColor[0.6,1,1],RGBColor[0.4,1,1],
RGBColor[0.2,1,1],RGBColor[0,0.8,1],RGBColor[0,0.6,1],
RGBColor[0,0.4,1],RGBColor[0,0.2,1],RGBColor[0,0,1],
RGBColor[0,0,0.5]}[[ 1+Floor[Mod[et,10]] ]] 

Blues[et_] :=
{List[0.8,1,1],List[0.6,1,1],List[0.4,1,1],
List[0.2,1,1],List[0,0.8,1],List[0,0.6,1],
List[0,0.4,1],List[0,0.2,1],List[0,0,1],
List[0,0,0.5]}[[ 1+Floor[Mod[et,10]] ]] 

(* 10 warm colors *)
(*Warms[et_] := With[{g=N[(16-Mod[et,16])/16]},RGBColor[g,g,0]]*)

Warms[et_] :=
{RGBColor[1,1,0.6],RGBColor[1,1,0.4],
RGBColor[1,1,0.2],RGBColor[1,1,0],RGBColor[1,0.8,0],
RGBColor[1,0.6,0],RGBColor[1,0.4,0],RGBColor[1,0.2,0],
RGBColor[1,0,0],RGBColor[0.5,0,0]}[[ 1+Floor[Mod[et,10]] ]] 

Warms[et_] :=
{List[1,1,0.6],List[1,1,0.4],
List[1,1,0.2],List[1,1,0],List[1,0.8,0],
List[1,0.6,0],List[1,0.4,0],List[1,0.2,0],
List[1,0,0],List[0.5,0,0]}[[ 1+Floor[Mod[et,10]] ]] 


(* White  *)
BlackAndWhite[et_] := List[1,1,1]

(* Gray and White Stripes *)
Stripes[et_] := If[EvenQ[Round[et]],List[1.0,1.0,1.0],
						List[0.5,0.5,0.5]]

Clear[FilledJuliaSet]
Options[FilledJuliaSet] = {
	EscapeRadius->2.0,
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Grays,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->100,
	PlotRange->{{-2.0,2.0},{-2.0,2.0}}
	};
fjsoptpattern = (EscapeRadius->_)|(MaxIterations->_)|
	(MaxIterationColor->_)|(ColorFunc->_)|(PlotPoints->_)|
	(ClearOut->_)|(ColorFunction->_)|(BandWidth->_);

FilledJuliaSet[f_DS,opts___] :=	Module[
	{showopts,er,mi,mic,cf,cf2,cp,bw,co,pp,pr,etfunc,escapetimes,dx,dy,critorb},
	er = EscapeRadius/.{opts}/.Options[FilledJuliaSet];
	mi = MaxIterations/.{opts}/.Options[FilledJuliaSet];
	mic = MaxIterationColor/.{opts}/.Options[FilledJuliaSet]/.ColorTable;
	cf2 = ColorFunction/.{opts}/.Options[FilledJuliaSet];
	bw = BandWidth/.{opts}/.Options[FilledJuliaSet];
	co = ClearOut/.{opts}/.Options[FilledJuliaSet];
	pp = PlotPoints/.{opts}/.Options[FilledJuliaSet];
	pr = PlotRange/.{opts}/.Options[FilledJuliaSet];
	dx = N[(pr[[1,2]]-pr[[1,1]])/pp];
	dy = N[(pr[[2,2]]-pr[[2,1]])/pp];
	
	showopts = DeleteCases[{opts},fjsoptpattern]~Join~
			{Frame->True,AspectRatio->dy/dx};
	cf = Function[t,If[t>=mi,List@@mic,cf2[t/bw]]];					
	If[co,Unprotect[Out];Clear[Out];Protect[Out]];
(*	escapetimes = Table[
		EscapeTime[f,pr[[1,1]]+j*dx+I*(pr[[2,1]]+i*dy),er,mi],{i,0,pp},{j,0,pp}];  *)
		
		(* Compile an escape time function *)					
	etfunc = MakeEscapeTime[f[[1]],er,mi];
	escapetimes = Table[
		etfunc[pr[[1,1]]+j*dx+I*(pr[[2,1]]+i*dy)],{i,0,pp},{j,0,pp}];	
		
	
	
	Show[Graphics[
		Raster[Map[cf,escapetimes,{2}],Transpose[pr]]
		],showopts]
	]
	
FilledJuliaSet::usage = "FilledJuliaSet[f,opts] plots the filled
Julia set of a complex polynomial using the escape time algorithm.
There are many options to control the computation.  For example
if f[z] := z^2+c and the value of c has been set, then
FilledJuliaSet[f,
	EscapeRadius->2.0,
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Grays,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->100,
	PlotRange->{{-2.0,2.0},{-2.0,2.0}}]
will make a plot of escape times of a 100 x 100 grid of initial points
in the rectangle [-2,2]x[-2,2].  The points which do not escape in 
50 iterations will be colored black and the ones which do escape will
be colored using a function called Grays which selects a shade of
gray (other built in color function are Blues, Warms, Stripes,
BlackAndWhite). Previous outputs will be cleared to free memory. The 
BandWidth option can be increased to color groups of escape times with
the same color."

FilledJuliaSet[f_,opts___] :=	Module[
	{showopts,er,mi,mic,cf,cf2,cp,bw,co,pp,pr,etfunc,escapetimes,dx,dy,critorb},
	er = EscapeRadius/.{opts}/.Options[FilledJuliaSet];
	mi = MaxIterations/.{opts}/.Options[FilledJuliaSet];
	mic = MaxIterationColor/.{opts}/.Options[FilledJuliaSet]/.ColorTable;
	cf2 = ColorFunction/.{opts}/.Options[FilledJuliaSet];
	bw = BandWidth/.{opts}/.Options[FilledJuliaSet];
	co = ClearOut/.{opts}/.Options[FilledJuliaSet];
	pp = PlotPoints/.{opts}/.Options[FilledJuliaSet];
	pr = PlotRange/.{opts}/.Options[FilledJuliaSet];
	dx = N[(pr[[1,2]]-pr[[1,1]])/pp];
	dy = N[(pr[[2,2]]-pr[[2,1]])/pp];
	
	showopts = DeleteCases[{opts},fjsoptpattern]~Join~
			{Frame->True,AspectRatio->dy/dx};
	cf = Function[t,If[t>=mi,List@@mic,cf2[t/bw]]];					
	If[co,Unprotect[Out];Clear[Out];Protect[Out]];
	(*escapetimes = Table[
		EscapeTime[f,pr[[1,1]]+j*dx+I*(pr[[2,1]]+i*dy),er,mi],{i,0,pp},{j,0,pp}];*)
	(* Compile an escape time function *)					
	etfunc = MakeEscapeTime[f,er,mi];
	escapetimes = Table[
		etfunc[pr[[1,1]]+j*dx+I*(pr[[2,1]]+i*dy)],{i,0,pp},{j,0,pp}];	
	
		
	Show[Graphics[
		Raster[Map[cf,escapetimes,{2}],Transpose[pr]]
		],showopts]
	]

(* This function computes the Julia set by the inverse iteration 
algorithm. It returns a list of points which can be plotted 
together with other dynamical features using ComplexPlot. It follows a
random branch through the preimage tree with probabilities specified
by the Probabilities option. *)

Options[JuliaSet] = {
	InitialPoint->0.0+0.0*I,
	PlotPoints->1000,
	Probabilities->Automatic,
	TreeDepth->5,
	RandomDepth->5
	};
	
jsoptpattern = (InitialPoint->_)|(PlotPoints->_)|
	(Probabilities->_)|(TreeDepth->_)|(RandomDepth->_);

JuliaSet[fs_,opts___] :=	Module[
	{showopts,ip,pp,pr,td,rd,js={},z,ifs,gs},
	ip = InitialPoint/.{opts}/.Options[JuliaSet];
	pp = PlotPoints/.{opts}/.Options[JuliaSet];
	pr = Probabilities/.{opts}/.Options[JuliaSet];
	td = TreeDepth/.{opts}/.Options[JuliaSet];
	rd = RandomDepth/.{opts}/.Options[JuliaSet];
	
	
	pr = If[pr===Automatic,
			Table[1/Length[fs],{i,1,Length[fs]}],pr];
	
	(* Compile the functions *)
	gs = Map[Compile[{{u,_Complex}},Evaluate[#[u]]]&,fs];
	(* Make a function system out of them *)
	ifs = MakeFunctionSystem[fs];
	(* Make a compiled random index generator *)
	rif = RandomIndexFinder[pr];
	
	z=ip;
	While[Length[js]<pp,
		z = Iterate[gs[[rif[]]][#]&,z,rd];
		js = Flatten[js~Join~Orbit[ifs,z,td]];
		(*Print[{z,js}]*);
		];
	js
	]

(* Mandelbrot-like set for a one-parameter family of complex maps.
The function f should be a function of two complex variables, f[z,c]*)

Clear[MandelbrotSet]
Options[MandelbrotSet] = {
	EscapeRadius->2.0,
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Grays,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->100,
	PlotRange->{{-3.0,1.0},{-2.0,2.0}}
	};
mbssoptpattern = (EscapeRadius->_)|(MaxIterations->_)|
	(MaxIterationColor->_)|(ColorFunc->_)|(PlotPoints->_)|
	(ClearOut->_)|(ColorFunction->_)|(BandWidth->_);

MandelbrotSet[f_,cp_,opts___] :=	Module[
	{showopts,er,mi,mic,cf,cf2,bw,co,pp,pr,etfunc,escapetimes,dx,dy},
	er = EscapeRadius/.{opts}/.Options[FilledJuliaSet];
	mi = MaxIterations/.{opts}/.Options[FilledJuliaSet];
	mic = MaxIterationColor/.{opts}/.Options[FilledJuliaSet]/.ColorTable;
	cf2 = ColorFunction/.{opts}/.Options[FilledJuliaSet];
	bw = BandWidth/.{opts}/.Options[FilledJuliaSet];
	co = ClearOut/.{opts}/.Options[FilledJuliaSet];
	pp = PlotPoints/.{opts}/.Options[FilledJuliaSet];
	pr = PlotRange/.{opts}/.Options[FilledJuliaSet];
	dx = N[(pr[[1,2]]-pr[[1,1]])/pp];
	dy = N[(pr[[2,2]]-pr[[2,1]])/pp];
	
	showopts = DeleteCases[{opts},fjsoptpattern]~Join~
			{Frame->True,AspectRatio->dy/dx};
	cf = Function[t,If[t>=mi,List@@mic,cf2[t/bw]]];					
	If[co,Unprotect[Out];Clear[Out];Protect[Out]];
	
	(* Compile a critical orbit escape time function *)					
	etfunc = MakeCriticalEscapeTime[f,cp,er,mi];
	escapetimes = Table[
		etfunc[pr[[1,1]]+j*dx+I*(pr[[2,1]]+i*dy)],{i,0,pp},{j,0,pp}];
	
	Show[Graphics[
		Raster[Map[cf,escapetimes,{2}],Transpose[pr]]
		],showopts]
	]
	
MandelbrotSet::usage = "MandelbrotSet[f,cp,opts] plots the filled
Mandelbrot set of a family of complex polynomials using the 
critical escape time algorithm. There are many options to control 
the computation.  For example if f[z,c] := z^2+c and cp[c] = 0.0,
then: 
MandelbrotSet[f,cp,
	EscapeRadius->2.0,
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Grays,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->100,
	PlotRange->{{-2.0,2.0},{-2.0,2.0}}]
will make a plot of escape times of for the critical points of a
100 x 100 grid of parameter values,c,in the rectangle [-2,2]x[-2,2].
The points whose critical orbits do not escape in 
50 iterations will be colored black and the ones which do escape will
be colored using a function called Grays which selects a shade of
gray (other built in color function are Blues, Warms, Stripes,
BlackAndWhite). Previous outputs will be cleared to free memory. The 
BandWidth option can be increased to color groups of escape times with
the same color."

Dynamics of One-Dimensional Maps

In this section you can 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, bifurcations 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 Evaluation menu and select the item called Evaluate Initialization Cells.  It is not necessary to look in detail at these programs now.  You will see how to use them as we go along.  Then have a look at the first section below to compute some orbits and make some interesting plots.

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 used 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_] := x^2+c
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).

c=-0.5;
f^12[0.5]

You can edit previous commands and reexecute it.  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.5]

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.5];

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.366025. 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,{-2,2}]

Here is a cobweb plot with sliders for controlling a cobweb plot of the quadratic map. First a program to make a cobweb plot with given values of the parameter, c, initial state, x0, initial time t0 and orbit length, n.

QuadraticCobwebPlot[c_,x0_,t0_,n_]:= Module[{f,orb},
	f=DS[Function[x,x^2+c]];
	orb= f^{t0,t0+n}[x0];
	CobwebPlot[f,orb,{-2,2}]]

Then we can use Mathematica’s built-in Manipulate command to create sliders for an interactive plot.

Manipulate[QuadraticCobwebPlot[c,x0,t0,n],{{c,-0.5},-2,0.5},{{x0,.0},-2,2},{t0,0,100,1},{{n,50},1,100,1}]

Bifurcations and Orbit Diagrams

Saddle-node bifurcations

The saddle-node bifurcation involves the creation of two new fixed points out of nothing.  First we get one neutral fixed point with multiplier μ = 1 which then splits off into two fixed points, one attractor and one repeller.  We will investigate using the quadaratic map.

f1[x_] :=  x^2 + c
f = DS[f1]

The graph of f[x] meets the diagonal tangentially and passes through as c increases through 1/4.

c=0.3;
p1= Plot[{x,f[x]},{x,-1,1},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
c=0.25;
p2= Plot[{x,f[x]},{x,-1,1},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
c=0.2;
p3= Plot[{x,f[x]},{x,-1,1},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];

Show[GraphicsGrid[{{p1,p2,p3}}]]

Pairs of periodic points of periodic q can also be created in this way.  First we get one new crossing of the graph of f^q with the diagonal which then splits into two.  Here is a saddle node with q=3 for the quadratic map.   As the parameter increases we get six new intersections with the diagonal which are the points of the two newly created period-3 orbits

c=-1.6;
p1=Plot[{x,f^3[x]},{x,-2,2},PlotLabel->StringForm["f^3[x], f = ``",Evaluate[f["x"]]]];
c=-1.8;
p2=Plot[{x,f^3[x]},{x,-2,2},PlotLabel->StringForm["f^3[x], f = ``",Evaluate[f["x"]]]];
Show[GraphicsGrid[{{p1,p2}}]]

Here is an interactive plotter for the graphs of iterates f^q[x] with a slider to set the value of c.  You can look for bifurcations as intersections of the graphs with the diagonal.

MakeQuadraticGraph[c_,q_]:= Module[{f,orb},
	f=DS[Function[x,x^2+c]];
	Plot[{x,f^q[x]},{x,-2,2},PlotRange->{-2,2}]
	]
Manipulate[MakeQuadraticGraph[c,q],{{c,0.0},-2,0.5},{q,1,10,1}]

One can also plot the bifurcation curves for points of period q by plotting the contour g[c,x] = 0 where g is the appropriate factor of f^q[x] -x.  First the bifurcation curve for fixed points

c=.;
g= f[x]-x
p1=ContourPlot[g,{c,-3,0.5},{x,-2,2},Contours->{0},ContourShading->False,PlotPoints->30,ContourStyle->{Thick,Black}]

Here is the curve for period 3 points.  We need to eliminate the factor of f^3[x] which gives the fixed points.  You can see where the period three saddle node bifurcation occurs -- two new period three orbits (six points in all) are born.

f^3[x]-x//Factor
g=%[[2]]
p3=ContourPlot[g,{c,-3,0.5},{x,-2,2},Contours->{0},ContourShading->False,PlotPoints->30,ContourStyle->{Thick,Blue}];
Show[p1,p3]

Here is the curve for period 2 points.  There is a bifurcation where two new periodic points are created.  This is just one period two orbit so it can’t be a saddle node bifurcation.  It is a period-doubling bifurcation.

f^2[x]-x//Factor
g=%[[2]]
p2=ContourPlot[g,{c,-3,0.5},{x,-2,2},Contours->{0},ContourShading->False,PlotPoints->30,ContourStyle->{Thick,Green}];
Show[p1,p2,p3]

Period-doubling bifurcations

Here are some figures to illustrate the idea of a period-doubling bifurcation.  These occur near fixed points with multiplier μ = -1.  As a parameter varies a new period-two cycle is born.   We can see this by plotting the graphs of f and its second iterate and seeing how they cross the diagonal.  As usual we will look at the logistic map.

Here are the graphs for parameters near the period doubling at c = -3/4.  Note how the graph of the second iterate develops two extra crossings with the diagonal.

c=-0.5;
p1= Plot[{x,f[x],f^2[x]},{x,-2,2},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
c=-0.75;
p2= Plot[{x,f[x],f^2[x]},{x,-2,2},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
c=-1;
p3= Plot[{x,f[x],f^2[x]},{x,-2,2},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];

Show[GraphicsGrid[{{p1,p2,p3}}]]

The same thing can happen to get a bifurcation from period q to period 2q.  For this we can look at the graphs of f^3 and f^6 to watch the creation of period 6 orbits.  The graph is quite complicated.  We will focus in on the places where the period 3 and period 6 curves meet near the middle.  Note how the period 6 curve develops two new crossings with the diagonal.

c=-1.76;
Plot[{x,f^3[x],f^6[x]},{x,-2,2}]
p1=Plot[{x,f^3[x],f^6[x]},{x,-0.1,0.2},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
c=-1.78;
p2=Plot[{x,f^3[x],f^6[x]},{x,-0.1,0.2},PlotLabel->StringForm["f = ``",Evaluate[f["x"]]]];
Show[GraphicsGrid[{{p1,p2}}]]

Quadratic Map Critical Orbit Diagram

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.

c=.;
f1[x_] := x^2+c
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.0

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[f,ic,{-2,0},ParameterPoints->1000,OrbitPoints->200,Delay->100,PlotRange->{-2,2}]

The most interesting part of the diagram is from c=-2 to c=-1.  You can zoom in on that part by changing the range for c.

BifurcationPlot[f,ic,{-2,-1},ParameterPoints->1000,OrbitPoints->200,Delay->200,PlotRange->{-2,2},ClearOut->True]

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,{-1.635,-1.624},ParameterPoints->1000,OrbitPoints->200,Delay->100,PlotRange->{-2,2},ClearOut->True]

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,{-1.635,-1.624},ParameterPoints->2000,OrbitPoints->400,Delay->200,PlotRange->{-0.2,0.2},ClearOut->True]

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.

Periodic Orbits

The logistic map as a DS object:

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

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.85.In each case you should be able to identify the period of the attractor.  The values   a=4.0,3.6785,3.927737 have no attractors.

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 x100,...,x150 along the orbit, ignoring the first 100 iterates.

a = 3.1;
orb = f^{100,150}[0.1]
CobwebPlot[f,orb,{0,1}]

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

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]

You can try to solve f^k[x]=x by factorization for small values of k.

Factor[f^1[x]-x]
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^1[x]-x==0,x]
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 nonreal, complex numbers.  Try a=3.83, which seemed to have a period 3 attractor in the previous plots.  You will actually find two period 3 orbits -- one is a repeller. 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 graph s 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.

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.  From the graph a few steps back, it looks like both orbits have a point close to 0.5. So we will search around there.

a=3.83;
orb1= SeekPeriodicOrbit[f,3,0.5]
orb2= SeekPeriodicOrbit[f,3,0.55]

Using two different initial guesses lead 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 shows that orb1 is an attractor and orb2 is a repeller.  To illustrate some of the dynamics of this map we will compute a third orbit which starts near the repeller.  Subsequently, it gets repellled and ends up approaching the attractor.

orb3 = f^{100}[0.525]

Chaos and 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.

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

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^{100}[0.1];
orb2= f^{100}[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]]/100

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.

There is similar behavior for the tent map.

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

orb1= T^{50}[0.112];
orb2= T^{50}[0.1120001];
TimeSeriesPlot[{orb1,orb2},LineColor->{Red,Green}]
CobwebPlot[T,{orb1,orb2},{0,1},LineColor->{Red,Green}]

We can also compute a Lyapunov exponent for the tent map.

mu = Multiplier[T,orb1];
lyap = Log[ Abs[mu]]/50.0

Remarkably, this seems to be very close to the Lyapunov exponent for the logistic map with a=4.  In fact they are both equal to the natural log of 2.  It is easy to see this for the tent map and the logistic map is equivalent to the tent map by a change of coordinates.

Log[2.0]

The tent map exhibits an interesting phenomenon called “the computer can lie”.  Here are some longer orbits of the tent map.  After about 60 iterations both orbits land on the fixed point at x=0.  The cause is the binary arithmetic used in the computer.  For the tent map, every rational number initial condition of the form k/2^n eventually lands on the fixed point.  But the computer approximates all numbers in this way.

orb1= T^{100}[0.112];
orb2= T^{100}[0.1120001];
TimeSeriesPlot[{orb1,orb2},LineColor->{Red,Green}]

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, which is eventually periodic with period q=5.

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 !

Complex Dynamics and Iterated Function Systems

In this section we will look at some dynamical phenomena peculiar to complex and also plot fractals generated by

Iterated Function Systems and Fractals

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^1[triangle]
Show[Graphics[frac],AspectRatio->1,PlotRange->{{-.01,1},{-.01,1}}]

Here is an animation showing the convergence of the iterates to the attractor.  First a function to plot a certain iterate of a given IFS with a given initial shape (default is a point).

MakeFractalPlot[ifs_IFS,iterations_,shape_,opts___]:= Module[{frac},
frac=ifs^iterations[shape];
Show[Graphics[frac],AspectRatio->1,opts]
]

Now we can manipulate the number of iterations and watch the convergence.  We can also see what happens with different intial shapes.

Manipulate[ MakeFractalPlot[sierpinski,n,shape],{n,0,8,1,Appearance->"Open"},{{shape,triangle},{point->"Point",triangle->"Triangle",segment->"Segment",square->"Square",diamond->"Diamond"}},ControlPlacement->Left]

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 an animation:

Manipulate[ MakeFractalPlot[koch,n,shape,PlotRange->{{0,1},{-1/2,1/2}}],{n,1,8,1,Appearance->"Open"},{{shape,segment},{point->"Point",triangle->"Triangle",segment->"Segment",square->"Square",diamond->"Diamond"}},ControlPlacement->Left]

Complex Dynamics

Plotting Orbits in the Complex Plane

Define a complex dynamical system based on the quadratic map, just as for real maps. But here z and c are both complex numbers.

q1[z_]:= z^2+c
q=DS[q1]

We can compute iterates and orbits using the same notation as before.  But now we can use nonreal numbers for parameters and initial conditions.  Mathematics uses a capital I for the the imaginary number Sqrt[-1].

z0 = 0.5- I*0.5;
c= -0.5+I*0.2;
q^10[z0]
orb = q^{0,10}[z0]

The orbit is a sequence in the complex plane.  We can plot it using a function called ComplexPlot.  Evidently the orbit converges to an attracting fixed point near p = -0.37 + I*0.11  in the plane.

ComplexPlot[orb,PointSize->0.02, Arrows->True,PlotRange->{{-1,1},{-1,1}}]

Here is an interactive plotter with sliders for the complex map z^2 + c. First write a program to make a cobweb plot with given values of the parameter, c, initial state, z0, initial time t0 and orbit length, n.

QuadraticComplexPlot[c_,z0_,n_,opts___]:= Module[{f,orb},
	f=DS[Function[x,x^2+c]];
	orb= f^{0,n}[z0];
	ComplexPlot[orb,PlotRange->{{-2,2},{-2,2}},PointSize->0.01,opts]]

Then we can use Mathematica’s built-in Manipulate command to create sliders for an interactive plot.

Manipulate[QuadraticComplexPlot[a+I*b,z0[[1]]+I*z0[[2]],Ceiling[n],ImageSize->500,Arrows->arrows,MoreGraphics->{Graphics[Text[{a,b},{-1.8,1.8}]]}],{{z0,{0.0,0.0}},Locator},{{a,-0.5},-2.3,0.7,Appearance->"Open"},{{b,0.0},-1.5,1.5,Appearance->"Open"},{n,1,30,1,Appearance->"Open"},{arrows,{True,False}},ControlPlacement->Left]

Filled Julia Sets

Define a complex dynamical system based on the quadratic map.  z and c are both complex numbers.

q1[z_]:= z^2+c
q=DS[q1]

The filled Julia set is the complement of the basin of attraction of the fixed point at z = Infinity.  Any point which exceeds the escape radius of Max[ 2, Abs[c] ] is in the basin.  Other points can be colored according to how many interations are required to reach the escape radius.  The points which do not escape constitute the filled Julia set.  The function FilledJuliaSet implements this scheme.

Here are some intereresting ones.  For c values in the “main cardiod” in the c-plane, the map has an attracting fixed point.

c=-0.5;
FilledJuliaSet[q,
	EscapeRadius->Max[2.0,Abs[c]],
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Blues,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->400,
	PlotRange->{{-2.0,2.0},{-2.0,2.0}}]//Timing

There should be an attracting fixed point in the interior of this Julia set.  The theory says it must attract the critical orbit with initial conditions z=0. So we can find it by computing the tail of the critical orbit.  We will add that to the plot. You can try some of the other suggested values for c.

c=-0.5  ;  (* -0.5, -0.8,-1.0, -1.3, -0.75, 0.3-0.4*I, -0.1+0.8*I, 0.360284+I*0.100376 *)
fjs= FilledJuliaSet[q,
	EscapeRadius->Max[2.0,Abs[c]],
	MaxIterations->100,
	MaxIterationColor->Black,
	ColorFunction->Blues,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->400,
	PlotRange->{{-2.0,2.0},{-2.0,2.0}}];
q^{100,120}[0.0+I*0.0];
critorb = MakePoints[%,{1,0,0},0.01];
Show[fjs,critorb]

Julia Sets by Inverse Iteration

The Julia set is completely invariant, meaning that the preimages of any point in it are also in it.  Starting from any point in J we can compute many others by following the branches of the inverse function system.  Here we use the two inverse branches of z^2+c.

Here are the  branches of the complex square root and the roots of z-c

Sign2[x_]:=If[x≥0,1.0,-1.0]
ComplexSquareRoot1[z_]:= With[{x=Re[z],y=Im[z],r=Sqrt[Re[z]^2+Im[z]^2]},Sqrt[(r+x)/2]+I*N[Sign2[y]*Sqrt[(r-x)/2]]]
ComplexSquareRoot2[z_]:= With[{x=Re[z],y=Im[z],r=Sqrt[Re[z]^2+Im[z]^2]},-Sqrt[(r+x)/2]-I*N[Sign2[y]*Sqrt[(r-x)/2]]]

ComplexSquareRootC1[z_]:= With[{x=Re[z-c],y=Im[z-c],r=Sqrt[Re[z-c]^2+Im[z-c]^2]},Sqrt[(r+x)/2]+I*N[Sign2[y]*Sqrt[(r-x)/2]]]
ComplexSquareRootC2[z_]:= With[{x=Re[z-c],y=Im[z-c],r=Sqrt[Re[z-c]^2+Im[z-c]^2]},-Sqrt[(r+x)/2]-I*N[Sign2[y]*Sqrt[(r-x)/2]]]

Make an iterated function system out of the two branches.

ComplexSqrt=IFS[ComplexSquareRoot1,ComplexSquareRoot2]
ComplexSqrtC=IFS[ComplexSquareRootC1,ComplexSquareRootC2]

Now we can easily find a list of roots of  roots of roots ...

ComplexSqrt[1]
ComplexSqrt^5[1]

For an initial point in J we can use an unstable fixed point.

FindUnstableFixedPoint[c_]:= With[{fps =({1,1}+ ComplexSqrt[1-4*c])/2},
	If[Abs[2*fps[[1]]]>1,fps[[1]],fps[[2]]]]

Here is a function which finds 2^n of point in the Julia set for z^2+c

JuliaPoints[c0_,n_]:= With[{z0=FindUnstableFixedPoint[c0]},
c=c0;
ComplexSqrtC^n[z0]]

ComplexToXY[Complex[x_,y_]]:={x,y}

jp = JuliaPoints[-1.0+0.0*I,15];
jp2=Map[ComplexToXY,jp];
ListPlot[jp2,PlotStyle->PointSize->0.005,PlotRange->{{-2,2},{-2,2}},AspectRatio->1]

Mandelbrot Set

c=.;
f[z_,c_]:= z^2+c
cp[c_]:=0.0+0.0*I;

ms=MandelbrotSet[f,cp,
	EscapeRadius->2.0,
	MaxIterations->50,
	MaxIterationColor->Black,
	ColorFunction->Grays,
	BandWidth->1,
	ClearOut->True,
	PlotPoints->300,
	PlotRange->{{-2.25,0.75},{-1.5,1.5}}]