ch199 — Matrix Exponentials and Dynamical Systems

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-whitegrid')
rng = np.random.default_rng(42)

def matrix_exp(A, t):
    """Compute e^(At) via eigendecomposition. Valid for diagonalizable A."""
    eigvals, V = np.linalg.eig(A)
    exp_lambda = np.exp(eigvals * t)
    return (V * exp_lambda[None, :]) @ np.linalg.inv(V)

def simulate_ode(A, x0, t_end=8.0, n_steps=300):
    t_vals = np.linspace(0, t_end, n_steps)
    traj = np.array([matrix_exp(A, t).real @ x0 for t in t_vals])
    return t_vals, traj

# Three 2x2 systems: stable, neutral, unstable
A_stable   = np.array([[-1., 2.],[-2.,-1.]])
A_neutral  = np.array([[ 0., 1.],[-1., 0.]])
A_unstable = np.array([[ 0.5, 1.],[-1., 0.5]])

x0 = np.array([1., 0.])
fig, axes = plt.subplots(2, 3, figsize=(13, 7))

for col, (A, title) in enumerate([
    (A_stable,   'Stable'),
    (A_neutral,  'Neutral (oscillator)'),
    (A_unstable, 'Unstable'),
]):
    eigs = np.linalg.eigvals(A)
    t_vals, traj = simulate_ode(A, x0)

    axes[0, col].plot(t_vals, traj[:,0], color='steelblue', label='x₁')
    axes[0, col].plot(t_vals, traj[:,1], color='darkorange', label='x₂')
    axes[0, col].set_title(f'{title}\nλ={eigs[0]:.2f}', fontsize=9)
    axes[0, col].set_xlabel('t'); axes[0, col].set_ylabel('State')
    axes[0, col].legend(fontsize=7)

    axes[1, col].plot(traj[:,0], traj[:,1], color='steelblue')
    axes[1, col].plot(*x0, 'go', ms=8, label='Start')
    axes[1, col].plot(*traj[-1], 'rs', ms=8, label='End')
    axes[1, col].set_xlabel('x₁'); axes[1, col].set_ylabel('x₂')
    axes[1, col].set_title('Phase portrait', fontsize=9)
    axes[1, col].legend(fontsize=7)
    axes[1, col].set_aspect('equal')

plt.suptitle("Linear ODE x' = Ax: Stable, Neutral, Unstable", fontsize=10)
plt.tight_layout()
plt.show()

def matrix_exp_scaling_squaring(A, n_terms=15):
    """
    Matrix exponential via scaling-and-squaring + Taylor series.
    More numerically stable than naive Taylor for large ||A||.
    
    Idea: e^A = (e^(A/2^s))^(2^s), choosing s so ||A/2^s|| is small.
    """
    norm_A = np.linalg.norm(A, ord=1)
    s = max(0, int(np.ceil(np.log2(max(norm_A, 1e-10)))))
    A_sc = A / (2**s)

    # Taylor series for e^(A_sc)
    result = np.eye(A.shape[0])
    term   = np.eye(A.shape[0])
    for k in range(1, n_terms):
        term   = term @ A_sc / k
        result = result + term
        if np.linalg.norm(term) < 1e-15 * np.linalg.norm(result):
            break

    # Square back s times
    for _ in range(s):
        result = result @ result
    return result

# Validate against eigendecomposition
A_test = rng.normal(0, 0.5, (4, 4))
A_sym  = (A_test - A_test.T) / 2   # antisymmetric: eigenvalues purely imaginary

t_test = 2.5
eA_eig = matrix_exp(A_sym, t_test)
eA_ss  = matrix_exp_scaling_squaring(A_sym * t_test)
print(f'Scaling-squaring vs eig error: {np.linalg.norm(eA_eig - eA_ss):.2e}')

# Verify: e^A for antisymmetric A should be orthogonal (det=1, A.T A = I)
print(f'||R.T R - I|| (should be ~0): {np.linalg.norm(eA_ss.T @ eA_ss - np.eye(4)):.2e}')
print(f'det(e^A) (should be 1): {np.linalg.det(eA_ss):.6f}')

# --- Experiment: Stability boundary in parameter space ---
# For the system A = [[a, -1],[1, a]], the stability boundary is Re(λ) = 0 → a = 0.
# Try changing: A_VALS, T_SIM

A_VALS = [-0.5, -0.1, 0.0, 0.1, 0.5]   # <-- stability parameter
T_SIM  = 10.0

fig, axes = plt.subplots(1, len(A_VALS), figsize=(14, 3), sharey=False)
for ax, a_val in zip(axes, A_VALS):
    A_param = np.array([[a_val, -1.],[1., a_val]])
    eigs    = np.linalg.eigvals(A_param)
    t_vals, traj = simulate_ode(A_param, np.array([1., 0.]), T_SIM)
    ax.plot(t_vals, np.linalg.norm(traj, axis=1))
    color = 'steelblue' if a_val < 0 else ('gray' if a_val == 0 else 'crimson')
    ax.set_title(f'a={a_val}\nλ={eigs[0]:.2f}', fontsize=8, color=color)
    ax.set_xlabel('t'); ax.set_ylabel('||x(t)||')

plt.suptitle('Stability as a Function of Parameter a', fontsize=10)
plt.tight_layout()
plt.show()
print('a<0: stable (decays), a=0: marginally stable, a>0: unstable (grows)')

# --- Mini Project: Coupled Oscillator Network ---
# n coupled pendulums on a ring: x'' = -L x where L is the graph Laplacian.
# Reformulate as first-order system and simulate.

n_pend = 6   # number of pendulums

# 1D ring Laplacian
L_ring = 2*np.eye(n_pend) - np.roll(np.eye(n_pend), 1, axis=1) - np.roll(np.eye(n_pend), -1, axis=1)

# First-order system: state = [x, v]
# [x']   [ 0   I ] [x]
# [v'] = [-L   0 ] [v]
Z  = np.zeros((n_pend, n_pend))
I  = np.eye(n_pend)
A_osc = np.block([[Z, I], [-L_ring, Z]])

eigs = np.linalg.eigvals(A_osc)
print('Eigenvalues (should be purely imaginary for conservative system):')
print(np.round(eigs, 4))
print('Max real part:', np.real(eigs).max())

# Initial condition: perturb first pendulum
x0 = np.zeros(2 * n_pend)
x0[0] = 1.0   # first pendulum displaced

t_vals = np.linspace(0, 20, 400)
traj = np.array([matrix_exp(A_osc, t).real @ x0 for t in t_vals])

fig, axes = plt.subplots(2, 1, figsize=(12, 6))
for i in range(n_pend):
    axes[0].plot(t_vals, traj[:, i], label=f'Pendulum {i}')
axes[0].set_xlabel('t')
axes[0].set_ylabel('Displacement')
axes[0].set_title('Coupled Pendulum Oscillators (ring topology)')
axes[0].legend(fontsize=7, ncol=3)

im = axes[1].imshow(traj[:, :n_pend].T, aspect='auto', origin='lower',
                   extent=[0, 20, 0, n_pend-1], cmap='RdBu_r')
axes[1].set_xlabel('t')
axes[1].set_ylabel('Pendulum index')
axes[1].set_title('Displacement (spacetime diagram)')
plt.colorbar(im, ax=axes[1])
plt.tight_layout()
plt.show()
print('Normal modes visible as standing wave patterns in the spacetime diagram.')