# Nest for large values of n

I am working on a problem to find the Feigenbaum constant. In this problem, I have to get $f^{2^n}(1/2)$ with $f(x) = kx(1-x)$. I use Nest to get $f^{2^n}(1/2)$, but it takes Mathematica forever to run with $n>5$; my goal is to get to $n=14$.

Does anybody know any function or trick to get $f^{2^n}(1/2)$ with a large value of $n$?

• Ha! another stupid decision from SO 10K users stackoverflow.com/q/16261682/353410 – Dr. belisarius Apr 29 '13 at 14:38
• You do know that iterating the logistic equation eventually results in a periodic sequence, no? Find the period, and then you can extrapolate on what the $2^n$-th iterate ought to be. – J. M.'s torpor Apr 29 '13 at 14:43
• @belisarius Try k=2.9 or 3.1 or 3.4. You happen to be choosing a value of k right at a bifurcation, where the convergence to the fixed point is absurdly slow. But, of course, the iteration need not always converge to a periodic orbit; that's the point behind chaos. – Mark McClure Apr 29 '13 at 15:01
• @belisarius No, I don't think you can analytically find the attractive orbit in terms of k in a simple way. What I think J.M. is suggesting is extrapolate from when the bifurcations occur. In fact, the existence of the Feigenbaum constant is really an assertion that such an extrapolation works. – Mark McClure Apr 29 '13 at 15:10
• @thewanderer I assume that you mean when n>5, but note that for n=5, your polynomial has degree 2^2^5. and it's simply not feasible to solve this using simple techniques. You should look into some of the papers of Keith Briggs to learn his extrapolation techniques: keithbriggs.info – Mark McClure Apr 29 '13 at 21:26

This is from an old notebook of mine. Not claiming it's the best nor even optimal. Surprisingly, it's well commented. I wrote it in a step-by-step fashion for demonstrational purposes for some class.

$$\text{\delta =}\lim_{n\rightarrow \infty } \delta _n=\lim_{n\rightarrow \infty } \frac{r_{n+1}-r_n}{r_{n+2}-r_{n+1}}=\lim_{n\rightarrow \infty } \frac{\mu _{n+1}-\mu _n}{\mu _{n+2}-\mu _{n+1}}$$

where $r=4\mu$.

## Initial conditions

Calculating $\mu_0$, $\mu_1$ and $\mu_2$; defining F and delta[n] for higher n.

f[x_] := mu x (1 - x)


(* Numerically, find the values μ_n at which the 2^n cycle is superstable, for the first few values of n. *)(* Superstable orbits contain 1/2. *)

(* We choose an initial value for x_0 to find bifurcation points; we chose x_0=1/2, because why not *)

eq0 = f[1/2] == 1/2;
mu0 = mu /. N[Solve[eq0, mu]][[1]]


2.

m = {};
AppendTo[m, mu0]


{2.}

eq1 = f[f[1/2]] == 1/2;
mu /. N[Solve[eq1, mu]]


{2., -1.23607, 3.23607}

(* There are three roots here. We can use FindRoot to find a numerical solution; if we start near the expected solution, perhaps we'll get it. *)

max = First[Select[mu /. N[Solve[eq1, mu]], # > m[[1]] &]]


3.23607

FindRoot[f[f[1/2]] == 1/2, {mu, max}];
mu1 = mu /. %


3.23607

AppendTo[m, mu1]


{2., 3.23607}

(* To get the later bifurcation points μ[n], we need a function which iterates f many times. In Mathematica, the standard method (Nest) actually expands the composition out algebraically before plugging in: iterating f 2^12 times involves generating an expression of order 2^12-1 long in mu. We use old-fashioned iteration to speed things up. Unfortunately, Mathematica does things algebraically unless we insist otherwise. We insist that x and mu be real numbers before iterating! *)

eq2 = f[f[f[f[1/2]]]] == 1/2;
mu2 =
First[
Select[
Sort[
Cases[
mu /. N[Solve[eq2, mu]], _Real
]
], # > m[[2]] &
]
]


3.49856

AppendTo[m, mu2]


{2., 3.23607, 3.49856}

F[x_Real, mu_Real, NTimes_] := Block[{y = x}, Do[y = mu y (1 - y), {NTimes}]; y]


Example:

mu /. FindRoot[F[0.5, mu, 2^2] == 0.5, {mu, Last[m]}]


3.49856

or within an interval:

mu /. FindRoot[F[1/2, mu, 2^2] == 1/2, {mu, 3.45, mu1, 3.55}]


3.49856

(* As n increases, there are more and more roots to (f^(12))[1/2]=1/2. To find the right one, we have to be careful. What we can do is use the fact that the roots are spaced in an approximate geometric progression: μ[n]~μ[n-1]+(μ[n-1]-μ[n-2])/δ *)

delta[n_] := (m[[n]] - m[[n - 1]])/(m[[n + 1]] - m[[n]])

delta[2]


4.70894

## Further recurrences

c[n_] := m[[n]] + (m[[n]] - m[[n - 1]])/delta[n - 1]

mu[n_] := mu[n] = mu /. FindRoot[F[1/2, mu, 2^n] == 1/2, {mu, c[n]}]

Do[AppendTo[m, mu[n]], {n, 3, 11}]

m
m // Length


{2., 3.23607, 3.49856, 3.55464, 3.56667, 3.56924, 3.5698, 3.56991, 3.56994, 3.56994, 3.56995, 3.56995}

12

No need to go to bigger n as there's no additional convergence.

## Onset of chaos is at r=...

Last[m]


3.56995

## Feigenbaum constant

Table[delta[n], {n, 2, Length[m] - 1}]


{4.70894, 4.68077, 4.66296, 4.6684, 4.66895, 4.66916, 4.66919, 4.6692, 4.6692, 4.6692, 4.6692, 4.66919, 4.66919}

Last[%]


4.66919

## Plots

plot1 = ListPlot[Table[{n, m[[n]]}, {n, 1, Length[m]}], Frame -> True,
PlotStyle -> {Black, PointSize[Medium]}, PlotRange -> All,
AxesOrigin -> {0, 2},
FrameLabel -> {"n",
Rotate["\!$$\*SubscriptBox[\(r$$, $$n$$]\)", 270 Degree]}]
plot2 = ListPlot[Table[{n, delta[n]}, {n, 2, Length[m] - 1}],
Frame -> True, PlotStyle -> {Black, PointSize[Medium]},
PlotRange -> All,
FrameLabel -> {"n",
Rotate["\!$$\*SubscriptBox[\(\[Delta]$$, $$n$$]\)", 270 Degree]}]


(* δ is converging nicely to the universal value. *)