Experiment: Feigenbaum's Constant

Goal: In the previous lab we used the orbit diagram to see that the quadratic function (Subscript[Q, c])(x) = (x^2) + c, the logistic function (Subscript[F, c])(x) = c x(1 - x), and the sine function c sin(x) all undergo a sequence of period-doubling bifurcations as the parameter tends to the chaotic regime. We have also seen that magnifications of the orbit diagram tend to look "the same." In this experiment, we will see that there really is some truth to this: we will see that these period doubling bifurcations always occur at the same rate.
 
Procedure: In this experiment you will work with either the quadratic or the logistic family. We first need a definition:

Definition. Suppose (Subscript[x, 0])is a critical point for F, that is F' ((Subscript[x, 0])) = 0. If (Subscript[x, 0])is also a periodic point of F with period n then the orbit of (Subscript[x, 0])is called superstable. The reason for this terminology is that ((F^n))' ((Subscript[x, 0])) = 0. 

	In this experiment we will first determine the c-values at which either (Subscript[Q, c])(x) = (x^2) + cor (Subscript[F, c])(x) = c x(1 - x)have superstable cycles of periods 1, 2, 4, 8, 16, 32, and 64. Using the programs, experimentally determine the c-values at which your function has a superstable point of the given period. Be sure to check that this point has the correct prime period. You should be looking in the Period-1 window of the corresponding diagram for these points.
	There are a number of ways to do this. For example, you could list 1000 points on the orbit of a random initial condition. By changing the parameter, you should then search for the "exact" cfor which the critical point (0 for (Subscript[Q, c])(x) = (x^2) + c, 0.5 for (Subscript[F, c])(x) = c x(1 - x)) lies on a cycle. You will usually not be able to find this parameter value exactly. However, you should find the c value for which you come closest to having the critical point periodic. Usually, this means that you find a point on the orbit within (10^(-6))of the critical point. You should find c accurate to seven decimal places, which is more or less the accuracy of the computations.
	Another approach might be to compute the first (2^n)points on the orbit of the critical point, and then seeing how close you come to this point. Modify the parameter repeatedly to try to come closer to the value for which the critical point is periodic with the right period. After finding the seven c-values for your function, record these numbers in tabular form:
	
1.	(Subscript[c, 0])= c-value for period (2^0)
2.	(Subscript[c, 1])= c-value for period (2^1)
3.	(Subscript[c, 2])= c-value for period (2^2)
4.	(Subscript[c, 3])= c-value for period (2^3)
5.	(Subscript[c, 4])= c-value for period (2^4)
6.	(Subscript[c, 5])= c-value for period (2^5)
7.	(Subscript[c, 6])= c-value for period (2^6)

	Now use a calculator to compute the following ratios:
	
			(Subscript[f, 0]) = (((Subscript[c, 0]) - (Subscript[c, 1]))/((Subscript[c, 1]) - (Subscript[c, 2]))),  
          (Subscript[f, 1]) = (((Subscript[c, 1]) - (Subscript[c, 2]))/((Subscript[c, 2]) - (Subscript[c, 3]))),  ... ,  
          (Subscript[f, 4]) = (((Subscript[c, 4]) - (Subscript[c, 5]))/((Subscript[c, 5]) - (Subscript[c, 6]))).
List these numbers in tabular form, too. Do you notice any convergence? You should, at least if you have carried out the above search to enough decimal places.

Notes and Questions: 
1. The limit as n goes to infinity of these ratios is called the Feigenbaum's constant. It turns out that this number is "universal" -it appears whenever a typical family undergoes the period-doubling route to chaos. Compare the results of both the quadratic and logistic families. Are they the same? The fact that these numbers are always the same is a remarkable result due to Feigenbaum several years ago.

2. This lab takes a long time to perform. It is useful to work with friends and divide up the tasks.

Solving the equations

We can obtain some values of c by solving the equations using Reduce in the following way:

Equation = Function[x, (x^2) + c];
Reduce[Equation[0] == 0, c]

Then, (Subscript[c, 0]) = 0. Now we can continue with the second iterate (n = 2) to obtain (Subscript[c, 1]).

n = 2;
NIterate[x_] := Nest[Equation, x, n]
NIterate[x]

N[Reduce[NIterate[0] == 0, c], 20]

Since c = 0was a solution for f(x) (Equation), we know that (Subscript[c, 1]) = -1. To find (Subscript[c, 2])we need to solve the equation for the fourth iterate of f(x) which is a (2^3) = 8  th degree polynomial in c (this can be done by changing the value of n in the definition of NIterate. There are 8 roots: (Subscript[c, 0]) = 0, (Subscript[c, 1]) = -1, four complex roots and two real roots. The bigger of these two real roots is (Subscript[c, 2]). It is the bigger one (or the smallest in absolute value) because this is the one that is closer to the cascade of bifurcations we are following.

If we do this one more time for the 8th iterate of f(x) then we have to wait until the program finds the (2^7) = 128 roots! (this can be a lot of time depending on the machine we are using). After that, we can make a similar analysis to the one done before (discard complex roots and roots already found and keep the largest real root which in this case is (Subscript[c, 3])=-1.3815474844320614695. Though we found four of the c-values with this method we have to use a different approach for the other 3, since the time involved in the calculations for the 16th iterate of f(x) can be gigantic. (Mathematica has to approximate the (2^15)roots!). 

We can try to find the place where the c-axis cuts the bifurcation diagram at the given period by looking at portions of the orbit diagram. Using the code below, you will see that the line of x = 0 has been colored red to make it easier to visualize. The next value of c is on the intersection of this line with the bifurcation diagram at the period 16 bifurcation. We can start looking at a closer portion of the diagram in the interval -1.5 < c < -1 (defined by Initialc and Finalc), but we know that the next intersection should happen at a c smaller thand (Subscript[c, 3]). So we just zoom into the region -1.5<c<-1.38. 

After this first zoom we see that we can make a better approximation if we plot the diagram around -1.4 < c < -1.38. To get a better visualization and approximation to the value of c we can change the limits of the plot on the x-axis (vertical axis). To do this you have to change the values of xInf and xSup.  If we zoom in the region given by - .25 < x < 0.25,  -1.4 < 
        c < -1.38we'll see two intersections with the red line, i.e., the one corresponding to (Subscript[c, 3])and the one corresponding to (Subscript[c, 4])that lies in the interval -1.4 < c < -1.395. So now we know (Subscript[c, 4])to two decimal places and we can continue with this method to get an approximation with the required accuracy. The next c-values can be handled in the same way. 

When the diagram becomes unclear you should increase the number of Iterates and Divisions as you zoom in to have a more neat picture.

Using the orbit diagram

Note: when trying to plot the bifurcation diagram with Finalc = 0 we get an underflow error message because iterates of the Random[] seed (which lies between zero and one) under Q(x)  = (x^2)goes to values smaller than the limits assumed by Mathematica ($MinNumber) . If we want to avoid the message we can use Finalc = -(10^(-3)) instead of 0, for example.

Equation = Function[x, (x^2) + c];
FunctionName = "Quadratic equation";
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 14};
Initialc = -2; Finalc = - .001;
xInf = -2; xSup = 2;
Iterates = 500;
DropIterates = 400;
Divisions = 500;

OrbitFunction[c_] := 
        Take[NestList[Equation, Random[], 
            Iterates - 
              1], -(Iterates - 
                DropIterates)];
DataSet = 
        Table[{c, OrbitFunction[c]}, {c, Initialc, 
            Finalc, (Finalc - Initialc)/(Divisions + 
                  1)}];
DataPlot = 
        Flatten[Table[{DataSet[[n]][[1]], DataSet[[n]][[2]][[i]]}, {n, 1, Length[DataSet]}, {i, 1, 
              Iterates - DropIterates}], 1];
OrbitDiagram = 
        ListPlot[DataPlot, Frame ->  True, Axes -> False, 
          PlotRange -> {{Initialc, Finalc}, {xInf, xSup}}, 
          PlotJoined -> False, 
          PlotStyle -> 
            AbsolutePointSize[
              0], 
PlotLabel -> "Orbit Diagram for \
the " <> FunctionName <>  ",  " <> 
                ToString[Initialc] " ≤ c ≤ " <> 
                ToString[Finalc], FrameLabel -> {"c", "x"}, 
          RotateLabel -> False, ImageSize -> 500, 
          DisplayFunction -> Identity];
Show[{OrbitDiagram, 
          Graphics[{RGBColor[1, 0, 0], Line[{{-2, 0}, {0, 0}}]}, 
            DisplayFunction -> Identity]}, 
        DisplayFunction -> $DisplayFunction];