(* Functions for Interating functions and finding orbits*) Iterate[f_,x_,n_] := Nest[f,x,n] Iterate::usage = "Iterate[f,x0,n] gives the result of applying the n-th iterate of the function f to the initial value x0." Orbit[f_,x_,n_,nskip_:0] := NestList[f,Nest[f,x,nskip],n] Orbit::usage = "Orbit[f,x0,n,nskip] returns the part of the orbit of x0 from iterate nskip to iterate n+nskip. If no nskip is given, it defaults to 0 and the result is just the orbit segment{x0,x1,...xn}. The return value is a list." (* Functions for turning orbits into Mathematica graphics objects *) (* 2D graphics*) MakePoints[orb_,color_:{0,0,0},size_:0.01] := Graphics[{Apply[RGBColor,color],PointSize[size],Map[Point,orb]}] MakeLine[orb_,color_:{0,0,0},thickness_:0.005] := Graphics[{Apply[RGBColor,color],Thickness[thickness],Line[orb]}] MakePolygon[orb_,color_:{0,0,0},size_:0.01] := Graphics[{Apply[RGBColor,color],PointSize[size],Polygon[orb]}] MakeDiagonal[x_,color_:{0,0,0},thickness_:0.005] := Graphics[{Apply[RGBColor,color],Thickness[thickness], Line[{{x[[1]],x[[1]]},{x[[2]],x[[2]]}}]}] 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]}]]}]] (* 3D versions *) MakePoints3D[orb_,color_:{0,0,0},size_:0.01] := Graphics3D[{Apply[RGBColor,color],PointSize[size],Map[Point,orb]}] MakeLine3D[orb_,color_:{0,0,0},thickness_:0.005] := Graphics3D[{Apply[RGBColor,color],Thickness[thickness],Line[orb]}] MakePolygon3D[orb_,color_:{0,0,0},size_:0.01] := Graphics3D[{Apply[RGBColor,color],PointSize[size],Polygon[orb]}] (* 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}}; (* Functions which transform 1D orbits into 2D lists for plotting *) CobwebOrbit[orb_] := With[{dup=Flatten[Transpose[{orb,orb}]]}, Transpose[{Drop[dup,-1],Drop[dup,1]}]] GraphOrbit[orb_] := Map[Flatten,Transpose[{Range[0,Length[orb]-1],orb}]] EmbedOrbit[orb_,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}]] ] ] (* A generic plotting function for 1D systems to be used later for defining specialized plots *) Options[GenericPlot] = { Points->True, Lines->True, PointColor->Black, LineColor->Black, PointSize->0.01, LineSize->0.002, Transformation->GraphOrbit, MoreGraphics->{} }; genoptpattern = (Points->_)|(Lines->_)|(PointColor->_)| (LineColor->_)|(PointSize->_)|(LineSize->_)| (Transformation->_)|(MoreGraphics->_); (* A version for plotting a single orbit *) GenericPlot[orb_List/;VectorQ[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] ] (* And a version for plotting a list of orbits *) GenericPlot[orbs:{{___}..},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] ] (* Here are 3D versions *) GenericPlot3D[orb_List/;VectorQ[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,MakeLine3D[trorb,lc,ls],{}], If[pf,MakePoints3D[trorb,pc,ps],{}] }~Join~mg; Show[gr,showopts] ] GenericPlot3D[orbs:{{___}..},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,MakeLine3D[#1,#3,#4],{}]&, {trorbs,lf,lc,ls}], MapThread[If[#2,MakePoints3D[#1,#3,#4],{}]&, {trorbs,pf,pc,ps}] }]~Join~mg; Show[gr,showopts] ] (* Cobweb Plots *) CobwebPlot[f_,orb_,xrange_,opts___] := GenericPlot[orb,opts,Transformation->CobwebOrbit,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 already defined and orb1 and orb2 have already been computed using Orbit[...], then CobwebPlot[f,{orb1,orb2},{0.0,1.0},LineColor->{Black,Red}] will make a plot showing the two orbits in the given colors." (* Time Series Plots *) TimeSeriesPlot[orb_,opts___] := GenericPlot[orb,opts,Transformation->GraphOrbit, 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." (* Delay Embedding Plots *) DelayPlot[orb_,dim_Integer:2,delay_Integer:1,opts___] := Switch[dim, 2,GenericPlot[orb,opts,Lines->False, Transformation->(EmbedOrbit[#,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 TimeSeriesPlot[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." (* Graph of iterate(s) of a map *) IteratePlot[f_,n_Integer,xrange_] := Plot[{x,Iterate[f,x,n]},{x,xrange[[1]],xrange[[2]]}, AspectRatio->1,PlotRange->{xrange,xrange}, PlotLabel->StringForm["f = `` n = ``", Evaluate[f["x"]],n]] IteratePlot[f_,n_List,xrange_] := Plot[Evaluate[{x}~Join~Map[Iterate[f,x,#]&,n]], {x,xrange[[1]],xrange[[2]]}, AspectRatio->1,PlotRange->{xrange,xrange}, PlotLabel->StringForm["f = `` n = ``", Evaluate[f["x"]],n]] IteratePlot::usage = "IteratePlot[f,n,xrange] produces a graph of the n-th iterate of f over the given range of x values. Giving a list of n values produces several super- imposed graphs. For example, if f is already defined IteratePlot[f,{1,2,7},{0,1}] plots the first, second and seventh iterates." (* Find the minimal period of a sequence to a given tolerance. It returns Infinity if no period is found. *) MinimalPeriod[orb_,tol_:10^(-6)] := With[{diff = Map[Max[Abs[#]]&, Drop[orb-Table[orb[[1]],{Length[orb]}],1]]}, Min[Position[diff,_?(Function[d,d<=tol])]] ] MinimalPeriod::usage = "MinimalPeriod[orb,tol] tries to find the minimal period of the sequence orb by searching for the an item in the sequence which agrees with the first item to the given tolerance. The entries of orb may be vectors." (* Random point in a given range *) (* 1D version returns a number *) RandomPoint[r_List/;VectorQ[r]] :=Random[Real,r]; (* Vector version returns a list. The range is a pair of lists of vectors r ={{xmin,ymin,...},{xmax,ymax,...}} *) 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,...}}." Options[SeekPeriodicOrbit] = { Radius->Automatic, MaxTries->10, Tolerance->0.000001, MaxIterations->100 }; spooptpattern = (Radius->_)|(MaxTries->_)|(Tolerance->_); SeekPeriodicOrbit[f_,q_,xguess_,opts___] := Module[{fq,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}} ]; fq = Iterate[f,x,q]-x//Simplify; sol = FindRoot[fq,##]& @@ Flatten[{ic,fropts},1]; x0 = x/.sol; po = Orbit[f,x0,q]; If[MinimalPeriod[po]===q,Return[po]]; (* Failing this, try random guesses nearby *) If[ra===Automatic, ra=2*Max[Abs[xguess-Iterate[f,xguess,q]]]]; 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[fq,##]& @@ Flatten[{ic,fropts},1]); po = Orbit[f,x0,q]; If[MinimalPeriod[po]===q,Return[po]]; mt--]; (* Failing this, return the empty list *) Return[{}]; ] SeekPeriodicOrbit::usage = "SeekPeriodicOrbit[f,q,xguess,opts] tries to find a periodic orbit of period q for the map f near the point xguess using Mathematica's FindRoot function. Works well only for small q and low dimensional maps. If it succeeds, it returns the list {x0,x1,...,xq} of length q+1. If it fails, it returns the empty list {}." (* Computing multipliers of periodic points *) Der[f_,x0_] := Module[{x},D[f[x],x]/.x->x0] Jac[f_,x0_] := Module[{x,xx}, x = Array[xx,Length[x0]]; Transpose[ Map[ D[f[x],#]&, x]/. Table[xx[i]->x0[[i]],{i,1,Length[x0]}] ] ] OrbitJacobian[f_,orb_List/;VectorQ[orb]] := If[Length[orb]===1,1, Apply[Times,Map[Der[f,#]&,Drop[Reverse[orb],1]]] ] OrbitJacobian[f_,orb:{{___}..}] := If[Length[orb]===1,IdentityMatrix[Length[orb[[1]]]], Apply[Dot,Map[Jac[f,#]&,Drop[Reverse[orb],1]]] ] Multiplier[f_,orb_List/;VectorQ[orb]] := OrbitJacobian[f,orb] Multiplier[f_,orb:{{___}..}] := With[{A=OrbitJacobian[f,orb]},Eigenvalues[A]] Multiplier::usage = "Multiplier[f,orb] computes the multiplier or multipliers of a periodic orbit, orb = {x0,...,xq} 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." "Iteration. \n A Mathematica package for iterating maps. \n (C) 1997 Richard Moeckel.\n Licensed to Dr. Hector Lomeli"