The Logistic Map and Chaos

python
math
The logistic map: Animations and properties.
Author

Adrian Cederberg

Published

August 1, 2024

Modified

August 27, 2024

Keywords

logistic map, logistic, map, chaos, python, seaborn, numpy, diagram, cobweb, bifurcation, polynomial

This is my first quarto blog post!

For the most part, I wanted to make this because the graphs are interesting to look at. Further, I wanted to test the abilities of quarto to use latex, markdown, and python.

The logistic map was used to explain chaos occurring in simple situations in a physics class I took back in college. This class went far deeper into the subject matter than this blog post does at the moment, but I might come back and write more in this blog post later.

For now, this will demonstrate some of the basic properties of the recurrence relation, provide some definitions, and showcase some cool animations.

The Logistic Map and Recurrence Relations

The logistic map is a function defined as

\[\begin{equation} l_r: x \mapsto rx(1-x) \end{equation}\]

where \(r \in [0, 4]\) and \(x \in [0, 1]\). A logistic sequence is defined as

\[\begin{equation} \lambda_r = \{l_r^k(x_0): k \in N\} \end{equation}\]

for some \(x_0\) and \(r\) and is defined by the recurrence relation

\[\begin{equation} x_{k+1} = rx_k(1-x_k) \end{equation}\]

In this case, the exponent on \(\lambda_r\) is the degree of composition; more precisely, \(\lambda^3_r(x) = \lambda_r(\lambda_r(\lambda_r(x)))\) for instance.

It is easy to write out some python functions that will compute the values of the logistic sequence:

import itertools
import json
import pathlib
from os import environ
from typing import Generator, List

import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from matplotlib.animation import FuncAnimation


# NOTE: Constants. These will determine how long it will take for this document
#       to render. Further, changing these documents can result in an error 
#       from deno with the message ``ERROR: RangeError: Maximum call stack size 
#       exceeded``.
LOGISTIC_X0_DEFAULT = float(environ.get("LOGISTIC_X0_DEFAULT", "0.9"))
LOGISTIC_STEP_SIZE = float(environ.get("LOGISTIC_STEP_SIZE", "0.01"))
LOGISTIC_STEP_SIZE_RR = float(environ.get("LOGISTIC_STEP_SIZE_RR", "0.02"))
LOGISTIC_ITERATIONS = int(environ.get("LOGISTIC_ITERATIONS", "1000 "))
LOGISTIC_ITERATIONS_KEPT = int(environ.get("LOGISTIC_ITERATIONS_KEPT", "100"))

DIR = pathlib.Path(".").resolve()
LOGISTIC_DATA_BIFURCATION = DIR / f"logistic-bifurcation-{LOGISTIC_ITERATIONS}-{LOGISTIC_STEP_SIZE}.csv"
LOGISTIC_DATA_COBWEB_RR = DIR / f"cobweb-{LOGISTIC_STEP_SIZE}-{LOGISTIC_STEP_SIZE_RR}.json"


def logistic_point(
  rr: float, 
  xx: float, 
  *, 
  power: int = 1,
  _power: int = 0,
) -> float:
  """Compute the value of the logistic map for some coefficient `rr` and 
  number `xx`.
  """

  _power += 1
  
  xx_next = rr * xx * (1 - xx) 
  if _power == power:
    return xx_next

  return logistic_point(rr, xx_next, power=power, _power=_power)


def _logistic_array(rr: float, xx, *, power = 1):
  """Helper for ``logistic_array``."""

  # NOTE: Numpy multiplies elements stepwise.

  yy, _power = xx, 0
  while _power < power:
    yy = rr * yy * (1 - yy)
    _power = _power + 1

  return yy


def logistic_array(
  rr: float,
  *,
  step: float = 0.001,
  power: int = 1,
):
  """Compute the graph of the logistic map on the closed interval $[0, 1]$."""

  xx = np.arange(0, 1 + step, step=step)
  yy = _logistic_array(rr, xx, power=power)

  return pd.DataFrame(dict(xx=xx, yy=yy))


def logistic_point_sequence(
  rr: float,
  xx: float,
  *,
  iterations: int = LOGISTIC_ITERATIONS,
) -> Generator[float, None, None]:
  """Iterate the logistic sequence starting at `xx` for coefficient `rr` for
  `iterations` number of steps."""

  yield xx 

  iteration = 0
  while (iteration := iteration + 1) < iterations:
    yield (xx := logistic_point(rr, xx))


def logistic_point_attractor(
  rr: float,
  xx: float,
  *,
  iterations: int = LOGISTIC_ITERATIONS,
  iterations_kept = LOGISTIC_ITERATIONS_KEPT,
) -> List[float]:
  """Calculate the attractors for the logistic sequence at a point."""

  seq = enumerate(logistic_point_sequence(rr, xx,  iterations=iterations)) 
  pos_min = iterations - iterations_kept 
  return [value for pos, value in seq if pos >= pos_min]

Bifurcation Diagram

Below is the bifurcation diagram for the logistic map. A bifurcation diagram is used to show the regions of attraction for the sequence depending on rr.

This graph is generated by iterating over the values of rr on the closed interval [0, 4] and completing the following process for each rr:

  1. Generate a random starting point x0 between \(0\) and \(1\).
  2. Do a large number of iterations.
  3. Choose some final number of terminating sequence elements to plot.
def create_bifurcation(
  *,
  start: float = float(2), 
  stop: float = float(4),
  step: float = 0.05,
  iterations: int = LOGISTIC_ITERATIONS,
  iterations_kept: int = LOGISTIC_ITERATIONS_KEPT,
):
  """Returns a ``DataFrame`` containing coefficients, attractors, and initial
  conditions.
  """

  rs = np.arange(start, stop, step)
  x0s =  list( np.random.rand() for _ in range(len(rs)))

  return pd.DataFrame(
    dict(
      rr=rs,
      x0=x0s,
      attractors=list(
        logistic_point_attractor(
          rr, 
          x0,
          iterations=iterations,
          iterations_kept=iterations_kept,
        )
        for rr, x0 in zip(rs, x0s)
      ) 
    )
  )


def load_bifurcation():
  """Load data if it has already been computed. This is especially useful when
  developing this document."""

  if not LOGISTIC_DATA_BIFURCATION.exists():
    data_raw = create_bifurcation(step=0.0001)

    rr = [rr for rr in data_raw["rr"] for _ in range(LOGISTIC_ITERATIONS_KEPT)] 
    xx =list(itertools.chain(*data_raw["attractors"]) )
    data = pd.DataFrame(dict(rr=rr, xx=xx))
    data.to_csv(LOGISTIC_DATA_BIFURCATION)
  else:
    data = pd.read_csv(LOGISTIC_DATA_BIFURCATION)

  return data


if not (_ := DIR / "./bifurcation.png").exists():
  plot_bifurcation = sns.scatterplot(data=load_bifurcation(), x="rr", y="xx", s=0.1)
  plot_bifurcation.figure.savefig(_)
  plt.close()

Logistic Bifurcation Diagram

A few well known facts about the logistic map are apparent from the figure:

  • Below a value of three, the logistic sequence actually converges
  • For \(r\) between \(3\) and roughly \(3.4\), the logistic map will have two distinct convergent subsequences.
  • For values of roughly \(r=3.4\) up until roughly \(3.6\) there will be four distinct convergent subsequences.
  • The bifurcation cascade blows up after approximately \(r=3.6\).

These facts can be proven by solving the recurrence relation. However doing so is not something I want to do in this particular blog post as the algebra is nasty and results in the solving for roots cubic equations or worse (as the result of solving a system of polynomial equations).

Cobweb Diagram

A cobweb diagram can be used to visualize the convergence (or lack there of) of a logistic sequence. The procedure for creating the diagram is as follows:

  1. An initial input \(x_0\) is chosen. This can be written as the point \((x_0, x_0)\) on the graph.
  2. The output \(x_1 = \lambda_r(x_0)\) is computed. This can be represented as the point \((x_0, x_1)\).
  3. The next point will be \((x_1, x_1)\), directly on the line \(x = y\).
  4. Since the next input into the logistic map is \(x_1\), we follow the \(x\) projection onto the logistic maps graph to get \((x_1, x_2 := lambda(x_1))\). This step is just step \(2\).
  5. Next step three is repeated for \(x_2\).

In the most explicit form this is

\[\begin{multline} \{ \\ (x_0, x_0), \\ (x_0, \lambda_r(x_0)), (\lambda_r(x_0), \lambda(x_0)), \\ (\lambda(x_0), \lambda^2(x_0)), (\lambda^2(x_0), \lambda^2(x_0))), \\ (\lambda^2(x_0), \lambda^3(x_0)), (\lambda^3(x_0), \lambda^3(x_0)), \\ ... \\ (\lambda^{k-1}(x_0), \lambda^{k}(x_0)), (\lambda^k(x_0), \lambda^k(x_0)), \\ ... \} \end{multline}\]

Moreover, if \(l_r = \{a_k := \lambda^k_r(x_0), k \in N\}\) is a logistic sequence, then we are just plotting

\[\begin{equation} \{ (a_0, a_0), (a_0, a_1), (a_1, a_1), (a_1, a_2), (a_2, a_2), ..., (a_k, a_{k+1}), (a_{k+1}, a_{k+1}), ... \} \end{equation}\]

If \(\lambda_r\) is plotted along with \(x \mapsto x\), the elements of the cobweb sequence should lie on either line, alternating between them. This visualization can be plotted using the following python code:

def iter_cobweb(
  rr: float,
  x0: float = LOGISTIC_X0_DEFAULT,
  *, 
  iterations: int = LOGISTIC_ITERATIONS,
):
  """Iterates the elements of the cobweb sequence."""
  iteration, xx = 0, x0
  while (iteration := iteration + 1) < iterations:
    yield (xx, xx)
    yield (xx, xx := logistic_point(rr, xx))


def iter_cobweb_segments(
  rr: float,
  x0: float = LOGISTIC_X0_DEFAULT,
  *,
  iterations: int = LOGISTIC_ITERATIONS,
):
  """Iterates tuples representing the line segments making up the cobweb 
  diagram"""

  points = iter_cobweb(rr, x0, iterations=iterations)
  xx, yy = zip(*points) 
  for xx_coords, yy_coords in zip(
      zip(xx[:-1], xx[1:]),
      zip(yy[:-1], yy[1:])
  ):
    yield (xx_coords, yy_coords)


def diagram_cobweb(
  rr: float, 
  x0: float = LOGISTIC_X0_DEFAULT, 
  *, 
  frames: int,
  interval: int = 200,
):
  """Creates the cobweb diagram animation."""

  # Do not fear, this just uses the values from the scope below when called.
  def update_cobweb(frame):
    line, = ax.plot(
      data_cobweb[frame][0],
      data_cobweb[frame][1],
      color="red",
    )
    return line


  data_logistic = logistic_array(rr)
  data_cobweb = list(iter_cobweb_segments(rr, x0, iterations=frames))

  plt.xlim(0, 1)
  plt.ylim(0, 1)
  fig, ax = plt.subplots()

  # Plot the logistic map.
  ax.plot(
    xx := data_logistic["xx"],
    data_logistic["yy"],
    color = "blue",
  )
  ax.plot(xx, xx, color="green")

  # Animate cobweb.
  animation = FuncAnimation(
    fig, 
    update_cobweb, 
    frames=frames,
    interval=interval
  )

  # Why do I have to do this twice?
  plt.close()
  plt.close()
  return animation

Starting with a large value we can see that sequence probably does not converge

if not (_ := DIR / "cobweb-3.8-0.9.gif" ).exists():
  animation = diagram_cobweb(3.8, 0.9, frames=200, interval=50)
  animation.save(_, writer="pillow")

Logistic Cobweb Diagram for \(r = 3.8\).

The lack of convergence for large \(r\) values will also be demonstrated in the section about the bifurcation diagram. With a smaller value of \(r\) the convergence easy to confirm visually:

if not (_ := DIR / "cobweb-2-0.9.gif" ).exists():
  animation = diagram_cobweb(2, 0.9, frames=10)
  animation.save(_, writer="pillow")

Logistic Cobweb Diagram for \(r = 2\).

The following code can be used to animate the sequence starting at the same point for varying values of \(r\). This will show how convergence varies as the value of \(r\) varies.

def load_cobweb_rr(x0):
  """Load the data if it exists already instead of generating it."""

  if LOGISTIC_DATA_COBWEB_RR.exists():
    with open(LOGISTIC_DATA_COBWEB_RR, "r") as file:
      raw = json.load(file)

    return {
      float(rr): df_raw
      for rr, df_raw in     raw["data_cobweb"].items()
    }, {
      float(rr): pd.DataFrame(json.loads(df_raw))
      for rr, df_raw in raw["data_logistic"].items()
    }

  rrs = np.arange(0, 4, LOGISTIC_STEP_SIZE_RR)
  data_cobweb = {
    rr: list(zip(*iter_cobweb(rr, x0, iterations=len(rrs))))
    for rr in rrs
  }
  data_logistic = {
    rr: logistic_array(rr)
    for rr in rrs
  }

  with open(LOGISTIC_DATA_COBWEB_RR, "w") as file:
    json.dump(
      {
        "data_cobweb": data_cobweb,
        "data_logistic": {
          rr: df.to_json() 
          for rr, df in data_logistic.items()
        }
      }, 
      file
    )

  return (data_cobweb, data_logistic)


def diagram_cobweb_rr(interval: int = 200):
  """Create the cobweb diagram animation for a fixed starting point and a 
  varied $r$."""

  x0 = LOGISTIC_X0_DEFAULT
  data_cobweb, data_logistic = load_cobweb_rr(x0)

  fig, ax = plt.subplots()

  rrs = np.array(sorted(data_cobweb))
  rr0 = min(rrs)


  line_cobweb, = ax.plot(
    data_cobweb[rr0][0],
    data_cobweb[rr0][1],
    color="blue",
  )
  line_logistic, = ax.plot(
    data_logistic[rr0]["xx"],
    data_logistic[rr0]["yy"],
    color="red",
  )
  ax.plot(xx := data_logistic[rr0]["xx"], xx, color="green")
  tx = ax.text(0, 1.1, f"r = {rr0}")

  def update_cobweb_rr(frame_rr):
    line_cobweb.set_xdata(data_cobweb[frame_rr][0])
    line_cobweb.set_ydata(data_cobweb[frame_rr][1])

    line_logistic.set_xdata(data_logistic[frame_rr]["xx"])
    line_logistic.set_ydata(data_logistic[frame_rr]["yy"])

    tx.set_text(f"r = {round(frame_rr, 2)}")

    return (line_cobweb, line_logistic)

  animation = FuncAnimation(
    fig,
    update_cobweb_rr,
    frames=rrs,
    interval=interval,
  )
  plt.close()
  return animation

The animation below shows cobweb diagrams with the value of \(r\) varied along \([0, 4]\).

if not (_ := DIR / "cobweb-rr-0.9.gif").exists():
  animation = diagram_cobweb_rr(interval=100)
  animation.save(_, writer="pillow")

Cobweb Diagram with varying \(r\) values and \(x=0.9\)

Logistic Compositions

In this context, logistic compositions refers to the sequence of functions

\[\begin{equation} \{\lambda_r^k: k \in N\} \end{equation}\]

for some \(r\). This sequence exists in the space of continuously differentiable functions on the closed unit interval, which is notated as \(C^\infty([0, 1])\). Here I will not bother to proof that polynomials are continuously differentiable as this is easy to find in any analysis book or online. In fact the only point of mentioning that this sequence exists in this space (in this current version of the blog post) is drive home that this is a sequence of functions, not numbers.

It is possible to make some nice graphs of the logistic function and its composition with itself. The following code is used to generate the subsequent diagrams:

def diagram_composition(rr: float, frames: int):
  """Create the composition animation for some ``rr``."""

  def update_composition(frame):
    line.set_ydata(data_composition[frame]["yy"])

    tx.set_text(f"k = {frame}")
    return (line,)

  data_composition = [logistic_array(rr, power=k + 1) for k in range(frames)]
  fig, ax = plt.subplots()

  tx = ax.text(0, 1.1, "k = 0")
  line, = ax.plot(
    xx := data_composition[0]["xx"],
    data_composition[0]["yy"],
    color = "blue",
  )
  ax.plot(xx, xx, color="red")
  animation = FuncAnimation(
    fig, 
    update_composition, 
    frames=frames
  )

  plt.close()
  return animation

The plot below shows \(\lambda_4^k\) (the compositions):

animation = diagram_composition(4, frames=64)
if not (output_path := DIR / "composition-4.gif").exists():
  animation.save(output_path, writer="pillow")

Compositions for \(r=4\)

Seeing the evolutions of the above plot it is clear why the logistic sequence does not converge for higher values of \(r\). In the case that \(r = 3\) the sequence does converge (as we will see in the next section) and the compositions have maxima and minima that are in a smaller range.

animation = diagram_composition(3, frames=9)
if not (output_path := DIR / "composition-3.gif").exists():
  animation.save(output_path, writer="pillow")

Compositions for \(r=3\)

Finally, when \(r=2\) the compositions appears to only have \(3\) extrema no matter the degree of composition:

animation = diagram_composition(1, frames=11)
if not (output_path := DIR / "composition-1.gif").exists():
  animation.save(output_path, writer="pillow")

Compositions for \(r=1\)