Dynamical Systems

Dynamical Systems Analysis

Overview

Qualitative and quantitative analysis of nonlinear dynamical systems. Phase portraits, fixed point classification, stability analysis, bifurcation diagrams, Poincare sections, and Lyapunov exponent computation.

When to Use

• Visualizing flow in phase space (2D and 3D systems)
• Finding and classifying fixed points (stable/unstable nodes, spirals, saddles, centers)
• Bifurcation analysis (how qualitative behavior changes with parameters)
• Detecting chaos (Lyapunov exponents, sensitivity to initial conditions)
• Poincare sections for periodicity analysis
• Limit cycle detection and characterization

Core Workflows

1. Phase Portrait (2D System)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def system(t, y, mu=1.0):
    """Van der Pol oscillator"""
    x, v = y
    return [v, mu * (1 - x**2) * v - x]

# Vector field on a grid
x_range = np.linspace(-4, 4, 20)
v_range = np.linspace(-6, 6, 20)
X, V = np.meshgrid(x_range, v_range)
U = V
W = 1.0 * (1 - X**2) * V - X

fig, ax = plt.subplots(figsize=(10, 8))

# Streamlines
ax.streamplot(X, V, U, W, density=1.5, color='gray', linewidth=0.5, arrowsize=1)

# Sample trajectories from different ICs
colors = plt.cm.viridis(np.linspace(0, 1, 6))
for i, ic in enumerate([[0.1, 0], [3, 0], [0, 4], [-2, -3], [1, -5], [4, 4]]):
    sol = solve_ivp(system, (0, 30), ic, t_eval=np.linspace(0, 30, 3000),
                    rtol=1e-10, atol=1e-12)
    ax.plot(sol.y[0], sol.y[1], color=colors[i], linewidth=1.2)
    ax.plot(ic[0], ic[1], 'o', color=colors[i], markersize=6)

# Fixed points
ax.plot(0, 0, 'rx', markersize=12, markeredgewidth=3, label='Unstable fixed point')

ax.set_xlabel('x', fontsize=13)
ax.set_ylabel('dx/dt', fontsize=13)
ax.set_title('Van der Pol Oscillator Phase Portrait (μ=1)', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.savefig('phase_portrait.png', dpi=150, bbox_inches='tight')

2. Fixed Point Analysis

from scipy.optimize import fsolve

def rhs(y, mu=1.0):
    """Autonomous system: dy/dt = f(y)"""
    x, v = y
    return [v, mu * (1 - x**2) * v - x]

def jacobian(y, mu=1.0):
    """Jacobian matrix J_ij = df_i/dy_j"""
    x, v = y
    return np.array([
        [0, 1],
        [-2*mu*x*v - 1, mu*(1 - x**2)]
    ])

# Find fixed points
y_guess = [0, 0]
fp = fsolve(rhs, y_guess)
print(f"Fixed point: ({fp[0]:.4f}, {fp[1]:.4f})")

# Classify via eigenvalues of Jacobian
J = jacobian(fp)
eigenvalues = np.linalg.eigvals(J)
print(f"Eigenvalues: {eigenvalues}")

# Classification logic
def classify_fixed_point(eigenvalues):
    real_parts = eigenvalues.real
    imag_parts = eigenvalues.imag

    if np.all(np.abs(imag_parts) < 1e-10):
        # Real eigenvalues
        if np.all(real_parts < 0):
            return "Stable node"
        elif np.all(real_parts > 0):
            return "Unstable node"
        else:
            return "Saddle point"
    else:
        # Complex eigenvalues
        if np.all(real_parts < 0):
            return "Stable spiral"
        elif np.all(real_parts > 0):
            return "Unstable spiral"
        else:
            return "Center"

print(f"Classification: {classify_fixed_point(eigenvalues)}")

3. Bifurcation Diagram

def logistic_map_bifurcation(r_range, n_discard=500, n_plot=200):
    """Bifurcation diagram for the logistic map x_{n+1} = r * x_n * (1 - x_n)"""
    r_values = []
    x_values = []

    for r in r_range:
        x = 0.5  # initial condition
        # Discard transients
        for _ in range(n_discard):
            x = r * x * (1 - x)
        # Collect steady-state iterates
        for _ in range(n_plot):
            x = r * x * (1 - x)
            r_values.append(r)
            x_values.append(x)

    return np.array(r_values), np.array(x_values)

r_range = np.linspace(2.5, 4.0, 2000)
r_vals, x_vals = logistic_map_bifurcation(r_range)

fig, ax = plt.subplots(figsize=(12, 7))
ax.scatter(r_vals, x_vals, s=0.01, c='black', alpha=0.5)
ax.set_xlabel('r', fontsize=13)
ax.set_ylabel('x*', fontsize=13)
ax.set_title('Logistic Map Bifurcation Diagram', fontsize=14)
ax.grid(True, alpha=0.2)
plt.savefig('bifurcation.png', dpi=200, bbox_inches='tight')

For continuous systems, use numerical continuation:

def continuous_bifurcation(rhs_func, param_range, y_init, param_name='mu'):
    """
    Simple parameter continuation for fixed point bifurcations.
    Tracks fixed points as a parameter varies.
    """
    fixed_points = []
    stabilities = []
    y_guess = y_init

    for param in param_range:
        fp = fsolve(lambda y: rhs_func(y, param), y_guess, full_output=True)
        if fp[2] == 1:  # converged
            y_guess = fp[0]  # use as next guess
            # Compute stability
            J = np.zeros((len(y_guess), len(y_guess)))
            eps = 1e-8
            f0 = np.array(rhs_func(y_guess, param))
            for j in range(len(y_guess)):
                y_pert = y_guess.copy()
                y_pert[j] += eps
                J[:, j] = (np.array(rhs_func(y_pert, param)) - f0) / eps
            eigs = np.linalg.eigvals(J)
            stable = np.all(eigs.real < 0)
            fixed_points.append(y_guess.copy())
            stabilities.append(stable)
        else:
            fixed_points.append(np.full_like(y_guess, np.nan))
            stabilities.append(False)

    return np.array(fixed_points), np.array(stabilities)

4. Lyapunov Exponent (Variational Method)

def largest_lyapunov(rhs, y0, t_total=100, dt_renorm=1.0, rtol=1e-10):
    """
    Compute largest Lyapunov exponent via variational equations.
    Method: Benettin et al. (1980), Wolf et al. (1985).
    """
    n = len(y0)

    def rhs_with_variational(t, Y):
        y = Y[:n]
        delta = Y[n:2*n]
        dydt = np.array(rhs(t, y))
        # Numerical Jacobian
        eps = 1e-8
        J = np.zeros((n, n))
        f0 = dydt
        for j in range(n):
            y_pert = np.array(y)
            y_pert[j] += eps
            J[:, j] = (np.array(rhs(t, y_pert)) - f0) / eps
        ddelta_dt = J @ delta
        return np.concatenate([dydt, ddelta_dt])

    # Random initial tangent vector
    np.random.seed(42)
    delta0 = np.random.randn(n)
    delta0 /= np.linalg.norm(delta0)
    Y = np.concatenate([y0, delta0])

    lyap_sum = 0.0
    n_steps = int(t_total / dt_renorm)
    running = []

    for i in range(n_steps):
        t0, t1 = i * dt_renorm, (i + 1) * dt_renorm
        sol = solve_ivp(rhs_with_variational, (t0, t1), Y,
                        method='RK45', rtol=rtol, atol=rtol*1e-2)
        Y = sol.y[:, -1]
        norm = np.linalg.norm(Y[n:2*n])
        lyap_sum += np.log(norm)
        Y[n:2*n] /= norm
        running.append(lyap_sum / ((i+1) * dt_renorm))

    lambda_max = lyap_sum / t_total
    return lambda_max, running

# Example: Lorenz system
def lorenz(t, y, sigma=10, rho=28, beta=8/3):
    x, y_, z = y
    return [sigma*(y_-x), x*(rho-z)-y_, x*y_-beta*z]

lam, running = largest_lyapunov(lorenz, [1, 1, 1], t_total=200)
print(f"Largest Lyapunov exponent: λ_max = {lam:.4f} s⁻¹")
print(f"Expected for Lorenz: ~0.91 s⁻¹")

5. Poincare Section

def poincare_section(rhs, y0, t_total, section_var=2, section_val=None,
                     direction='positive'):
    """
    Compute Poincare section by detecting crossings of a hyperplane.
    section_var: index of variable defining the section
    section_val: value at which to take the section (default: mean of trajectory)
    """
    sol = solve_ivp(rhs, (0, t_total), y0,
                    t_eval=np.linspace(0, t_total, int(t_total * 1000)),
                    rtol=1e-10, atol=1e-12)

    if section_val is None:
        section_val = np.mean(sol.y[section_var])

    # Find crossings
    y_section = sol.y[section_var] - section_val
    crossings = []

    for i in range(1, len(y_section)):
        if direction == 'positive' and y_section[i-1] < 0 and y_section[i] >= 0:
            # Linear interpolation for precise crossing
            frac = -y_section[i-1] / (y_section[i] - y_section[i-1])
            point = sol.y[:, i-1] + frac * (sol.y[:, i] - sol.y[:, i-1])
            crossings.append(point)
        elif direction == 'negative' and y_section[i-1] >= 0 and y_section[i] < 0:
            frac = y_section[i-1] / (y_section[i-1] - y_section[i])
            point = sol.y[:, i-1] + frac * (sol.y[:, i] - sol.y[:, i-1])
            crossings.append(point)

    return np.array(crossings)

# Example: Lorenz attractor Poincare section at z = 27
crossings = poincare_section(lorenz, [1, 1, 1], t_total=500,
                             section_var=2, section_val=27.0)
plt.scatter(crossings[:, 0], crossings[:, 1], s=0.5, c='black')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Poincaré Section of Lorenz Attractor (z = 27)')
plt.savefig('poincare_section.png', dpi=150)

Fixed Point Classification Reference

Common Bifurcation Types

Troubleshooting