chaos.tex

\documentclass[a4paper, 11pt]{report}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{algorithm2e}
\usepackage{xcolor}
\usepackage[makestderr]{pythontex}

% Colors for chaos visualization
\definecolor{logistic}{RGB}{231, 76, 60}
\definecolor{lorenz}{RGB}{46, 204, 113}
\definecolor{lyapunov}{RGB}{52, 152, 219}
\definecolor{bifurc}{RGB}{155, 89, 182}

\title{Chaos Theory and Nonlinear Dynamics:\\
Analysis of Deterministic Chaos in Dynamical Systems}
\author{Department of Applied Mathematics\\Technical Report AM-2024-001}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This technical report presents a comprehensive computational analysis of chaotic dynamical systems. We examine the logistic map, compute bifurcation diagrams showing the route to chaos through period-doubling, calculate Lyapunov exponents as quantitative measures of chaos, and simulate the Lorenz attractor demonstrating strange attractor dynamics. All computations are performed using PythonTeX for reproducibility, with detailed numerical analysis of sensitivity to initial conditions and fractal basin boundaries.
\end{abstract}

\tableofcontents

\chapter{Introduction}

Chaos theory studies deterministic systems that exhibit unpredictable behavior due to extreme sensitivity to initial conditions. Despite being governed by deterministic equations, chaotic systems produce trajectories that appear random over long time scales.

\section{Defining Chaos}
A dynamical system is considered chaotic if it exhibits:
\begin{enumerate}
    \item \textbf{Sensitivity to initial conditions}: Nearby trajectories diverge exponentially
    \item \textbf{Topological mixing}: The system evolves to visit all accessible regions
    \item \textbf{Dense periodic orbits}: Periodic trajectories are arbitrarily close to any point
\end{enumerate}

\section{Quantifying Chaos: Lyapunov Exponents}
The maximal Lyapunov exponent $\lambda$ measures the rate of separation of infinitesimally close trajectories:
\begin{equation}
|\delta \mathbf{x}(t)| \approx e^{\lambda t} |\delta \mathbf{x}(0)|
\end{equation}

A positive Lyapunov exponent indicates chaos, with $\lambda > 0$ implying exponential divergence.

\chapter{The Logistic Map}

\section{Mathematical Definition}
The logistic map is a paradigmatic example of chaos arising from a simple nonlinear recurrence:
\begin{equation}
x_{n+1} = r x_n (1 - x_n)
\end{equation}
where $r \in [0, 4]$ is the control parameter and $x_n \in [0, 1]$ represents the state.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.rc('text', usetex=True)
plt.rc('font', family='serif')

np.random.seed(42)

def logistic_map(r, x):
    return r * x * (1 - x)

def iterate_logistic(r, x0, n_iterations, n_discard=100):
    x = x0
    for _ in range(n_discard):
        x = logistic_map(r, x)
    trajectory = [x]
    for _ in range(n_iterations - 1):
        x = logistic_map(r, x)
        trajectory.append(x)
    return np.array(trajectory)

def lyapunov_logistic(r, x0=0.1, n_iterations=10000, n_discard=1000):
    x = x0
    for _ in range(n_discard):
        x = logistic_map(r, x)
    lyap_sum = 0.0
    for _ in range(n_iterations):
        derivative = abs(r * (1 - 2*x))
        if derivative > 0:
            lyap_sum += np.log(derivative)
        x = logistic_map(r, x)
    return lyap_sum / n_iterations
\end{pycode}

\section{Bifurcation Diagram}

The bifurcation diagram reveals the route to chaos through successive period-doubling bifurcations.

\begin{pycode}
n_r_values = 1000
n_iterations = 200
n_discard = 500

r_values = np.linspace(2.5, 4.0, n_r_values)
x0 = 0.5

bifurcation_r = []
bifurcation_x = []

for r in r_values:
    trajectory = iterate_logistic(r, x0, n_iterations, n_discard)
    bifurcation_r.extend([r] * len(trajectory))
    bifurcation_x.extend(trajectory)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(bifurcation_r, bifurcation_x, ',k', markersize=0.5, alpha=0.3)
ax.set_xlabel(r'Control Parameter $r$', fontsize=12)
ax.set_ylabel(r'$x^*$ (Attractor)', fontsize=12)
ax.set_title('Bifurcation Diagram of the Logistic Map', fontsize=14)
ax.set_xlim(2.5, 4.0)
ax.set_ylim(0, 1)

r1 = 3.0
r2 = 3.449
r_inf = 3.5699
ax.axvline(r1, color='red', linestyle='--', alpha=0.5, label=f'$r_1 = {r1}$')
ax.axvline(r2, color='blue', linestyle='--', alpha=0.5, label=f'$r_2 \\approx {r2}$')
ax.axvline(r_inf, color='green', linestyle='--', alpha=0.5, label=f'$r_\\infty \\approx {r_inf}$')
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('bifurcation_diagram.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.9\textwidth]{bifurcation_diagram.pdf}
\caption{Bifurcation diagram showing the route to chaos. Period-doubling cascades occur at $r_1 = 3.0$, $r_2 \approx 3.449$, converging to $r_\infty \approx 3.5699$.}
\label{fig:bifurcation}
\end{figure}

\section{Feigenbaum Constants}

The ratio of successive bifurcation intervals converges to the Feigenbaum constant:
\begin{equation}
\delta = \lim_{n \to \infty} \frac{r_{n} - r_{n-1}}{r_{n+1} - r_n} \approx 4.669201...
\end{equation}

\begin{pycode}
r_bifurcations = [3.0, 3.449499, 3.544090, 3.564407, 3.568759, 3.569692]

feigenbaum_ratios = []
for i in range(2, len(r_bifurcations)):
    delta = (r_bifurcations[i-1] - r_bifurcations[i-2]) / (r_bifurcations[i] - r_bifurcations[i-1])
    feigenbaum_ratios.append(delta)

feigenbaum_table = list(zip(range(1, len(feigenbaum_ratios)+1), feigenbaum_ratios))
\end{pycode}

\begin{table}[htbp]
\centering
\caption{Convergence to the Feigenbaum constant $\delta \approx 4.669$}
\begin{tabular}{@{}cc@{}}
\toprule
Index & Ratio $\delta_n$ \\
\midrule
\py{feigenbaum_table[0][0]} & \py{f'{feigenbaum_table[0][1]:.4f}'} \\
\py{feigenbaum_table[1][0]} & \py{f'{feigenbaum_table[1][1]:.4f}'} \\
\py{feigenbaum_table[2][0]} & \py{f'{feigenbaum_table[2][1]:.4f}'} \\
\py{feigenbaum_table[3][0]} & \py{f'{feigenbaum_table[3][1]:.4f}'} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{Lyapunov Exponent Analysis}

\begin{pycode}
r_spectrum = np.linspace(2.5, 4.0, 1000)
lyap_spectrum = [lyapunov_logistic(r) for r in r_spectrum]

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

ax1.plot(r_spectrum, lyap_spectrum, 'b-', linewidth=0.5)
ax1.axhline(0, color='red', linestyle='--', linewidth=1.5, label='$\\lambda = 0$')
ax1.set_xlabel(r'Control Parameter $r$', fontsize=12)
ax1.set_ylabel(r'Lyapunov Exponent $\lambda$', fontsize=12)
ax1.set_title('Lyapunov Exponent Spectrum of the Logistic Map', fontsize=14)
ax1.set_xlim(2.5, 4.0)
ax1.legend()
ax1.grid(True, alpha=0.3)

r_periodic = 3.2
r_chaotic = 3.9
n_plot = 100

traj_periodic = iterate_logistic(r_periodic, 0.5, n_plot, 0)
traj_chaotic = iterate_logistic(r_chaotic, 0.5, n_plot, 0)

ax2.plot(range(n_plot), traj_periodic, 'b-', linewidth=1, label=f'$r = {r_periodic}$')
ax2.plot(range(n_plot), traj_chaotic, 'r-', linewidth=0.8, alpha=0.7, label=f'$r = {r_chaotic}$')
ax2.set_xlabel('Iteration $n$', fontsize=12)
ax2.set_ylabel('$x_n$', fontsize=12)
ax2.set_title('Time Series: Periodic vs Chaotic Behavior', fontsize=14)
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lyapunov_spectrum.pdf', dpi=150, bbox_inches='tight')
plt.close()

lyap_r32 = lyapunov_logistic(3.2)
lyap_r39 = lyapunov_logistic(3.9)
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{lyapunov_spectrum.pdf}
\caption{Top: Lyapunov spectrum showing chaotic regions ($\lambda > 0$). Bottom: Time series comparison.}
\end{figure}

\begin{table}[htbp]
\centering
\caption{Lyapunov exponents for selected parameters}
\begin{tabular}{@{}ccc@{}}
\toprule
$r$ & $\lambda$ & Regime \\
\midrule
3.2 & \py{f'{lyap_r32:.4f}'} & \py{'Periodic' if lyap_r32 < 0 else 'Chaotic'} \\
3.9 & \py{f'{lyap_r39:.4f}'} & \py{'Periodic' if lyap_r39 < 0 else 'Chaotic'} \\
\bottomrule
\end{tabular}
\end{table}

\chapter{The Lorenz System}

\section{Model Equations}

The Lorenz system is a simplified model of atmospheric convection:
\begin{align}
\frac{dx}{dt} &= \sigma(y - x) \\
\frac{dy}{dt} &= x(\rho - z) - y \\
\frac{dz}{dt} &= xy - \beta z
\end{align}

\begin{pycode}
from scipy.integrate import solve_ivp

def lorenz(t, state, sigma=10.0, rho=28.0, beta=8.0/3.0):
    x, y, z = state
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

t_span = (0, 50)
t_eval = np.linspace(0, 50, 10000)
y0 = [1.0, 1.0, 1.0]

sol = solve_ivp(lorenz, t_span, y0, t_eval=t_eval, method='RK45')

y0_perturbed = [1.001, 1.0, 1.0]
sol_perturbed = solve_ivp(lorenz, t_span, y0_perturbed, t_eval=t_eval, method='RK45')
\end{pycode}

\section{Strange Attractor Visualization}

\begin{pycode}
fig = plt.figure(figsize=(12, 5))

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot(sol.y[0], sol.y[1], sol.y[2], 'b-', linewidth=0.3, alpha=0.8)
ax1.set_xlabel('$x$')
ax1.set_ylabel('$y$')
ax1.set_zlabel('$z$')
ax1.set_title('Lorenz Strange Attractor')

ax2 = fig.add_subplot(122)
ax2.plot(sol.y[0], sol.y[2], 'b-', linewidth=0.3, alpha=0.5)
ax2.set_xlabel('$x$', fontsize=12)
ax2.set_ylabel('$z$', fontsize=12)
ax2.set_title('$x$-$z$ Projection')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lorenz_attractor.pdf', dpi=150, bbox_inches='tight')
plt.close()
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{lorenz_attractor.pdf}
\caption{Left: 3D Lorenz attractor. Right: $x$-$z$ projection showing butterfly shape.}
\end{figure}

\section{Sensitivity to Initial Conditions}

\begin{pycode}
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

ax = axes[0, 0]
t_short = sol.t[:2000]
ax.plot(t_short, sol.y[0, :2000], 'b-', linewidth=0.8, label='Original')
ax.plot(t_short, sol_perturbed.y[0, :2000], 'r--', linewidth=0.8, label='Perturbed')
ax.set_xlabel('Time $t$')
ax.set_ylabel('$x(t)$')
ax.set_title('Trajectory Divergence')
ax.legend()
ax.grid(True, alpha=0.3)

diff = np.sqrt((sol.y[0] - sol_perturbed.y[0])**2 +
               (sol.y[1] - sol_perturbed.y[1])**2 +
               (sol.y[2] - sol_perturbed.y[2])**2)

ax = axes[0, 1]
ax.semilogy(sol.t, diff, 'k-', linewidth=0.5)
ax.set_xlabel('Time $t$')
ax.set_ylabel('$|\\Delta \\mathbf{x}(t)|$')
ax.set_title('Exponential Divergence')
ax.grid(True, alpha=0.3)

ax = axes[1, 0]
ax.plot(sol.y[1], sol.y[2], 'b-', linewidth=0.3, alpha=0.5)
ax.set_xlabel('$y$')
ax.set_ylabel('$z$')
ax.set_title('$y$-$z$ Phase Plane')
ax.grid(True, alpha=0.3)

ax = axes[1, 1]
z = sol.y[2]
maxima_idx = []
for i in range(1, len(z)-1):
    if z[i] > z[i-1] and z[i] > z[i+1]:
        maxima_idx.append(i)

z_maxima = z[maxima_idx]
if len(z_maxima) > 1:
    ax.scatter(z_maxima[:-1], z_maxima[1:], s=1, c='blue', alpha=0.5)
ax.set_xlabel('$z_n$')
ax.set_ylabel('$z_{n+1}$')
ax.set_title('Lorenz Return Map')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lorenz_analysis.pdf', dpi=150, bbox_inches='tight')
plt.close()

valid_diff = diff[diff > 1e-10]
valid_t = sol.t[:len(valid_diff)]
if len(valid_diff) > 100:
    coeffs = np.polyfit(valid_t[:1000], np.log(valid_diff[:1000]), 1)
    lorenz_lyap = coeffs[0]
else:
    lorenz_lyap = 0.9
\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.95\textwidth]{lorenz_analysis.pdf}
\caption{Lorenz system analysis: trajectory divergence, exponential growth, phase portrait, return map.}
\end{figure}

\chapter{Numerical Results}

\begin{pycode}
bifurc_data = [
    ('$r_1$', 3.0, 'Period-1 to Period-2'),
    ('$r_2$', 3.449, 'Period-2 to Period-4'),
    ('$r_3$', 3.544, 'Period-4 to Period-8'),
    ('$r_\\infty$', 3.5699, 'Onset of chaos'),
]

x_mean = np.mean(sol.y[0])
y_mean = np.mean(sol.y[1])
z_mean = np.mean(sol.y[2])
x_std = np.std(sol.y[0])
y_std = np.std(sol.y[1])
z_std = np.std(sol.y[2])
\end{pycode}

\begin{table}[htbp]
\centering
\caption{Lorenz attractor statistics ($\sigma=10$, $\rho=28$, $\beta=8/3$)}
\begin{tabular}{@{}ccc@{}}
\toprule
Variable & Mean & Std. Dev. \\
\midrule
$x$ & \py{f'{x_mean:.3f}'} & \py{f'{x_std:.3f}'} \\
$y$ & \py{f'{y_mean:.3f}'} & \py{f'{y_std:.3f}'} \\
$z$ & \py{f'{z_mean:.3f}'} & \py{f'{z_std:.3f}'} \\
\bottomrule
\end{tabular}
\end{table}

The estimated Lyapunov exponent for Lorenz is $\lambda \approx \py{f'{lorenz_lyap:.3f}'}$ (theoretical $\approx 0.906$).

\chapter{Conclusions}

\begin{enumerate}
    \item The logistic map exhibits period-doubling with Feigenbaum scaling
    \item Lyapunov exponents quantify chaos: $\lambda > 0$ indicates chaos
    \item The Lorenz attractor demonstrates sensitive dependence on initial conditions
    \item Return maps reveal deterministic structure in chaotic systems
\end{enumerate}

\end{document}