```
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``.
= float(environ.get("LOGISTIC_X0_DEFAULT", "0.9"))
LOGISTIC_X0_DEFAULT = float(environ.get("LOGISTIC_STEP_SIZE", "0.01"))
LOGISTIC_STEP_SIZE = float(environ.get("LOGISTIC_STEP_SIZE_RR", "0.02"))
LOGISTIC_STEP_SIZE_RR = int(environ.get("LOGISTIC_ITERATIONS", "1000 "))
LOGISTIC_ITERATIONS = int(environ.get("LOGISTIC_ITERATIONS_KEPT", "100"))
LOGISTIC_ITERATIONS_KEPT
= pathlib.Path(".").resolve()
DIR = DIR / f"logistic-bifurcation-{LOGISTIC_ITERATIONS}-{LOGISTIC_STEP_SIZE}.csv"
LOGISTIC_DATA_BIFURCATION = DIR / f"cobweb-{LOGISTIC_STEP_SIZE}-{LOGISTIC_STEP_SIZE_RR}.json"
LOGISTIC_DATA_COBWEB_RR
def logistic_point(
float,
rr: float,
xx: *,
int = 1,
power: int = 0,
_power: -> float:
) """Compute the value of the logistic map for some coefficient `rr` and
number `xx`.
"""
+= 1
_power
= rr * xx * (1 - xx)
xx_next 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.
= xx, 0
yy, _power while _power < power:
= rr * yy * (1 - yy)
yy = _power + 1
_power
return yy
def logistic_array(
float,
rr: *,
float = 0.001,
step: int = 1,
power:
):"""Compute the graph of the logistic map on the closed interval $[0, 1]$."""
= np.arange(0, 1 + step, step=step)
xx = _logistic_array(rr, xx, power=power)
yy
return pd.DataFrame(dict(xx=xx, yy=yy))
def logistic_point_sequence(
float,
rr: float,
xx: *,
int = LOGISTIC_ITERATIONS,
iterations: -> Generator[float, None, None]:
) """Iterate the logistic sequence starting at `xx` for coefficient `rr` for
`iterations` number of steps."""
yield xx
= 0
iteration while (iteration := iteration + 1) < iterations:
yield (xx := logistic_point(rr, xx))
def logistic_point_attractor(
float,
rr: float,
xx: *,
int = LOGISTIC_ITERATIONS,
iterations: = LOGISTIC_ITERATIONS_KEPT,
iterations_kept -> List[float]:
) """Calculate the attractors for the logistic sequence at a point."""
= enumerate(logistic_point_sequence(rr, xx, iterations=iterations))
seq = iterations - iterations_kept
pos_min return [value for pos, value in seq if pos >= pos_min]
```

# The Logistic Map and Chaos

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:

## 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`

:

- Generate a random starting point
`x0`

between \(0\) and \(1\). - Do a large number of iterations.
- Choose some final number of terminating sequence elements to plot.

```
def create_bifurcation(
*,
float = float(2),
start: float = float(4),
stop: float = 0.05,
step: int = LOGISTIC_ITERATIONS,
iterations: int = LOGISTIC_ITERATIONS_KEPT,
iterations_kept:
):"""Returns a ``DataFrame`` containing coefficients, attractors, and initial
conditions.
"""
= np.arange(start, stop, step)
rs = list( np.random.rand() for _ in range(len(rs)))
x0s
return pd.DataFrame(
dict(
=rs,
rr=x0s,
x0=list(
attractors
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():
= create_bifurcation(step=0.0001)
data_raw
= [rr for rr in data_raw["rr"] for _ in range(LOGISTIC_ITERATIONS_KEPT)]
rr =list(itertools.chain(*data_raw["attractors"]) )
xx = pd.DataFrame(dict(rr=rr, xx=xx))
data
data.to_csv(LOGISTIC_DATA_BIFURCATION)else:
= pd.read_csv(LOGISTIC_DATA_BIFURCATION)
data
return data
if not (_ := DIR / "./bifurcation.png").exists():
= sns.scatterplot(data=load_bifurcation(), x="rr", y="xx", s=0.1)
plot_bifurcation
plot_bifurcation.figure.savefig(_) plt.close()
```

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:

- An initial input \(x_0\) is chosen. This can be written as the point \((x_0, x_0)\) on the graph.
- The output \(x_1 = \lambda_r(x_0)\) is computed. This can be represented as the point \((x_0, x_1)\).
- The next point will be \((x_1, x_1)\), directly on the line \(x = y\).
- 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\)*. - 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(
float,
rr: float = LOGISTIC_X0_DEFAULT,
x0: *,
int = LOGISTIC_ITERATIONS,
iterations:
):"""Iterates the elements of the cobweb sequence."""
= 0, x0
iteration, xx while (iteration := iteration + 1) < iterations:
yield (xx, xx)
yield (xx, xx := logistic_point(rr, xx))
def iter_cobweb_segments(
float,
rr: float = LOGISTIC_X0_DEFAULT,
x0: *,
int = LOGISTIC_ITERATIONS,
iterations:
):"""Iterates tuples representing the line segments making up the cobweb
diagram"""
= iter_cobweb(rr, x0, iterations=iterations)
points = zip(*points)
xx, yy 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(
float,
rr: float = LOGISTIC_X0_DEFAULT,
x0: *,
int,
frames: int = 200,
interval:
):"""Creates the cobweb diagram animation."""
# Do not fear, this just uses the values from the scope below when called.
def update_cobweb(frame):
= ax.plot(
line, 0],
data_cobweb[frame][1],
data_cobweb[frame][="red",
color
)return line
= logistic_array(rr)
data_logistic = list(iter_cobweb_segments(rr, x0, iterations=frames))
data_cobweb
0, 1)
plt.xlim(0, 1)
plt.ylim(= plt.subplots()
fig, ax
# Plot the logistic map.
ax.plot(:= data_logistic["xx"],
xx "yy"],
data_logistic[= "blue",
color
)="green")
ax.plot(xx, xx, color
# Animate cobweb.
= FuncAnimation(
animation
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():
= diagram_cobweb(3.8, 0.9, frames=200, interval=50)
animation ="pillow") animation.save(_, writer
```

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():
= diagram_cobweb(2, 0.9, frames=10)
animation ="pillow") animation.save(_, writer
```

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:
= json.load(file)
raw
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()
}
= np.arange(0, 4, LOGISTIC_STEP_SIZE_RR)
rrs = {
data_cobweb list(zip(*iter_cobweb(rr, x0, iterations=len(rrs))))
rr: 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$."""
= LOGISTIC_X0_DEFAULT
x0 = load_cobweb_rr(x0)
data_cobweb, data_logistic
= plt.subplots()
fig, ax
= np.array(sorted(data_cobweb))
rrs = min(rrs)
rr0
= ax.plot(
line_cobweb, 0],
data_cobweb[rr0][1],
data_cobweb[rr0][="blue",
color
)= ax.plot(
line_logistic, "xx"],
data_logistic[rr0]["yy"],
data_logistic[rr0][="red",
color
):= data_logistic[rr0]["xx"], xx, color="green")
ax.plot(xx = ax.text(0, 1.1, f"r = {rr0}")
tx
def update_cobweb_rr(frame_rr):
0])
line_cobweb.set_xdata(data_cobweb[frame_rr][1])
line_cobweb.set_ydata(data_cobweb[frame_rr][
"xx"])
line_logistic.set_xdata(data_logistic[frame_rr]["yy"])
line_logistic.set_ydata(data_logistic[frame_rr][
f"r = {round(frame_rr, 2)}")
tx.set_text(
return (line_cobweb, line_logistic)
= FuncAnimation(
animation
fig,
update_cobweb_rr,=rrs,
frames=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():
= diagram_cobweb_rr(interval=100)
animation ="pillow") animation.save(_, writer
```

## 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):
"yy"])
line.set_ydata(data_composition[frame][
f"k = {frame}")
tx.set_text(return (line,)
= [logistic_array(rr, power=k + 1) for k in range(frames)]
data_composition = plt.subplots()
fig, ax
= ax.text(0, 1.1, "k = 0")
tx = ax.plot(
line, := data_composition[0]["xx"],
xx 0]["yy"],
data_composition[= "blue",
color
)="red")
ax.plot(xx, xx, color= FuncAnimation(
animation
fig,
update_composition, =frames
frames
)
plt.close()return animation
```

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

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

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.

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

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

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