Feigenbaum’s universal constant in the period-doubling cascade

As the parameter \(r\) grows, the logistic map \(x \mapsto rx(1-x)\) doubles its period again and again – 1, 2, 4, 8, … – at parameter values whose spacing shrinks by a constant factor, Feigenbaum’s \(\delta \approx 4.6692\). Universality means the same \(\delta\) appears for any smooth one-humped map. This script measures \(\delta\) for the logistic map and for the unrelated sine map \(x \mapsto r\sin(\pi x)\) and shows both sequences of ratios converging to the one constant.

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import brentq, newton

from mathematicskit.ode_dynamics.systems.logistic_map import bifurcation_diagram, estimate_feigenbaum_delta
from mathematicskit.ode_dynamics.visualizers.plots import plot_bifurcation_diagram

FEIGENBAUM_DELTA = 4.669201609


def logistic(x, r):
    return r * x * (1.0 - x)


def sine_map(x, r):
    return r * np.sin(np.pi * x)

Library estimate from the first bifurcations

estimate_feigenbaum_delta() locates the first three period doublings of the logistic map by scanning \(r\) – a rough first look at the constant.

print(f"estimate_feigenbaum_delta() -> {estimate_feigenbaum_delta():.4f}  (true delta = {FEIGENBAUM_DELTA:.4f})")

estimate_feigenbaum_delta() -> 4.7453  (true delta = 4.6692)

Precise ratios from superstable cycles

Between consecutive doublings sits a superstable parameter \(R_n\) where the map’s maximum \(x_c = 1/2\) lies on the period-\(2^n\) cycle, i.e. \(f_R^{2^n}(x_c) = x_c\). These roots are easy to find accurately, and their spacings shrink by the same ratio: \(\delta_n = (R_{n} - R_{n-1}) / (R_{n+1} - R_n) \to \delta\).

def superstable_parameters(f, R0, R1_bracket, n_max=8):
    def g(R, n):
        x = 0.5
        for _ in range(2**n):
            x = f(x, R)
        return x - 0.5

    R = [R0, brentq(g, *R1_bracket, args=(1,), xtol=1e-14)]
    for n in range(2, n_max + 1):
        delta_guess = 4.67 if n == 2 else (R[-2] - R[-3]) / (R[-1] - R[-2])
        guess = R[-1] + (R[-1] - R[-2]) / delta_guess
        R.append(newton(g, guess, args=(n,), tol=1e-14, maxiter=200))
    R = np.array(R)
    return R, (R[1:-1] - R[:-2]) / (R[2:] - R[1:-1])


R_log, delta_log = superstable_parameters(logistic, 2.0, (3.0, 3.44))
R_sin, delta_sin = superstable_parameters(sine_map, 0.5, (0.72, 0.83))
print(" n   delta_n (logistic)   delta_n (sine map)")
for n, (dl, ds) in enumerate(zip(delta_log, delta_sin), start=1):
    print(f"{n:2d}   {dl:17.5f}   {ds:17.5f}")

n   delta_n (logistic)   delta_n (sine map)
1             4.70894             4.04574
2             4.68077             4.55585
3             4.66296             4.64518
4             4.66840             4.66407
5             4.66895             4.66811
6             4.66916             4.66897
7             4.66919             4.66915

Two different maps, one cascade, one constant

r_values = np.linspace(2.8, 4.0, 1500)
r_plot, x_plot = bifurcation_diagram(r_values, n_transient=500, n_keep=100)

s_values = np.linspace(0.65, 1.0, 1500)
x = np.full_like(s_values, 0.5)
for _ in range(500):
    x = sine_map(x, s_values)
s_orbit = np.empty((100, s_values.size))
for k in range(100):
    x = sine_map(x, s_values)
    s_orbit[k] = x

fig, axes = plt.subplots(1, 3, figsize=(14, 4.5))
plot_bifurcation_diagram(r_plot, x_plot, ax=axes[0])
axes[0].set_title(r"Logistic map $rx(1-x)$")
axes[1].plot(np.tile(s_values, 100), s_orbit.ravel(), ",", color="black", alpha=0.5)
axes[1].set_xlabel("r")
axes[1].set_ylabel("x (after transient)")
axes[1].set_title(r"Sine map $r\sin(\pi x)$")
for ax, R in ((axes[0], R_log), (axes[1], R_sin)):
    for i, Rn in enumerate(R[1:5]):
        ax.axvline(Rn, color="tab:red", lw=0.6, alpha=0.7, label=r"superstable $R_n$" if i == 0 else None)
    ax.legend(loc="lower left", fontsize=8)

n = np.arange(1, delta_log.size + 1)
axes[2].plot(n, delta_log, "o-", label="logistic map")
axes[2].plot(n, delta_sin, "s--", label="sine map")
axes[2].axhline(FEIGENBAUM_DELTA, color="k", lw=1, ls=":", label=r"Feigenbaum $\delta = 4.6692$")
axes[2].set_xlabel("n")
axes[2].set_ylabel(r"$\delta_n = (R_n - R_{n-1}) / (R_{n+1} - R_n)$")
axes[2].set_title(r"Universality: both ratios $\to \delta$")
axes[2].legend(fontsize=8)
fig.tight_layout()

plt.show()