The Feigenbaum Constant

Math 242 Modern Computational Mathematics

Bifurcation values of the logistic map

First we define the logistic map as usual:

f[r_][x_]:=r*x(1-x)

Here is a version of our bifurcation diagram code from last time. Unlike last time, this module returns the list of points, rather than the graphics object. Below, we will add some more graphics elements before displaying the graphics.

bifurcationDiagramPoints[f_,start_,rStart_,rStop_,rStep_,numIter_:1000,numTail_:100]:=Module[{traj,tail,points,allPoints={}},
Do[
(* compute a trajectory of (Subscript[f, r]) *)
traj=NestList[f[r],start,numIter];

(* take only the tail of the trajectory *)
tail=traj[[-numTail;;]];

(* convert the tail to a list of points, and join them to the list of points *)
points=Table[{r,x},{x,tail}];
allPoints=Join[allPoints,points],

(* repeat for r from rStart to rStop in steps of rStep *)
{r,rStart,rStop,rStep}];

(* return the list of all points *)
allPoints
]

Next we will display a portion of the logistic map bifurcation diagram, with some of the bifurcation value marked.

bifDiagram=bifurcationDiagramPoints[f,0.3,2.95,3.6,0.0002]; (* the bifurcation diagram *)

(* here are some bifurcation values, computed to several decimal places *)
bifVals={3.0,3.4495,3.54409,3.56441,3.568759,3.5696916,3.5698913,3.5699340};
lines=Table[Line[{{r,0},{r,1}}],{r,bifVals}];   (* draw lines at some bifurcation values *)
labels=Table[Text[r,{r,0},{-1,-1},{0,1}],{r,bifVals[[;;4]]}];  (* label the lines *)

(* now display the graphics *)
Graphics[ {PointSize[0.001],Opacity[0.1], Point[bifDiagram],Red,Opacity[1],lines,labels}, AspectRatio->1/2,Frame->True]

Today we will explore the ratios of the distances between successive bifurcation values.

First, compute the distances between successive bifurcation values:

diffs=Table[bifVals[[i]]-bifVals[[i-1]],{i,2,Length[bifVals]}]

Now compute the ratios between successive distances:

ratios=Table[diffs[[i-1]]/diffs[[i]],{i,2,Length[diffs]}]

Notice that these ratios are almost the same! In fact, they are converging to a particular value: 4.669201..., which is known as the Feigenbaum constant.

The computations above are based on the first period-doubling cascade that occurs for the logistic map as the parameter r increases. However, the logistic map exhibits many other period-doubling cascades. What will we find if we compute the ratios of distances between bifurcation values for other period-doubling cascades?

Computing bifurcation values precisely

In order to compute bifurcation values to several decimal places, we must be able to accurately determine the cycle type for any given trajectory. Specifically, given a trajectory and a positive integer k, we want to determine whether the trajectory converges to a k-cycle. We’ve done this before: our basic strategy is to compute the trajectory to several hundred terms and then see if the last k terms (approximately) repeat. Today we’ll do this in a slightly more systematic way than before.

Comparing vectors up to some tolerance value

Here is a trajectory:

traj=NestList[f[3.3],0.3,20]

Take the last k terms of the trajectory, along with the k terms before those:

k=2;
vec1=traj[[-k;;]]
vec2=traj[[-2k;;-k-1]]

We want to determine whether vec1 and vec2 are approximately equal. In other words, we want to determine whether the difference of the vectors is close to zero — specifically, whether the absolute value of each component of the difference is less than some tolerance value.

Here is the absolute value of the vector difference:

Abs[vec1-vec2]

Abs[vec1-vec2]<0.02

We want to determine whether these values are each less than some tolerance, such as 0.01 or 0.0001. In Mathematica, the usual less-than and greater-than operators don’t compare vectors. Instead, we need to use Mathematica’s VectorLess comparison operator. This operator can be typed by the key sequence Escv<Esc, which Mathematica converts to \[VectorLess]. Here are some examples:

{1,2,3}\[VectorLess]{2,3,4}
{1,2,3}\[VectorLess]{2,2,2}

Now let's use the \[VectorLess] operator to compare |vec1-vec2| with a vector of tolerance values.

tolerance=0.1;
Abs[vec1-vec2]\[VectorLess]Table[tolerance,Length[vec1]]

Norm[vec1-vec2]<tolerance

Now try a smaller tolerance value:

tolerance=0.001;
Abs[vec1-vec2]\[VectorLess]Table[tolerance,Length[vec1]]

Detecting whether or not a trajectory has a k-cycle

Use the vector comparison method above to write a module that determines whether a given trajectory has a k-cycle. The module requires three parameters: a trajectory, a positive integer k, and a tolerance value. The module should return True if the last k values of the trajectory are within the tolerance value of the previous k values; otherwise the module returns False.

hasKCycle[traj_,k_,tolerance_]:=Module[{vec1,vec2},
vec1=traj[[-k;;]];
vec2=traj[[-2k;;-k-1]];
Abs[vec1-vec2]\[VectorLess]Table[tolerance,Length[vec1]]
]

Test your code below:

traj=NestList[f[3.5],0.3,30]

hasKCycle[traj,4,0.1]

hasKCycle[traj,4,0.0000001]

hasKCycle[traj,2,0.1]

hasKCycle[traj,8,0.1]

Searching for a bifurcation value

Suppose we know that the bifurcation from k-cycles to 2k-cycles occurs in some interval (Subscript[r, left])<r<(Subscript[r, right]). We want to find the largest value of r for which we detect a k-cycle. We will use a bisection method to zero in on this bifurcation value. Specifically, we test whether we observe a k-cycle at the midpoint of our interval: if so, we know that the bifurcation value is in the right half of the interval; if not, we know the bifurcation value is in the left half of the interval. By repeatedly subdividing the interval and testing the midpoints, we can estimate the bifurcation value quite well.

Here is pseudocode for the algorithm:

Algorithm: findBifurcation
Input: function f; parameters numIter, k, tolerance, (Subscript[r, left]), and (Subscript[r, right])
Set left=(Subscript[r, left]) and right=(Subscript[r, right]).
Repeat while right-left>tolerance:
	Set mid=(left+right)/2.
	Compute trajectory at r=mid for numIter iterations.
	If the trajectory has a k-cycle, then set left=mid; otherwise set right=mid.
Return: mid

Implement the algorithm below:

findBifurcation[f_,numIter_,k_,tolerance_,rLeft_,rRight_]:=Module[{left=rLeft,mid,right=rRight,start=0.3,traj},
While[right-left>tolerance,
mid=(right+left)/2;
traj=NestList[f[mid],start,numIter];
(*Print[mid,traj]*);
If[hasKCycle[traj,k,tolerance],left=mid,right=mid]
];
(*return mid*)
mid
]

Try it out. Detect where the 2-cycle bifurcates to a 4-cycle:

findBifurcation[f,1000,2,0.001,3.4,3.6]

Detect where the 16-cycle bifurcates to a 32-cycle:

findBifurcation[f,1000,16,0.00001,3.5,3.6]

If we compute enough terms of the trajectory and choose a tolerance that is small enough (but not too small), we can determine the cycle type.

Examine the period-doubling cascade starting with the 3-cycle

We observed that the logistic map trajectories exhibit 3-cycles for r near 3.83. Compute the bifurcation values of the period-doubling cascade that starts with a 3-cycle. Then compute the ratios between the differences of these bifurcation values. What do you observe?

Determine where 3-cycle bifurcates into 6-cycle:

s[1]=findBifurcation[f,10000,3,0.00001,3.83,3.9]

Determine where 6-cycle bifurcates into 12-cycle:

s[2]=findBifurcation[f,10000,6,0.00001,3.84,3.9]

Determine where 12-cycle bifurcates into 24-cycle:

s[3]=findBifurcation[f,30000,12,0.000001,3.83,3.9]

Determine where 24-cycle bifurcates into 48-cycle:

s[4]=findBifurcation[f,30000,24,0.000001,3.83,3.9]

Determine where 48-cycle bifurcates into 96-cycle:

s[5]=findBifurcation[f,30000,48,0.0000001,3.83,3.9]

Estimate the ratios:

diffs=Table[s[i]-s[i-1],{i,2,5}]

ratios=Table[diffs[[i-1]]/diffs[[i]],{i,2,4}]

Explore other functions

Many other families of functions can be iterated to produce trajectories and bifurcation diagrams. For each of the following functions, create a bifurcation diagram. Then identify period-doubling cascades and compute bifurcation values. Finally, compute the ratios of successive distances between successive bifurcation values. What do you observe?

(Subscript[g, r])(x)=r sin(π x), for 0<=r<=1 and 0<=x<=1

g[r_][x_]:=r*Sin[Pi*x]

Manipulate[Plot[g[r][x],{x,0,1},PlotRange->{0,1}],{r,0,1}]

bifDiagram=bifurcationDiagramPoints[g,0.3,0,1,0.0002]; (* the bifurcation diagram *)

(* now display the graphics *)
Graphics[ {PointSize[0.001],Opacity[0.1], Point[bifDiagram]}, AspectRatio->1/2,Frame->True]

(Subscript[h, r])(x)=x (e^(r(1-x))) for r>0 and x>=0

See Table 3.1 in the text for even more functions!