iterated maps on interval [0,1]

xinwei
2 aprile 1989

(* hard coded to plot f[x] only in unit square [[0,1]x[0,1]]  *)

BeginPackage["Iter`"]

Feigenbaum::usage = "Feigenbaum[f,{x0:0,x1:1,dx:.1},orbitLength:5] computes orbits of a few points in the interval [x0,x1] < [0,1] under iterates of f[x]"

forbit::usage = "forbit[f,x,n] generates a list of the first n iterates of f on x."

Staircase::usage = "Staircase[A_List] makes a staircase of points interpolated into a sequence of (x,y) points."

DrawStaircase::usage = "DrawStaircase[A_List,blue_Real] draws a staircase of points interpolated into a sequence of (x,y) points, in a color specified by blue"

Begin["`Private`"]

Attributes[Staircase] = {HoldFirst}

Staircase[A_List] = Flatten[
					Table[ {{A[[i]],A[[i]]},{A[[i]],A[[i+1]]}}, 
					{i,1,Length[A] - 1}],
					1
					]

Attributes[DrawStaircase] = {HoldFirst}

DrawStaircase[A_List,blue_Real] := ListPlot[Staircase[A], 
					PlotStyle->{RGBColor[0.3,.5,Sin[blue]^2],PointSize[.01]},
					PlotJoined->True]

forbit[f_,x_,n_] := Block[{fn,u,nn},
				fn[0,u_] := u;
				fn[1,u_] := f[#1]&[u];
				fn[nn_,u_] := fn[nn,u] = f[#1]&[fn[nn - 1,u]];
				Table[fn[nn,x],{nn,0,n}]
				]

Attributes[Feigenbaum] = {HoldFirst}

Feigenbaum[f_,{x0_:0.,x1_:1.,dx_:0.1},orbitLength_Integer:5] :=
	Feigenbaum[f,{0.,1.,0.1}]

Feigenbaum[f_,{x0_:0.,x1_:1.,dx_:0.1},orbitLength_Integer:5] :=

	Block[ {z,fn,forbit,g,flow,glow},

(* define iterates *)

	fn[0,u_] := u;
	fn[1,u_] := f[#1]&[u];
	fn[n_,u_] := fn[n,u] = f[#1]&[fn[n - 1,u]];

	forbit[u_,n_] := Table[fn[nn,u],{nn,0,n}];

	dc =  dx (Pi/2)/(x1 - x0);

(* debug *)
		
	Print["x0 = ", x0];
	Print["x1 = ", x1];
	Print["dx = ", dx];
	Print["dc = ", dc];
	Print["orbitLength = ", orbitLength];
	Print["f[z] == ", f[#1]&[z]];

	g = Plot[{x,f[#1]&[x]},{x,N[0],N[1]},
		PlotRange->{{0,1},{0,1}},			(* HARD CODED *)
		PlotStyle->{RGBColor[.7,.7,.7],RGBColor[1,0,.2]},
		AxesLabel->{"x","f[x]"}
		];

(* compute orbits *)
	glow = {};
	
	Do[Append[glow,
		Block[{k = 1,gk},
		flow[k] := forbit[xx,orbitLength];
		Print["k = ", k, ";	xx = ",xx,";	orbit = ", N[flow[k]]];
		gk = DrawStaircase[flow[k], N[Pi/4 + dc] ];
		Show[{g,gk}];
		k++;
		Return[gk] ]
		],
	{xx,N[x0],N[x1],N[dx]}
	];
(*
	Print["Show glow(1) "];

	Show[glow[1]];

	Print["DoÉShow[glow[kk]]"];
	
	Do[Show[{g,Table[glow[kk]]}]
	,{kk,1,Length[glow]}
	];
*)
	Return[glow]  (* returns array of graphs *)
]

End[]
EndPackage[]

Null