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 plotted in the bifurcation diagram, 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]

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]tolerance

Now try a smaller tolerance value:

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

Detecting whether or not a trajectory appears to have a k-cycle (Algorithm 3.21 in the text)

Use the vector comparison method above to write a module that determines whether a given trajectory appears to have 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. (We will implement Algorithm 3.21 from the text.)

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

Test your code here:

hasKCycle[traj,2,0.01]

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

hasKCycle[traj,4,0.1]

hasKCycle[traj,4,0.01]

hasKCycle[traj,4,0.000001]

Searching for a bifurcation value (Algorithm 3.22 in the text)

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, which is Algorithm 3.22 in the text:

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],0.3,numIter];
If[hasKCycle[traj,k,tolerance],left=mid,right=mid]
];
mid
]

Now the following code will search for the largest value of r where a 2-cycle appears:

findBifurcation[f,10000,2,0.00001,3,4]//N

The following code will search for the largest value of r where a 4-cycle appears:

findBifurcation[f,100000,4,0.000001,3.44,4]//N

Similarly, find more bifurcation values for the logistic map.

Examine the period-doubling cascade that starts with 3-cycles

Recall (from the homework) the values of r that produce 3-cycles. Use the code above to locate the bifurcations in the period-doubling cascade that produces 6-cycles, 12-cycles, 24-cycles, etc. Then use these bifurcation values to estimate the Feigenbaum constant.

Explore bifurcations arising from another function!

The function (Subscript[f, r])(x)=r x(1-x) is not the only function that produces behavior from which we can observe period-doubling cascades and chaos. Replace (Subscript[h, r])(x)=x (e^(r(1-x))) for r>0 and x>=0. Create a bifurcation diagram and find bifurcation values. Can you approximate the Feigenbaum constant from this function?

h[r_][x_]:=x Exp[r(1-x)]
Manipulate[
Plot[h[r][x],{x,0,10},PlotRange->All],
{r,0,10}]

pts=bifurcationDiagramPoints[h,0.3,0,5,0.001]//Quiet;
Graphics[ {PointSize[0.001],Opacity[0.2], Point[pts]}, AspectRatio->1/2,Frame->True]