Detect intervals in which period-doubling biffurcations occur

Question

Let $F_\mu (x) = \mu x (1-x)$ be the logistic map. Our goal is to find -numerically- intervals $[\mu_i, \mu_{i+1}]$ in which a period-doubling biffurcation occurs, for some unique $\mu \in [\mu_i, \mu_{i+1}]$.

My first approach to the problem has been to generate a table with the points of the periodic sink of $F_\mu$, as a function of $\mu$ (i.e., after a suitable transient, I have saved about 500 images by $F_\mu$ of some initial condition point), using the following code:

def logistic_iter(mu, x, trans, iter, four):

    for i in range(trans):
        x = mu * x * (1.-x)

    for i in range(iter):
        x = mu * x * (1.-x)
        fout.write("%+.24e, %+.24e\n" % (mu, x))

for $\mu \in (0,4]$.

However, I am not guessing an approach to compute the intervals in which period-doubling bifurcations do occur; I don't know if I could calculate them using the table I have generated, or maybe I should use the conditions on $F_\mu$ and its derivatives stated by the period-doubling bifurcation theorem. Thanks in advance for your help.

Answer 1

Close to bifurcation points and in general for large $N$, the minimal distance between elements of a stable $N$ cycle may become rather small. To avoid to detect a cycle closure too early, one can also use the stability condition.

A fixed point equation $x_{+1}=g(x)$ is stable as some fixed point $x^*$ if $|g'(x^*)|<1$.

The derivative $(F_\mu^N)'(x)$ can be computed alongside the fixed-point iteration. If $x_0$ is the point obtained after the initial iteration, one can compute

x=x0
dx = 1
for k in range(...):
    dx =mu*(1-2*x)*dx
    x = mu*x*(1-x)

For a candidate for the cycle closure we can now test that both $x-x_0$ is small and the approximate cycle is stable.

    if abs(x-x0)<1e-6 and abs(dx)<1: N=k+1; break
x0 = x 

Now in most cases this will only give an approximate cycle. One has to employ Newton's method or a Newton-like method to improve the accuracy. Close to a bifurcation this can lead to the cycle actually turning out to be unstable. Correcting for that leads into fiddly details...

In the region where there is a stable $N$-cycle, $|(F_\mu^N)'(x_\mu)|<1$, this derivative value falls almost linearly from $1$ to $-1$ for increasing $\mu$, the crossing of $-1$ signifying the location of the next bifurcation to a cycle of length $2N$. One can easily find the points where the derivative is $0$ and $-1$ via the secant method. The curve at $1$ has a parabolic maximum, at least for even $N$, as $(F_\mu^{2N})'(x_\mu)=[(F_\mu^N)'(x_\mu)]^2$ implies a parabolic maximum.

Starting from $\mu=3.2$ in the region of a 2-cycle, one can thus progress from bifurcation point to the next. With $N=2^{13}$ this starts to take some time. This gives values \begin{array}{r|lll} N&\mu:(F_\mu^N)'(x_\mu)\approx 0&\mu:(F_\mu^N)'(x_\mu)=-1& a:\Delta\mu=a\Delta x \\ \hline 2 & 3.24047619890513 & 3.44948974278318 & 0.20412414139881577 \\ 4 & 3.49861586518712 & 3.54409035955192 & 0.04392655500239954 \\ 8 & 3.55465275837676 & 3.56440726609543 & 0.009410300867575267 \\ 16 & 3.56666992883554 & 3.56875941954383 & 0.002015741040611005 \\ 32 & 3.56924407788826 & 3.56969160980140 & 0.0004317211496063787 \\ 64 & 3.56979541068891 & 3.56989125937812 & 9.246144305075316e-05 \\ 128 & 3.56991349041541 & 3.56993401837398 & 1.980214101486517e-05 \\ 256 & 3.56993877953466 & 3.56994317604841 & 4.240729755916831e-06 \\ 512 & 3.56994419569071 & 3.56994513734217 & 9.079502953480227e-07 \\ 1024 & 3.56994535566589 & 3.56994555739125 & 1.9417135274185843e-07 \\ 2048 & 3.56994560409669 & 3.56994564735291 & 4.130240115155148e-08 \\ 4096 & 3.56994565730132 & 3.56994566661991 & 8.565943654997352e-09 \\ \end{array} The middle values are the bifurcation points for the dyadic periods, the last column is the inverse slope, or approximately half the interval length of the period $N$ segment.

For the non-dyadic cycles one needs either a starting point or some better idea of how to locate them. Using the gap locations for instance in WP diagram, this procedure results in a shorter tables as the narrow windows increase numerical instability \begin{array}{r|lll} N&\mu:(F_\mu^N)'(x_\mu)\approx 0&\mu:(F_\mu^N)'(x_\mu)=-1& a:\Delta\mu=a\Delta x \\ \hline 6 & 3.62690666691446 & 3.63037957812559 & 0.0034729112111314495 \\ 12 & 3.63129386571999 & 3.63219721958149 & 0.0009033538615054917 \\ 24 & 3.63240148551568 & 3.63261417127399 & 0.00021268575831604147 \\ 48 & 3.63265081026205 & 3.63269676101533 & 4.595075328671703e-05 \\ 96 & 3.63270498983568 & 3.63271489743776 & 9.907602071416682e-06 \\ 192 & 3.63271661879498 & 3.63271873854568 & 2.119750696781524e-06 \\ \\ 7 & 3.70171104892548 & 3.70214110016727 & 0.00043005124178910234 \\ 14 & 3.70228342657636 & 3.70241190326916 & 0.0001284766928026278 \\ 28 & 3.70244285267251 & 3.70247372689324 & 3.087422073228223e-05 \\ 56 & 3.70247901911882 & 3.70248569961662 & 6.680497797708733e-06 \\ 112 & 3.70248690310852 & 3.70248834356056 & 1.4404520400778342e-06 \\ \\ 5 & 3.73867876928011 & 3.74099037061598 & 0.0023116013358749003 \\ 10 & 3.74184900649612 & 3.74257705325005 & 0.0007280467539310246 \\ 20 & 3.74274607265424 & 3.74291925537727 & 0.00017318272302537012 \\ 40 & 3.74294922011552 & 3.74298672717352 & 3.750705800854383e-05 \\ 80 & 3.74299344528601 & 3.74300153468941 & 8.089403397565585e-06 \\ 160 & 3.74300294083368 & 3.74300467177685 & 1.7309431729890515e-06 \\ 320 & 3.74300497614715 & 3.74300534690940 & 3.7076225170432725e-07 \\ \\ 3 & 3.82976051137067 & 3.84145908621077 & 0.011698574840095875 \\ 6 & 3.84457154070478 & 3.84763967570011 & 0.00306813499533137 \\ 12 & 3.84834717612974 & 3.84907380702609 & 0.0007266308963513067 \\ 24 & 3.84919810905193 & 3.84935495108779 & 0.0001568420358588566 \\ 48 & 3.84938313482655 & 3.84941697715290 & 3.384232634927233e-05 \\ 96 & 3.84942284877557 & 3.84943008801832 & 7.239242751146846e-06 \\ 192 & 3.84943136109451 & 3.84943291293459 & 1.5518400786578022e-06 \\ \end{array}