# Detecting chaotic dynamics with the fast Fourier transformÂ¶

Preamble: Run the cells below to import the necessary Python packages

*This notebook created by William Gilpin. Consult the course website for all content and GitHub repository for raw files and runnable online code.*

```
import numpy as np
# Wipe all outputs from this notebook
from IPython.display import Image, clear_output, display
clear_output(True)
# Import local plotting functions and in-notebook display functions
import matplotlib.pyplot as plt
%matplotlib inline
```

# Revisiting period doubling and the route to chaosÂ¶

In an earlier notebook, we analyzed the logistic map, a deterministic discrete-time dynamical system:

$$x_{n+1} = f(x_n) = r x_n (1 - x_n)$$

where $x_n$ is the population at time $n$ and $r$ is a parameter that controls the growth rate of the population in each generation. Recall that, as $r$ increases, the system undergoes a series of bifurcations, from a stable fixed point to a period-2 cycle to a period-4 cycle, period-8, and so on, until it reaches a critical value of $r$ at which the system becomes chaotic.

This mechanism is observed in many other physical systems, including, famously, in the velocity fluctuations of a rotating fluid flow as the rotation speed is increased. The flow fluctuations gradually double in frequency, then double again, and so on, until the flow becomes fully turbulent.

*This notebook was partly adapted from Jake VanderPlas's excellent FFT tutorial*

### What does the period doubling look like in the frequency domain?**Â¶

*Image from Wikipedia*

*Image from source*

*Image from source*

```
```

## Period doubling in continuous timeÂ¶

Rather than simulating a full turbulent flow, we will simulate a continuous-time deterministic dynamical system that exhibits period-doubling bifurcations leading to chaos. This three-dimensional system is much simpler than an infinite-dimensional simulation of the Navier-Stokes equations, but it still exhibits the same qualitative behavior.

```
# solve_ivp
from scipy.integrate import solve_ivp
class Rossler:
"""
A continuous-time dynamical system that exhibits some properties associated with
velocity fluctuations in turbulent fluid flow. The parameter c can b thought of as
increasing the rotary driving force on the sequence.
This system exhibits period doubling in continuous time
"""
def __init__(self, a=0.2, b=0.2, c=5.7):
self.a = a
self.b = b
self.c = c
def rhs(self, t, x):
return np.array([-x[1] - x[2],
x[0] + self.a * x[1],
self.b + x[2] * (x[0] - self.c)])
def solve(self, t_span, x0, **kwargs):
return solve_ivp(self.rhs, t_span, x0, method="Radau", **kwargs)
ic = np.array([6.55, -0.317, 0.284])
t_span = [0, 500]
model = Rossler()
sol = model.solve(t_span, ic, max_step=0.01)
# discard a transient
plt.figure(figsize=(9, 5))
plt.plot(sol.t, sol.y.T)
plt.xlabel('t')
plt.ylabel('x')
plt.figure(figsize=(5, 5))
plt.plot(sol.y[0], sol.y[1])
plt.xlabel('x')
plt.ylabel('y')
```

Text(0, 0.5, 'y')

We will vary the value of a parameter $c$ in the system in order to observe period doubling leading to full chaos, which is turbulence.

```
ic = np.array([6.55, -0.317, 0.284])
t_span = [0, 500]
model = Rossler()
model.c = 2.0
sol = model.solve(t_span, ic, max_step=0.1)
y2 = sol.y.T[-2048:]
model.c = 3.8
sol = model.solve(t_span, ic, max_step=0.1)
y4 = sol.y.T[-2048:]
model.c = 4.1
sol = model.solve(t_span, ic, max_step=0.1)
y8 = sol.y.T[-2048:]
model.c = 5.7
sol = model.solve(t_span, ic, max_step=0.1)
yc = sol.y.T[-2048:]
t = sol.t[:2048]
plt.figure(figsize=(9, 4))
plt.subplot(4, 1, 1)
plt.plot(t, y2[:, 0])
plt.xlim(0, 200)
plt.subplot(4, 1, 2)
plt.plot(t, y4[:, 0])
plt.xlim(0, 200)
plt.subplot(4, 1, 3)
plt.plot(t, y8[:, 0])
plt.xlim(0, 200)
plt.subplot(4, 1, 4)
plt.plot(t, yc[:, 0])
plt.xlim(0, 200)
plt.xlabel('t')
plt.ylabel('x')
plt.figure(figsize=(9, 9))
plt.subplot(2, 2, 1)
plt.plot(y2[:, 0], y2[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('c = 2.0')
plt.subplot(2, 2, 2)
plt.plot(y4[:, 0], y4[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('c = 4.0')
plt.subplot(2, 2, 3)
plt.plot(y8[:, 0], y8[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('c = 4.35')
plt.subplot(2, 2, 4)
plt.plot(yc[:, 0], yc[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('c = 5.7')
# plt.figure(figsize=(9, 5))
# plt.plot(sol.t, sol.y.T)
# plt.xlabel('t')
# plt.ylabel('x')
```

Text(0.5, 1.0, 'c = 5.7')

```
```

# The discrete Fourier TransformÂ¶

Suppose that we want to analyze a signal $x(t)$ that is a function of time. We can think of this signal as a sum of sinusoids of different frequencies, amplitudes, and phases:

$$ x(t) = \sum_{k=0}^{\infty} A_k \sin(2\pi f_k t + \phi_k) $$

where $A_k$ is the amplitude of the $k$th sinusoid, $f_k$ is its frequency, and $\phi_k$ is its phase. The discrete Fourier transform (DFT) is an algorithm that allows us to compute the amplitudes, frequencies, and phases of the sinusoids that make up a signal. It is an information-preserving integral transformation that maps a signal from the time domain to the frequency domain.

The power spectrum of a signal is the square of the amplitudes of the sinusoids that make up the signal. We can think of it as a histogram showing the relative amount of power in each frequency bin.

```
time_seriesc = yc[:, 0]
time_series2 = y2[:, 0]
time_series4 = y4[:, 0]
time_series8 = y8[:, 0]
plt.figure(figsize=(8, 4))
plt.plot(time_seriesc, 'k')
n = len(time_seriesc)
basis_functions = np.exp(-2j * np.pi * np.arange(n)[:, None] * np.arange(n) / n)
plt.figure(figsize=(8, 4))
plt.plot(basis_functions[3])
plt.figure(figsize=(8, 4))
plt.plot(basis_functions[5])
plt.figure(figsize=(8, 4))
plt.plot(basis_functions[7])
```

/Users/william/miniconda3/envs/cphy/lib/python3.10/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)

[<matplotlib.lines.Line2D at 0x16b666680>]