3
$\begingroup$

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$?

$\endgroup$
17
  • 3
    $\begingroup$ Ha! another stupid decision from SO 10K users stackoverflow.com/q/16261682/353410 $\endgroup$ – Dr. belisarius Apr 29 '13 at 14:38
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M.'s torpor Apr 29 '13 at 14:43
  • 2
    $\begingroup$ @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. $\endgroup$ – Mark McClure Apr 29 '13 at 15:01
  • 1
    $\begingroup$ @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. $\endgroup$ – Mark McClure Apr 29 '13 at 15:10
  • 1
    $\begingroup$ @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 $\endgroup$ – Mark McClure Apr 29 '13 at 21:26
6
$\begingroup$

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

enter image description here

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.