Program 2 - Cobweb

Program Code

Off[General::"spell1", General::"spell"];

BeginPackage["Cobweb`"];

Unprotect[Cobweb];

Cobweb::usage = "Cobweb[f, x0, numtokeep, numtodrop] creates a cobweb diagram for the function f 
	starting at x0 using numtokeep iterations after dropping numtodrop values.";
	
Options[Cobweb] = {
	LinePlotStyle -> {Hue[2/3]},
	FPlotStyle -> {},
	PlotRange->All
};

Begin["`Private`"];

Cobweb::badfunc = "`2` detected in function `1`.";
Cobweb::badinit = "Initial value (`1`) is not a real number."; 
Cobweb::badits = "Number of iterations to keep (`1`) is not a positive integer."; 
Cobweb::badign = "Number of iterations to ignore (`1`) is not a non-negative integer."; 
Cobweb::badopt = "Options expected (instead of `1`) beyond position 4 in argument list."; 
Cobweb::badval = "An error occurred after `1` iterations.  Using values up to this point.";
Cobweb::baderr = "An error occurred while computing the first `1` iterations.  Unable to continue.";

RuleQ[r_] := Head[r] === Rule;
RealQ[r_] := NumericQ[r] && Im[r] == 0;
RealMachineQ[r_] := MachineNumberQ[r] && Im[r] == 0 || Abs[r] < 1.0;

Cobweb[f_, x0_, its_ , ign_, args___] := Module[{g, var, x, xmin, xmax, lineopts, fopts, pr, temp, points, i, plotopts},	
	Check[f[x0], Return[]];
	Switch[Head[f],
		Function, g = f,
		Symbol, If[MemberQ[Attributes[f], NumericFunction], g = f, g = #&],
		_, var = DeleteCases[Union[Flatten[Apply[List, {f}, {0, Infinity}]]], _?RealQ];
			Switch[Length[var],
				0, Message[Cobweb::badfunc, f, "No variables"]; Return[],
				1, If[Head[var[[1]]] === Symbol, 
					g[x_] := f /. var[[1]]->x,
					Message[Cobweb::badfunc, f, "Unknown symbol"]; Return[]
				],
				_Integer, Message[Cobweb::badfunc, f, "Too many variables"]; Return[]
			];
	];
	If[!RealQ[x0], Message[Cobweb::badinit,x0]; Return[]];
	If[Head[its] =!= Integer || its < 1, Message[Cobweb::badits, its]; Return[]];
	If[Head[ign] =!= Integer || ign < 0, Message[Cobweb::badign, ign]; Return[]];
	If[i = 0; Scan[(i++; If[Head[#1] =!= Rule, Return[True]])&, opts],
		Message[Cobweb::badopt, opts[[i]]]; Return;
	];
	
	{lineopts, fopts, pr} = {LinePlotStyle, FPlotStyle, PlotRange} /. {args} /. Options[Cobweb];
	plotopts = DeleteCases[{args}, (LinePlotStyle->_) | (FPlotStyle->_) | (PlotRange->_)];
	
	Off[General::"unfl"];
	x = Nest[g, N[x0], ign];
	If[!RealMachineQ[x], Message[Cobweb::baderr, ign]; Return[]];
	points = Flatten[Table[{{x, x = g[x]}, {x, x}}, {its}], 1];
	
	points = Select[points, VectorQ[#, RealMachineQ]&];
	If[!RealMachineQ[g[points[[-1,1]]]], points = Delete[points, -1]];
	
	If[Head[pr] === List && Length[pr[[1]]] == 2,
		{xmin, xmax} = pr[[1]],
		temp = Transpose[points];
		{xmin, xmax} = {Min[temp[[1]]], Max[temp[[1]]]}
	];
	If[xmax - xmin < 10^-16, xmin -= 1; xmax += 1];
	On[General::"unfl"];
	
	Return[Show[
		Plot[g[x], {x, xmin, xmax}, PlotStyle->fopts, DisplayFunction->Identity],
		Plot[x, {x,xmin,xmax}, DisplayFunction->Identity],
		Graphics[{Sequence@@lineopts, Line[points]}],
		PlotRange->pr, plotopts, DisplayFunction->$DisplayFunction
	]];
];

End[];

SetAttributes[Cobweb, ReadProtected];
Protect[Cobweb];

EndPackage[];

On[General::"spell1", General::"spell"];

Examples

Cobweb[3.83  x (1 - x), .3, 100, 0];

It appears that the orbit is displaying some kind of erratic behavior. However, it may be that we are simply looking at the transient part of the orbit.

Cobweb[3.83  x (1 - x), .3, 100, 60];

As we see, the first plot contained the transient part of the orbit.  In the long run the orbit converges to a periodic orbit of period 3.

Cobweb[x^2, 0.99, 10, 0];

Convergence to the fixed point of 0 for f(x) = x^2

Cobweb[Cos, 0, 100, 0];

Fixed point of Cosine.