Dynamical Systems Analysis

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')

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)}")

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')

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)

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⁻¹")

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)