Here we will investigate some of the dynamics for one-dimensional maps.

If you have a printed version, please take a look at the coloured version on my website at www.will-farmer.com.

In [238]:
import numpy as np
import sympy as sp
from sympy import symbols
import matplotlib.pyplot as plt
import seaborn
import networkx as nx
import pydot
import ipywidgets
import pandas as pd

from IPython.display import display

sp.init_printing()
In [239]:
%matplotlib inline
In [240]:
%load_ext cython
The cython extension is already loaded. To reload it, use:
  %reload_ext cython

Beginning the Exploration - Cobweb Diagrams

We will use the graph of the logistic map

$$ x_{n+1} = \mu x_n (1 - x_n) $$

and the diagonal $y=x$. This map expresses chaotic behavior for certain values of $\mu$. We can examine "orbits" of this system by looking at what values the map bounces around to. A cobweb diagram is a good way to see these "orbits".

Task 1. Find the parameter values at which the stable period $2^1$, $2^2$, and $2^3$ orbits are first created. Label these $\mu_1$, $\mu_2$, $\mu_3$.

In [241]:
%%cython
import numpy as np
cimport numpy as np

def cobweb(f, int n=100, int start=0, float initial = 0.5):
    """ Generate the path for a cobweb diagram """
    cdef np.ndarray[np.float64_t, ndim=2] web = np.zeros((n, 2))
    web[0, 0] = initial
    web[0, 1] = initial
    cdef int state = 1
    for i in range(1, n):
        if state:
            web[i, 0] = web[i - 1, 0]
            web[i, 1] = f(web[i - 1, 0])
        else:
            web[i, 0] = web[i - 1, 1]
            web[i, 1] = web[i - 1, 1]
        state ^= True
    return web[start:]

Using the Mac tool, we can find these values to be as follows. Note however, that this is by no means an exact process, and reliability should be questioned. We can recreate the mac tool by using ipywidgets in the IPython Notebook.

In [243]:
x = np.linspace(0, 1, 100)

@ipywidgets.interact(mu=(1, 4, 0.01))
def plot(mu=3):
    f = lambda x: mu * x * (1 - x)
    web = cobweb(f, n=1000)
    plt.figure(figsize=(6, 6))
    plt.plot(x, x)
    plt.plot(x, f(x))
    plt.plot(web[:, 0], web[:, 1], linewidth=0.5)
    plt.show()

And here they are.

In [6]:
x = np.linspace(0, 1, 100)
web1 = cobweb(lambda x: 3.069946 * x * (1 - x), n=1000, start=100)
web2 = cobweb(lambda x: 3.449945 * x * (1 - x), n=1000, start=500)
web3 = cobweb(lambda x: 3.549946 * x * (1 - x), n=1000, start=500)

plt.figure(figsize=(6, 6))
plt.plot(x, x)
plt.plot(x, 3.069946 * x * (1 - x), 'b-')
plt.plot(x, 3.449945 * x * (1 - x), 'g-')
plt.plot(x, 3.549946 * x * (1 - x), 'r-')
plt.plot(web1[:, 0], web1[:, 1], 'b-', linewidth=0.5, label=r'$\mu_1$')
plt.plot(web2[:, 0], web2[:, 1], 'g-', linewidth=0.5, label=r'$\mu_2$')
plt.plot(web3[:, 0], web3[:, 1], 'r-', linewidth=0.5, label=r'$\mu_3$')
plt.legend(loc=0)
plt.show()

Bifurcation Diagrams

We've already noted that the accuracy of finding these bifurcation points was low, let's instead examine a bifurcation diagram. A bifurcation diagram is essentially a probabilistic view of our map for different values of $\mu$. For the following plots, the $x$-axis is differing values of $\mu$, and the $y$-axis is a large number of plotted values after the transient.

In [104]:
%%cython

import numpy as np
cimport numpy as np
cimport cython

@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function
def bifurcation(np.int64_t precision=1000, np.int64_t keep=500, np.int64_t num_compute=10000,
                np.float64_t xmin=0, np.float64_t xmax=4, np.float64_t ymin=0, np.float64_t ymax=1):
    """ Acquire bifurcation points for varying mu for logistic map """
    cdef np.ndarray[np.float64_t, ndim=1] mu = np.linspace(xmin, xmax, precision, dtype=np.float64)
    cdef np.float64_t x = 0.5
    cdef np.int64_t i, j, k
    cdef np.ndarray[np.float64_t, ndim=2] points = np.zeros((len(mu) * keep, 2), dtype=np.float64)
    k = 0
    for i in range(len(mu)):
        for j in range(num_compute):
            x = mu[i] * x * (1 - x)
            if j > (num_compute - keep): # we throw away the transient
                points[k, 0] = mu[i]
                points[k, 1] = x
                k += 1
    return points

The full bifurcation diagram for the logistic map follows.

In [8]:
points = bifurcation(xmin=1, xmax=4)
plt.figure(figsize=(12, 4))
plt.plot(points[:, 0], points[:, 1], ',', color='k', alpha=0.8)
plt.xlim(1, 4)
plt.show()

Now that we can plot the bifurcation diagram, let's examine the first several bifurcations.

In [9]:
mu_vals = np.array([3, 3.45, 3.544, 3.5645, 3.56875, 3.5697, 3.5698925, 3.569934, 3.56994316, 3.5699452, 3.569945646])
In [10]:
def plot_bifurcation(fig, axarr, index, x, y, xmin, xmax, ymin, ymax, precision, keep, num):
    points = bifurcation(precision=precision, xmin=xmin, xmax=xmax, keep=keep, num_compute=num)
    axarr[x, y].plot(points[:, 0], points[:, 1], ',', color='k', alpha=0.8)
    axarr[x, y].set_xlim(xmin, xmax)
    axarr[x, y].set_ylim(ymin, ymax)
    axarr[x, y].set_title(r'${1} < \mu_{0} < {2}$, $2^{0}$ cycle'.format(index, xmin, xmax))
    axarr[x, y].set_yticks([])
    for i, mu in enumerate(mu_vals):
        axarr[x, y].plot(np.ones(10) * mu, np.linspace(0, 1, 10), 'r-', alpha=0.25)
        axarr[x, y].annotate(i + 1, xy=(mu, ymax-(0.05 * (ymax - ymin))), color='red')

fig, axarr = plt.subplots(3, 3, figsize=(12, 12))
plot_bifurcation(fig, axarr, 1, 0, 0, 2.9, 3.5, 0.4, 1, 500, 100, 1000)
plot_bifurcation(fig, axarr, 2, 0, 1, 3.4, 3.6, 0.8, 0.9, 500, 500, 5000)
plot_bifurcation(fig, axarr, 2, 0, 2, 3.53, 3.6, 0.8, 0.9, 500, 500, 5000)
plot_bifurcation(fig, axarr, 3, 1, 0, 3.56, 3.58, 0.888, 0.896, 500, 500, 5000)
plot_bifurcation(fig, axarr, 4, 1, 1, 3.568, 3.575, 0.889, 0.894, 1000, 500, 5000)
plot_bifurcation(fig, axarr, 5, 1, 2, 3.5695, 3.571, 0.890, 0.892, 1200, 500, 5000)
plot_bifurcation(fig, axarr, 6, 2, 0, 3.5698, 3.57, 0.8903, 0.8905, 1500, 500, 10000)
plot_bifurcation(fig, axarr, 7, 2, 1, 3.56992, 3.56997, 0.8903, 0.89045, 1500, 2000, 10000)
plot_bifurcation(fig, axarr, 8, 2, 2, 3.56994, 3.56995, 0.89041, 0.89043, 2000, 5000, 50000)
plt.tight_layout()
plt.show()

So what are these values that we've found? Let's record them.

In [11]:
points = bifurcation(xmin=1, xmax=4, precision=3000, num_compute=20000, keep=100)
plt.figure(figsize=(12, 4))
plt.plot(points[:, 0], points[:, 1], ',', color='k', alpha=0.8)
for mu in mu_vals:
    plt.plot(np.ones(10) * mu, np.linspace(0, 1, 10), 'r-', alpha=0.25)
plt.xlim(2.9, 4)
plt.show()

We can plot our cobweb diagrams with these more accurate values.

In [12]:
x = np.linspace(0, 1, 100)
fig, axarr = plt.subplots(3, 3, figsize=(12, 12))
for i in range(9):
    axarr[int(i / 3), i % 3].plot(x, x)
    mu = mu_vals[i]
    f = lambda x: mu * x * (1 - x)
    axarr[int(i / 3), i % 3].plot(x, f(x), alpha=1)
    web = cobweb(f, n=1000, start=800)
    axarr[int(i / 3), i % 3].plot(web[:, 0], web[:, 1], linewidth=0.5)
    axarr[int(i / 3), i % 3].set_title(r'$\mu = {}$'.format(mu))
plt.show()

We could keep recording these values if we wanted, as this will keep going infinitely, and appears to be approaching $\mu_\infty = 3.57$.

Feigenbaum Constant Estimation

To be more precise about this estimation, we can use the Feigenbaum constant, estimated with

$$ \delta_i = \frac{\mu_i - \mu_{i - 1}}{\mu_{i + 1} - \mu_i} $$
In [13]:
nice_table = pd.DataFrame({'0 mu':np.zeros(len(mu_vals)),
                           '1 delta':np.zeros(len(mu_vals))})
nice_table['0 mu'] = mu_vals
nice_table['1 delta'] = ([0, 0] +
                         [((mu_vals[i - 1] - mu_vals[i - 2]) / (mu_vals[i] - mu_vals[i - 1]))
                          for i in range(2, len(mu_vals))])
nice_table
Out[13]:
0 mu 1 delta
0 3.000000 0.000000
1 3.450000 0.000000
2 3.544000 4.787234
3 3.564500 4.585366
4 3.568750 4.823529
5 3.569700 4.473684
6 3.569892 4.935065
7 3.569934 4.638554
8 3.569943 4.530568
9 3.569945 4.490196
10 3.569946 4.573991

We note that these values are not that close to the real Feigenbaum constant, which is probably due to error with the bifurcation point locating process. In other words, I'm not that good at finding the exact point at which it bifurcates.

Period 3 and 5 Orbits

We found all powers of $2^n$, but we haven't found any odd-numbered orbits. Let's track these down.

In [109]:
%%cython

import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport exp

cdef float map_func(float mu, float x):
    return mu * x * (1 - x)

@cython.boundscheck(False) # turn off bounds-checking for entire function
@cython.wraparound(False)  # turn off negative index wrapping for entire function
def bifurcation(np.int64_t precision=1000, np.int64_t keep=500, np.int64_t num_compute=10000,
                np.float64_t xmin=0, np.float64_t xmax=4, np.float64_t ymin=0, np.float64_t ymax=1):
    """ Acquire bifurcation points for varying mu for map """
    cdef np.ndarray[np.float64_t, ndim=1] mu = np.linspace(xmin, xmax, precision, dtype=np.float64)
    cdef np.float64_t x = 0.5
    cdef np.int64_t i, j, k
    cdef np.ndarray[np.float64_t, ndim=2] points = np.zeros((len(mu) * keep, 2), dtype=np.float64)
    k = 0
    for i in range(len(mu)):
        for j in range(num_compute):
            x = map_func(mu[i], x)
            if j > (num_compute - keep): # we throw away the transient
                points[k, 0] = mu[i]
                points[k, 1] = x
                k += 1
    return points

Now we can perform the same process to find the odd-numbered orbits. I'm not doing 9 different levels of this, because they're hard to find, and it's tedious work.

In [301]:
mu_vals = np.array([3.84, 3.74, 3.702, 3.68725])
In [302]:
def plot_bifurcation(fig, axarr, index, x, y, xmin, xmax, ymin, ymax, precision, keep, num):
    points = bifurcation(precision=precision, xmin=xmin, xmax=xmax, keep=keep, num_compute=num)
    axarr[x, y].plot(points[:, 0], points[:, 1], ',', color='k', alpha=0.8)
    axarr[x, y].set_xlim(xmin, xmax)
    axarr[x, y].set_ylim(ymin, ymax)
    axarr[x, y].set_title(r'${1} < \mu_{0} < {2}$, ${3}$ cycle'.format(index, xmin, xmax, 2 * index + 1))
    axarr[x, y].set_yticks([])
    for i, mu in enumerate(mu_vals):
        axarr[x, y].plot(np.ones(10) * mu, np.linspace(0, 1, 10), 'r-', alpha=0.25)

fig, axarr = plt.subplots(2, 2, figsize=(8, 8))
plot_bifurcation(fig, axarr, 1, 0, 0, 3.8, 3.9, 0, 1, 500, 200, 2000)
plot_bifurcation(fig, axarr, 2, 0, 1, 3.735, 3.75, 0.1, 1, 500, 500, 5000)
plot_bifurcation(fig, axarr, 2, 1, 0, 3.7, 3.705, 0.2, 1, 500, 500, 5000)
plot_bifurcation(fig, axarr, 3, 1, 1, 3.687, 3.6875, 0.2, 1, 500, 300, 5000)
plt.tight_layout()
plt.show()