Tools for Analyzing Discrete Dynamical Systems:
Bifurcation Diagrams
Written by Richard A. Holmgren
last modified June 18, 1996

These programs were developed for use with the text  A First Course in Discrete Dynamical Systems, 2nd edition by Richard A. Holmgren, which was published by Springer-Verlag in 1996.  They are much faster than the programs described in the appendix of the book.    Suggestions for improvement are gratefully accepted.  Please contact Professor Holmgren at rholmgre@alleg.edu with suggestions and comments.

Bifurcation Diagrams  --- Example 1.

Description

The following two sets of commands  will create a bifurcation diagram for a parametrized family of functions.  In this example, we are plotting a bifurcation diagram for Newton's function for the family of cubics q(x) = x3 + p x + 1.  (This is Figure 12-10 on page 148 of the text.)  Notice that the function defined in the program has two variables, x and the parmater value p.
	In the first set of code, the list of ordered pairs  data is created.  Each ordered pair represents a point on the bifurcation diagram.  To reduce computational time and to minimize the size of the list  data, the function bifurcList computes and temporarily stores the first retained iterate value as startPoint for each parameter value  If under continued iteration, we return to this starting value, creation of the list for that parameter value stops, since additional points would be redundant.  To ensure that startPoint is close to the attracting set (if such a set exists), one should choose  firstIt to be relatively large.  Note that there are tradeoffs in selecting firstIt. On the one hand, as  firstIt gets large, the amount of computation necessary to find startPoint increases.  On the other hand, as startPoint becomes larger, the list  data will become shorter.  It is worth experimenting with a small number of parameter values (controlled by psteps in the program) and iterates before using large values for psteps and iterates.   In the example shown, computation time can be decrease by about 30% by setting firstIt to 250.  On the other hand the length of data increased by 25%.  To decrease storage size I chose to use the higher value for firstIt.  In this example I creae a high resolution image of the bifurcation diagram.  A reasonably good set of data for a diagram can be created in about 20 seconds by setting firstIt equal to 200, iterates at 100, psteps at 150.  
	The second set of code  creates the bifurcation diagram by using ListPlot to plot the list data.

Creating the list data, the points in the bifurcation diagram.

(*  p is the parameter   *)

Clear[f,x,firstIt,iterates,seed,pmin,pmax,p,
	psteps,step,bifurcList,data];

f = Compile[{{x,_Real},{p,_Real}}, (2 x^3 - 1)/(3x^2 + p)];

(* firstIt is the first iteration shown, 
	iterates is the maximum number of iterations 
		calculated after firstIt,
	seed is the starting point of the iteration. *)

firstIt = 500;
iterates = 150;
seed = 0.0;

(*  pmin and pmax define the bounds of the parameters,
		psteps is the number of paramter values used.  *)

pmin =  -1.3;
pmax =  -1.25;
psteps = 400;

p = N[pmin];
step = N[(pmax - pmin)/psteps];

bifurcList[p_] := Module[
	{startPoint = Nest[f[#,p]&,seed,firstIt]},
	Map[{p,#}&,
			FixedPointList[f[#,p]&,
				startPoint,
				iterates,
				SameTest -> (#2 == startPoint &)]]]

Timing[
data = 	Partition[
			Flatten[
				Table[
					bifurcList[t],
					{t,pmin,pmax,N[(pmax-pmin)/psteps]}]
			],
		2];
]

Length[data]

The following code plots the data stored in the list data:

(*  verticalMin and verticalMax define 
	the bounds of the plotted values.  *)

verticalMin =  -.29;
verticalMax =  .29;

Timing[
	ListPlot[data,
		PlotRange -> {{pmin,pmax},{verticalMin,verticalMax}},
		AxesOrigin->{pmin,verticalMin},
		PlotStyle->{PointSize[.002]},
		Frame -> True]
	]

Bifurcation Diagrams --- Example 2.
	(avoiding overflow errors)

Description

This example uses the same code as the previous example to plot a bifurcation diagram for f(x,p) = p x (1-x).  The function Chop[...] is used in defining f to avoid overflow errors.

Finding the points in the bifurcation diagram.

(*  p is the parameter   *)

Clear[f,x,firstIt,iterates,seed,pmin,pmax,
	p,psteps,step,bifurcList,data];

f = Compile[{x,p},Chop[p x (1-x)]];

(* firstIt is the first iteration shown, 
	number of iterates calculated after firstIt,
	seed is the starting point of the iteration. *)

firstIt = 250;
iterates = 150;
seed = .5;

(*  pmin and pmax define the bounds of the parameters,
		psteps is the number of paramter values used.  *)

pmin =  2;
pmax =  4;
psteps = 400;

p = N[pmin];
step = N[(pmax - pmin)/psteps];

bifurcList[p_] := Module[
	{startPoint = Nest[f[#,p]&,seed,firstIt]},
	Map[{p,#}&,
			FixedPointList[f[#,p]&,
				startPoint,
				iterates,
				SameTest -> (#2 == startPoint &)]]]

Timing[
data = 	Partition[
			Flatten[
				Table[
					bifurcList[t],
					{t,pmin,pmax,N[(pmax-pmin)/psteps]}]
			],
		2];
]

Length[data]

The following code plots the data stored in the list data:

(*  verticalMin and verticalMax define 
		the bounds of the plotted values.  *)

verticalMin =  0;
verticalMax =  1;

Timing[
	ListPlot[data,
		PlotRange -> {{pmin,pmax},{verticalMin,verticalMax}},
		AxesOrigin->{pmin,verticalMin},
		PlotStyle->{PointSize[.002]},
		Frame -> True]
	]