Feigenbaum constant
Computing Feigenbaum constant
So I saw this on Stephen Wolfram’s blog where a simple recursive equation can yield chaotic behaviour and seems to have properties like Hopf bifurcation. I wanted to check out myself.
Reference: https://blog.stephenwolfram.com/2019/07/mitchell-feigenbaum-1944-2019-4-66920160910299067185320382/
from sympy import * ## I would love to do it like this. Unfortunately I cannot find a nice way to make a recursive ## sequence in sympy. So I went for the pythonic approach. If someone has a better idea let me ## know! # a, x = symbols('a x', real=True) # i = symbols('i', integer=True) # f, g = symbols('f g', cls=Function) init_printing()
from functools import lru_cache @lru_cache(maxsize=8) def f(i, a, x0): if i == 0: return x0 else: return a * f(i-1, a, x0) * (1 - f(i-1, a, x0)) # sequence(f, (x, 0, 10))
def series(n, a, x0=1/3): return [f(i, a, x0) for i in range(n)] series(30, 3.2)
from statistics import mean def series_mean(n, a, x0=1/3): s = series(n, a, x0) return s, mean(s[n//3:])
N = 50 fig, axes = plt.subplots(2, 2, sharex=True, sharey=True) for a, ax in zip( (2, 3.2, 3.4, 3.5), axes.ravel() ): s, smean = series_mean(N, a) ax.plot(s) ax.set_title(f'$a = {a}$') ax.hlines(smean, N//3, N, linestyles='dashed')
In the last subplot we begin to see period doublings.
Bifurcation diagram
In the blog the following Wolfram code is used to generate the calculate the bifurcation. The first 50 values of the series are ignored to avoid transients and upto 300 values are calculated vfor each value of a.
ListPlot[Flatten[ Table[{a, #} & /@ Drop[NestList[Compile[x, a x (1 - x)], N[1/3], 300], 50], {a, 0, 4, .01}], 1], Frame -> True, FrameLabel -> {"a", "x"}]
Before we do that, let us see if it makes any difference if we vary the initial condition.
import numpy as np N = 300 a_values = np.linspace(0, 4, 100) x0_values = np.linspace(0.1, 5, 100) plt.figure() for x0 in x0_values: smeans = [series_mean(N, a, x0)[1] for a in a_values] plt.scatter(a_values, smeans, s=1, c="r") plt.xlabel("a") plt.ylabel("series average")
/usr/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: overflow encountered in double_scalars
The average value where the series oscillates around does not seem to depend on the value of x0. So now instead of plotting the mean, we can plot the full distribution where of values where the series oscillates around.
import numpy as np N = 300 a_values = np.linspace(0, 4, 100) plt.figure() for a in a_values: s = series(N, a)[50:] plt.scatter(a * np.ones_like(s), s, s=1, c="r") plt.xlabel("a") plt.ylabel("series distribution")
There is our multiple Hopf bifurcation :) Let us see what happens for a > 3.5.
N = 300 fig, axes = plt.subplots(2, 2, sharex=True, sharey=True, figsize=(15, 5)) for a, ax in zip( (3.6, 3.7, 3.8, 3.9), axes.ravel() ): s, smean = series_mean(N, a) ax.plot(s) ax.set_title(f'$a = {a}$') ax.hlines(smean, N//3, N, linestyles='dashed')
Nice! Maybe I will do a follow up to compute the Lyapunov constant.
You candownloadthis notebook, or see a static viewon nbviewer.
About the author
Ashwin Vishnu Mohanan, Ph.D. in Fluid mechanics