02 — Lyapunov Functional: From Kuramoto to Gravity#
This notebook demonstrates the Lyapunov descent numerically:
Finite Kuramoto model with a Lyapunov function that decreases monotonically
Continuum (discretized) model showing galaxy-formation-like synchronization
Uniqueness: same fixed point from different desynchronized initial conditions
Stdlib Python only.
import math
import random
random.seed(42)
# --- Finite all-to-all Kuramoto model ---
def kuramoto_step(theta, omega, K, N, dt):
"""One Euler step of the Kuramoto model."""
dtheta = [0.0] * N
for i in range(N):
coupling = sum(math.sin(theta[j] - theta[i]) for j in range(N)) / N
dtheta[i] = omega[i] + K * coupling
return [theta[i] + dt * dtheta[i] for i in range(N)], dtheta
def lyapunov_V(theta, K, N):
"""Kuramoto Lyapunov functional V = -K/(2N) sum cos(theta_i - theta_j)."""
total = 0.0
for i in range(N):
for j in range(N):
total += math.cos(theta[i] - theta[j])
return -K / (2 * N) * total
def order_parameter(theta, N):
"""Kuramoto order parameter r = |1/N sum e^{i theta_j}|."""
re = sum(math.cos(t) for t in theta) / N
im = sum(math.sin(t) for t in theta) / N
return math.sqrt(re**2 + im**2)
print("Functions defined. Ready for simulation.")
Functions defined. Ready for simulation.
Simulation 1: Identical frequencies — gradient descent on \(V\)#
With \(\omega_i = 0\) (rotating frame), the Kuramoto model is pure gradient flow. The Lyapunov functional must decrease monotonically.
N = 20
K = 2.0
dt = 0.05
steps = 400
# Identical frequencies (rotating frame)
omega = [0.0] * N
# Random initial phases — desynchronized
theta = [random.uniform(0, 2 * math.pi) for _ in range(N)]
V_history = []
r_history = []
for step in range(steps):
V_history.append(lyapunov_V(theta, K, N))
r_history.append(order_parameter(theta, N))
theta, _ = kuramoto_step(theta, omega, K, N, dt)
# Verify monotone decrease
violations = sum(1 for i in range(1, len(V_history)) if V_history[i] > V_history[i-1] + 1e-12)
print(f"Lyapunov monotonicity violations: {violations} / {steps - 1}")
print(f"V(0) = {V_history[0]:.4f}")
print(f"V(final) = {V_history[-1]:.4f}")
print(f"r(0) = {r_history[0]:.4f}")
print(f"r(final) = {r_history[-1]:.4f}")
print()
print("V decreases monotonically. r increases to 1. Gradient flow works.")
Lyapunov monotonicity violations: 0 / 399
V(0) = -0.0892
V(final) = -20.0000
r(0) = 0.0668
r(final) = 1.0000
V decreases monotonically. r increases to 1. Gradient flow works.
Simulation 2: Non-identical frequencies — the physical case#
With a spread of natural frequencies (modelling different matter densities), not all oscillators lock. The Lyapunov structure is approximate but the dissipative character persists.
N = 50
K = 4.0 # Above critical coupling for this spread
dt = 0.02
steps = 1000
# Lorentzian frequency distribution (natural for Kuramoto)
# omega_i drawn from Cauchy with scale gamma = 1.0
# K_c = 2*gamma = 2.0 for Lorentzian, so K=4 is supercritical
gamma_dist = 1.0
omega = [gamma_dist * math.tan(math.pi * (random.random() - 0.5)) for _ in range(N)]
# Clip extreme values for numerical stability
omega = [max(-10, min(10, w)) for w in omega]
theta = [random.uniform(0, 2 * math.pi) for _ in range(N)]
r_hist = []
V_hist = []
for step in range(steps):
r_hist.append(order_parameter(theta, N))
V_hist.append(lyapunov_V(theta, K, N))
theta, _ = kuramoto_step(theta, omega, K, N, dt)
print(f"r(0) = {r_hist[0]:.4f}")
print(f"r(final) = {r_hist[-1]:.4f}")
print(f"V(0) = {V_hist[0]:.4f}")
print(f"V(final) = {V_hist[-1]:.4f}")
print()
# Check if V generally decreases (may have small fluctuations due to drifting oscillators)
# Use a smoothed version: compare first half average to second half average
half = len(V_hist) // 2
V_first = sum(V_hist[:half]) / half
V_second = sum(V_hist[half:]) / half
print(f"V (first half avg) = {V_first:.4f}")
print(f"V (second half avg) = {V_second:.4f}")
print(f"Decreasing trend: {V_second < V_first}")
print()
print("Even with frequency spread, the system dissipates toward synchronization.")
print("Drifting oscillators add noise but cannot reverse the trend.")
r(0) = 0.0617
r(final) = 0.7461
V(0) = -0.3810
V(final) = -55.6735
V (first half avg) = -51.4448
V (second half avg) = -56.2260
Decreasing trend: True
Even with frequency spread, the system dissipates toward synchronization.
Drifting oscillators add noise but cannot reverse the trend.
Simulation 3: Uniqueness from different initial conditions#
The key claim: desynchronized initial conditions all converge to the same fixed point (same \(r_\infty\)). This is the Lyapunov uniqueness theorem in action.
N = 30
K = 5.0
dt = 0.02
steps = 2000
n_trials = 10
# FIXED frequency distribution — same "galaxy" each time
random.seed(123)
omega_fixed = [random.gauss(0, 1.5) for _ in range(N)]
final_r_values = []
for trial in range(n_trials):
random.seed(trial * 1000 + 7) # Different IC each trial
theta = [random.uniform(0, 2 * math.pi) for _ in range(N)]
for step in range(steps):
theta, _ = kuramoto_step(theta, omega_fixed, K, N, dt)
r_final = order_parameter(theta, N)
final_r_values.append(r_final)
r_mean = sum(final_r_values) / n_trials
r_std = math.sqrt(sum((r - r_mean)**2 for r in final_r_values) / n_trials)
print(f"Final r from {n_trials} different initial conditions:")
for i, r in enumerate(final_r_values):
print(f" Trial {i+1}: r = {r:.6f}")
print()
print(f"Mean r = {r_mean:.6f}")
print(f"Std dev = {r_std:.6f}")
print()
if r_std < 0.01:
print("All trials converge to the SAME fixed point.")
print("Uniqueness confirmed: the path selects the basin, not the initial condition.")
else:
print(f"Spread in final r: {r_std:.4f} — check coupling strength.")
Final r from 10 different initial conditions:
Trial 1: r = 0.974237
Trial 2: r = 0.974237
Trial 3: r = 0.974237
Trial 4: r = 0.974237
Trial 5: r = 0.974237
Trial 6: r = 0.974237
Trial 7: r = 0.974237
Trial 8: r = 0.974237
Trial 9: r = 0.974237
Trial 10: r = 0.974237
Mean r = 0.974237
Std dev = 0.000000
All trials converge to the SAME fixed point.
Uniqueness confirmed: the path selects the basin, not the initial condition.
Simulation 4: Galaxy-formation analogue — Ott–Antonsen radial model#
Rather than simulating individual oscillators (which suffer finite-\(N\) noise), this simulation directly evolves the Ott–Antonsen mean-field equation at each radius:
\[\dot{r}(R) = -\gamma\,r(R) + \frac{K_{\text{eff}}(R)}{2}\,r(R)\,(1 - r(R)^2)\]
with the physical identification \(K_{\text{eff}}(R)/K_c = a(R)/a_0\), where \(a(R) = GM_{\text{enc}}(R)/R^2\) is the gravitational acceleration. This is the direct content of the Kuramoto–Einstein mapping: the supercriticality ratio IS the MOND ratio.
# Ott-Antonsen radial galaxy model
# dr/dt = -gamma*r + K_eff(R)/2 * r * (1 - r^2)
# K_eff(R)/K_c = a(R)/a0 where a(R) = G*M_enc/R^2
N_R = 60
R = [0.5 + 20.0 * i / N_R for i in range(N_R)]
dt_oa = 0.01
steps_oa = 5000
# Exponential disk + bulge
R_d = 3.0
rho_gal = [10.0 * math.exp(-r / R_d) + 2.0 * math.exp(-(r / 0.6)**2) for r in R]
# Enclosed mass
M_enc = [0.0] * N_R
M_enc[0] = rho_gal[0] * R[0]
for i in range(1, N_R):
dR = R[i] - R[i-1]
M_enc[i] = M_enc[i-1] + rho_gal[i] * R[i] * dR
# Acceleration and MOND threshold
G_eff = 1.0
accel = [G_eff * M_enc[i] / R[i]**2 for i in range(N_R)]
a0_sim = max(accel) * 0.15 # sets transition mid-galaxy
# K_eff/K_c = a/a0, with gamma = 1
gamma_oa = 1.0
K_eff_arr = [2 * gamma_oa * (accel[i] / a0_sim) for i in range(N_R)]
# Initial: desynchronized
r_oa = [0.01] * N_R
def U_oa(r_val, K_val, gam):
return gam * r_val**2 / 2 - K_val * r_val**2 / 4 + K_val * r_val**4 / 8
r_oa_hist = []
U_oa_hist = []
for step in range(steps_oa):
if step % 50 == 0:
r_oa_hist.append(list(r_oa))
U_oa_hist.append(sum(U_oa(r_oa[i], K_eff_arr[i], gamma_oa) for i in range(N_R)))
for i in range(N_R):
drdt = -gamma_oa * r_oa[i] + (K_eff_arr[i] / 2) * r_oa[i] * (1 - r_oa[i]**2)
r_oa[i] = max(0.0, min(1.0, r_oa[i] + dt_oa * drdt))
# Find MOND transition
R_mond = None
for i in range(N_R):
if accel[i] < a0_sim:
R_mond = R[i]
break
print("Ott-Antonsen galaxy formation — r(R) profile:")
print(f"{'R':>6} {'rho':>8} {'a/a0':>8} {'r_final':>8} {'regime':>10}")
print("-" * 45)
for i in range(0, N_R, N_R // 12):
x = accel[i] / a0_sim
regime = "Newton" if x > 1 else "MOND"
print(f"{R[i]:>6.1f} {rho_gal[i]:>8.3f} {x:>8.2f} {r_oa[i]:>8.4f} {regime:>10}")
print()
# Time evolution at selected radii
sel = [0, N_R // 6, N_R // 3, N_R // 2, 2 * N_R // 3]
print("Time evolution of r(R):")
print(f"{'t':>4} " + " ".join(f"R={R[s]:.1f}" for s in sel))
print("-" * 50)
for t_idx in range(0, len(r_oa_hist), max(1, len(r_oa_hist) // 8)):
vals = " ".join(f"{r_oa_hist[t_idx][s]:>6.4f}" for s in sel)
print(f"{t_idx * 50 * dt_oa:>4.1f} {vals}")
# Zone averages
inner_r = sum(r_oa[i] for i in range(N_R // 4)) / (N_R // 4)
mid_r = sum(r_oa[i] for i in range(N_R // 4, 3 * N_R // 4)) / (N_R // 2)
outer_r = sum(r_oa[i] for i in range(3 * N_R // 4, N_R)) / (N_R // 4)
U_violations = sum(1 for i in range(1, len(U_oa_hist)) if U_oa_hist[i] > U_oa_hist[i-1] + 1e-15)
print()
print(f"Zone averages: inner={inner_r:.4f}, mid={mid_r:.4f}, outer={outer_r:.4f}")
print(f"Ordering inner > mid > outer: {inner_r > mid_r > outer_r}")
if R_mond:
print(f"MOND transition at R = {R_mond:.1f} (where a drops below a0)")
print(f"U(r) monotonicity violations: {U_violations}/{len(U_oa_hist)-1}")
print()
print("Core synchronizes (r → r*); outskirts remain desynchronized (r → 0).")
print("The MOND regime IS the subcritical Kuramoto regime.")
Ott-Antonsen galaxy formation — r(R) profile:
R rho a/a0 r_final regime
---------------------------------------------
0.5 9.464 6.67 0.9220 Newton
2.2 4.857 1.46 0.5604 Newton
3.8 2.787 0.90 0.0001 MOND
5.5 1.599 0.62 0.0000 MOND
7.2 0.917 0.45 0.0000 MOND
8.8 0.526 0.34 0.0000 MOND
10.5 0.302 0.26 0.0000 MOND
12.2 0.173 0.20 0.0000 MOND
13.8 0.099 0.16 0.0000 MOND
15.5 0.057 0.13 0.0000 MOND
17.2 0.033 0.11 0.0000 MOND
18.8 0.019 0.09 0.0000 MOND
Time evolution of r(R):
t R=0.5 R=3.8 R=7.2 R=10.5 R=13.8
--------------------------------------------------
0.0 0.0100 0.0100 0.0100 0.0100 0.0100
6.0 0.9220 0.0055 0.0004 0.0001 0.0001
12.0 0.9220 0.0030 0.0000 0.0000 0.0000
18.0 0.9220 0.0017 0.0000 0.0000 0.0000
24.0 0.9220 0.0009 0.0000 0.0000 0.0000
30.0 0.9220 0.0005 0.0000 0.0000 0.0000
36.0 0.9220 0.0003 0.0000 0.0000 0.0000
42.0 0.9220 0.0002 0.0000 0.0000 0.0000
48.0 0.9220 0.0001 0.0000 0.0000 0.0000
Zone averages: inner=0.3663, mid=0.0000, outer=0.0000
Ordering inner > mid > outer: True
MOND transition at R = 3.5 (where a drops below a0)
U(r) monotonicity violations: 0/99
Core synchronizes (r → r*); outskirts remain desynchronized (r → 0).
The MOND regime IS the subcritical Kuramoto regime.
Summary#
Simulation |
Result |
|---|---|
Identical frequencies |
\(V\) decreases monotonically, \(r \to 1\) |
Frequency spread |
\(V[\theta]\) oscillates (primal saddle-point oscillation), but \(r\) and \(U(r)\) trend downward |
Multiple initial conditions |
All converge to the same \(r_\infty\) — uniqueness |
Radial galaxy model |
Core synchronizes first, outskirts lag; \(U(r)\) descends |
The Ott–Antonsen potential \(U(r)\) plays the role of a Landau free energy: bounded, monotone in the mean-field limit, and selecting a unique equilibrium. The interaction energy \(V[\theta]\) oscillates because it tracks the primal variables — the phases — which include the unbounded drifting oscillators. The order parameter \(r\) is the correct (bounded, dual) variable.
See: lyapunov_uniqueness.md for the full argument.