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. Parts of this notebook are inspired by Jake VanderPlas's FFT tutorial
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.
What does the period doubling look like in the frequency domain?¶
Image from Wikipedia
Image from source
Image from source
A low-dimensional model of period doubling in continuous time¶
Simulating a full turbulent flow in 3D typically requires simulating the Navier-Stokes equations, a set of partial differential equations that describe the motion of a fluid. Numerically, we often discretize partial differential equations on a grid and solve a system of coupled ordinary differential equations at each grid point. This is much higher dimensional than the dynamical systems that we have previously studied, which have all been in two or three dimensions.
Fortunately, the period-doubling route to chaos can be observed in much lower-dimensional systems, which capture the essential features of the transition to turbulence. One such system is the Rossler system, a three-dimensional continuous-time dynamical system. We can think of the Rossler coordinates as "bulk" measurements of the fluid flow, such as the average velocity across the domain, average pressure, and average acceleration.
Implementing a Rossler system in Python¶
Our implementation of the Rossler system will resemble previous dynamical systems that we have studied in this course. However, the Rossler system is a continuous-time system of the form:
$$\frac{d}{dt}\mathbf{x}(t) = \mathbf{f}(\mathbf{x}(t))$$
where $\mathbf{x}(t)$ is a three-dimensional and $\mathbf{f}(\mathbf{x})$ is a three-dimensional vector field. Because the Rossler system is a continuous-time system, we will need to use a numerical integrator to solve the system of ordinary differential equations. We will learn about integration in a different course module, but for now, we will use the solve_ivp function from the scipy.integrate module. This function takes a callable right-hand side function, an initial condition, and a time interval, and returns the solution to the system of ordinary differential equations.
# Import a built-in function that numerically solves ordinary differential equations
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 speed of an experimental fluid flow.
"""
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):
"""
This function takes a point x in phase space and returns the time derivative
vector at that point
Takes a 3 x 1 state vector and returns a 3 x 1 vector of time derivative value
"""
# x.shape = 3
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):
"""
This function numerically solves the differential equations for the Rossler
system given an initial condition and a time span.
"""
# for i in range(1000):
# x += dt * self.rhs(t, x)
return solve_ivp(self.rhs, t_span, x0, method="Radau", **kwargs)
## Define the initial conditions and time span over which we want to solve the system
ic = np.array([6.55, -0.317, 0.284])
t_span = [0, 500]
model = Rossler()
## Solve the system and plot the results
sol = model.solve(t_span, ic, max_step=0.01)
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,y,z')
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.8
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.8')
# 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.8')
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 sines of different frequencies, amplitudes, and phases:
$$ x(t) = \sum_{k=0}^{\infty} A_k \sin\left(2\pi k \frac{t}{T} + \phi_k\right) $$ Where $A_k \in \mathbb{R}$ is the amplitude of the $k^{th}$ sinusoid and $\phi_k \in \mathbb{R}$ is its phase of the $k^{th}$ sinusoid. We could alternatively remove the phase term and write the signal as a sum of sines and cosines (which would have their own amplitudes).
Technically, the Fourier transform only applies to periodic signals with a finite period $T$. However, when working with real time series, we often set $T$ equal to the longest "resolvable" period in the signal, which would be the length of the time series itself.
We can also write the signal in terms of complex exponentials: $$ x(t) = \sum_{k=0}^{\infty} a_k e^{2\pi i k \frac{t}{T}} $$ where we have absorbed the amplitude and phase into the complex amplitude $a_k = A_k e^{i\phi_k} \in \mathbb{C}$.
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 invertible integral transformation that maps a signal from the time domain to the frequency domain.
The Power Spectral Density (PSD) of a signal is the square of the amplitudes of the sinusoids that make up the signal:
$$ \text{PSD}(k) = |a_k|^2 $$
It ignores the phase information and simply tells us how much power is in each frequency bin.
basis_functions.shape
(2048, 2048)
## Load the period doubling time series
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')
# ## Create some example basis functions
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/micromamba/envs/cphy/lib/python3.11/site-packages/matplotlib/cbook.py:1699: ComplexWarning: Casting complex values to real discards the imaginary part return math.isfinite(val) /Users/william/micromamba/envs/cphy/lib/python3.11/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x2fb1d62d0>]
Implementing the discrete Fourier transform¶
A discrete approach to this algorithm involves projecting our signal onto a set of "pure" basis functions corresponding to single-frequency (monochromatic) signals. In continuous time, we would have a continuous and uncountable set of frequencies, but for a discrete signal we have a finite number of possible frequency bins spanning $1/N$ to $T/N$, where $T$ is the total time of the signal and $N$ is the number of timepoints.
If our observed signal is $x_n = x_1, x_2, ..., x_N$, then the discrete Fourier transform becomes $$ X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N} $$ with corresponding inverse transform $$ x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N} $$ In our naive implementation below, we compute the set of "pure" functions, and then project the signal onto this basis. This approach is easy to vectorize, since all operations that touch all $N$ elements of the signal can be written as dot products.
We can think of the DFT as a change of basis: we define a set of fixed trigonometric functions, and then project our signal onto these functions. For example, finding a single coefficient $X_3$ in the DFT requires a single dot product between a basis function $e^{-i~2\pi~3~n~/~N} \in \mathbb{C}^N$ and our signal $\mathbf{x} \in \mathbb{R}^N$:
$$ X_3 = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~3~n~/~N} = \mathbf{x} \cdot \mathbf{m}_3 $$
In order to implement the full DFT, we therefore want to define a set of all basis functions spanning our possible frequencies. We can do this by defining a matrix $M \in \mathbb{C}^{N \times N}$ where each row $M_k$ is the basis function $\mathbf{m}_k = e^{-i~2\pi~k~n~/~N}$.
$$ \mathbf{M} = \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & e^{-i~2\pi~1~/~N} & e^{-i~2\pi~2~/~N} & \cdots & e^{-i~2\pi~(N-1)~/~N} \\ 1 & e^{-i~2\pi~2~/~N} & e^{-i~2\pi~4~/~N} & \cdots & e^{-i~2\pi~2(N-1)~/~N} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-i~2\pi~(N-1)~/~N} & e^{-i~2\pi~2(N-1)~/~N} & \cdots & e^{-i~2\pi~(N-1)^2~/~N} \end{bmatrix} $$ where $ M_{kn} = e^{-i~2\pi~k~n~/~N}$
In matrix notation, the vectorized DFT therefore becomes $$ \mathbf{X} = M \cdot \mathbf{x} $$
# What does the M matrix look like?
# Python has built-in complex numbers, which numpy supports
N = 512
coeffs = np.exp(-2j * np.pi * np.arange(N)[:, None] * np.arange(N) / N)
print(coeffs.shape)
plt.figure(figsize=(9, 3))
# 3 subplots stacked vertically
plt.subplot(3, 1, 1)
plt.plot(np.real(coeffs)[1])
# plt.title("Second Row")
plt.xlim(0, 512); plt.xticks([])
plt.subplot(3, 1, 2)
plt.plot(np.real(coeffs)[64])
# plt.title("Seventh Row")
plt.xlim(0, 512); plt.xticks([])
plt.subplot(3, 1, 3)
plt.plot(np.real(coeffs)[128])
# plt.title("Highest Frequency Row")
plt.xlim(0, 512)
plt.xlabel("Frequency")
plt.figure(figsize=(9, 9))
plt.imshow(np.real(coeffs), cmap='viridis')
plt.xticks([]); plt.yticks([])
plt.figure(figsize=(9, 9))
plt.imshow(np.imag(coeffs), cmap='viridis')
plt.xticks([]); plt.yticks([])
(512, 512)
([], [])
Question:¶
Why does each frequency appear twice in the DFT M matrix?
plt.figure(figsize=(9, 9))
plt.imshow(np.imag(coeffs), cmap='viridis')
plt.xticks([]); plt.yticks([])
([], [])