Chaos

READme and animation!

Phase Map: Finding Islands of Stability in Chaos

import matplotlib.pyplot as plt
from pendulum_equations import _lambdifygenerated as accel# 
# Number of pixels per axis (e.g., 200x200 = 40,000 pendulums simulated at once)
GRID_RES = 200  
L1 = L2 = L3 = 1.0
m1 = m2 = m3 = 1.0g = 9.81
theta1_range = np.linspace(-np.pi, np.pi, GRID_RES)
theta2_range = np.linspace(-np.pi, np.pi, GRID_RES)
T1, T2 = np.meshgrid(theta1_range, theta2_range)
# Initialize a massive state grid: Shape is (GRID_RES, GRID_RES, 6)
state_grid = np.zeros((GRID_RES, GRID_RES, 6), dtype=float)
state_grid[:, :, 0] = T1           
# X-axis of pixels maps to unique theta1 valuesstate_grid[:, :, 1] = T2           
# Y-axis of pixels maps to unique theta2 valuesstate_grid[:, :, 2] = np.pi / 4    
# Keep theta3 constant at 45 degrees for allstate_grid[:, :, 3] = 0.0          
# omega1 = 0state_grid[:, :, 4] = 0.0          
# omega2 = 0state_grid[:, :, 5] = 0.0          
# omega3 = 0# 

def vectorized_derivatives(state):    
# Unpack the 3D grid axes    
t1 = state[:, :, 0]    
t2 = state[:, :, 1]    
t3 = state[:, :, 2]    
w1 = state[:, :, 3]    
w2 = state[:, :, 4]    
w3 = state[:, :, 5]    
a1, a2, a3 = accel(t1, t2, t3, w1, w2, w3, L1, L2, L3, m1, m2, m3, g)    
#Runge-Kutta
return np.stack([w1, w2, w3, a1, a2, a3], axis=-1)
def rk4_step(state, dt):    
k1 = vectorized_derivatives(state)    
k2 = vectorized_derivatives(state + 0.5 * dt * k1)    
k3 = vectorized_derivatives(state + 0.5 * dt * k2)    
k4 = vectorized_derivatives(state + dt * k3)    
return state + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
dt = 0.0005
total_steps = 15000  
# Total time window to monitor flipping / divergence
# To measure chaos vs stability, track the max velocity reached by each pendulum.
# Chaotic pendulums will violently whip around (high max velocity/flipping).
# Stable islands will stay tightly bound (very low max velocity).
max_omega_observed = np.zeros((GRID_RES, GRID_RES))
print(f"Simulating {GRID_RES}x{GRID_RES} ({GRID_RES**2}) pendulums simultaneously snoopsawggggg...")
for step in range(total_steps):    state_grid = rk4_step(state_grid, dt)        
# Extract instantaneous absolute speeds of the middle and bottom bobs    
current_speeds = np.abs(state_grid[:, :, 4]) + np.abs(state_grid[:, :, 5])        
# Store the peak (peak, get it, peak) velocity hit by each specific coordinate pixel    
max_omega_observed = np.maximum(max_omega_observed, current_speeds)        if step % 3000 == 0 and step > 0:        
print(f"  -> Progress: {int((step/total_steps)*100)}% complete...")
print("Simulation finished! Puzzah! Rendering Stability Map...")# 
plt.figure(figsize=(8, 8))# Logarithmic scaling brings out hidden geometric structures and fractal edges
stability_map = np.log1p(max_omega_observed)
plt.imshow(    stability_map,     
extent=[-np.pi, np.pi, -np.pi, np.pi],     origin='lower',     cmap='inferno')
plt.colorbar(label='Chaos Metric (Log Peak Angular Velocity)')
plt.xlabel(r'Initial $\theta_1$ (rad)')
plt.ylabel(r'Initial $\theta_2$ (rad)')plt.title('Triple Pendulum Phase Map: Finding Islands of Stability')
plt.show()

Simulating the triple pendulum

import sympy as sp import inspectt=sp.symbols('t')theta1=sp.Function('theta1')(t)theta2=sp.Function('theta2')(t)theta3=sp.Function('theta3')(t)L1, L2, L3 = sp.symbols('L1 L2 L3')m1, m2, m3 = sp.symbols('m1 m2 m3')x1= L1*sp.sin(theta1)y1= -L1*sp.cos(theta1)x2= x1 + L2*sp.sin(theta2)y2= y1 - L2*sp.cos(theta2)x3= x2 + L3*sp.sin(theta3)y3= y2 - L3*sp.cos(theta3)T1,T2,T3 = sp.symbols('T1 T2 T3')V1,V2,V3 = sp.symbols('V1 V2 V3')x1dd= sp.diff(x1,t,2)x2dd= sp.diff(x2,t,2)x3dd= sp.diff(x3,t,2)y1dd= sp.diff(y1,t,2)y2dd= sp.diff(y2,t,2)y3dd= sp.diff(y3,t,2)vx1=sp.diff(x1,t)vy1=sp.diff(y1,t)vx2=sp.diff(x2,t)vy2=sp.diff(y2,t)vx3=sp.diff(x3,t)vy3=sp.diff(y3,t)T = (sp.Rational(1,2)*m1*(vx1**2 + vy1**2) + sp.Rational(1,2)*m2*(vx2**2 + vy2**2) + sp.Rational(1,2)*m3*(vx3**2 + vy3**2))g = sp.symbols('g')V = m1*g*y1 + m2*g*y2 + m3*g*y3 Lagrangian = T - Vprint("Getting Lagrangian!!!")Lagrangian = sp.simplify(Lagrangian)print(Lagrangian)print("Hey, we got the Lagrangian!")#Euler-Lagrange equations:Euler1= sp.diff(sp.diff(Lagrangian, sp.diff(theta1,t)), t) - sp.diff(Lagrangian, theta1)Euler2= sp.diff(sp.diff(Lagrangian, sp.diff(theta2,t)), t) - sp.diff(Lagrangian, theta2)Euler3= sp.diff(sp.diff(Lagrangian, sp.diff(theta3,t)), t) - sp.diff(Lagrangian, theta3)theta1dd=sp.diff(theta1,t,2)theta2dd=sp.diff(theta2,t,2)theta3dd=sp.diff(theta3,t,2)derivs=[theta1dd, theta2dd, theta3dd]solution=sp.solve((Euler1, Euler2, Euler3), derivs, dict=True)alpha1 = solution[0][derivs[0]]alpha2 = solution[0][derivs[1]]alpha3 = solution[0][derivs[2]]#now i want to simplify the solution and then convert it to a function that can be used in numerical simulations.print("alpha1=", alpha1)print("******************************************************************************************")print("alpha2=", alpha2)print("******************************************************************************************")print("alpha3=", alpha3)accel_func = sp.lambdify( (theta1, theta2, theta3, sp.diff(theta1,t), sp.diff(theta2,t), sp.diff(theta3,t), L1, L2, L3, m1, m2, m3, g), (solution[0][theta1dd], solution[0][theta2dd], solution[0][theta3dd]), "numpy")print("Writing equations with proper NumPy math imports...")with open("pendulum_equations.py", "w") as f: f.write("import numpy as np\n") # This injects the missing sin and cos tools directly into the math module! f.write("from numpy import sin, cos\n\n") f.write(inspect.getsource(accel_func))print("Successfully saved math script to 'pendulum_equations.py'!")# ... Your existing derivation code that creates 'accel' ...# #x1dd= -L1*((theta1.diff()*theta1.diff())*sp.sin(theta1)) + (theta1.diff().diff()*(L1*sp.cos(theta1)))# x2dd=x1dd - (L2*theta2.diff()*theta2.diff()*sp.sin(theta2)) + (theta2.diff().diff()*L2*sp.cos(theta2))# x3dd=x2dd - (theta3.diff()*theta3.diff()*L3*sp.sin(theta3)) + (theta3.diff().diff()*sp.sin(theta3)*L3)# y1dd= (theta1.diff()*theta1.diff()*L1*sp.cos(theta1))+(theta1.diff().diff()*L1*sp.sin(theta1))# y2dd=y1dd + (theta2.diff()*theta2.diff()*L2*sp.cos(theta2))+(theta2.diff().diff()*L2*sp.sin(theta2))# y3dd=y2dd + (theta3.diff()*theta3.diff()*L3*sp.cos(theta3))+(theta3.diff().diff()*L3*sp.sin(theta3))# vx1=sp.diff(x1,t)# vy1=sp.diff(y1,t)# vx2=sp.diff(x2,t)# vy2=sp.diff(y2,t)# vx3=sp.diff(x3,t)# vy3=sp.diff(y3,t)# T= # #Equation1 # eq1 =sp.Eq(m1*x1dd, (-T1*sp.sin(theta1))+(T2*sp.sin(theta2)))# #Equation2# eq2 = sp.Eq(m1*y1dd, (T1*sp.cos(theta1))-(T2*sp.cos(theta2))-(m1*9.81))# #Equation3# eq3 = sp.Eq(m2*x2dd, (-T2*sp.sin(theta2))+(T3*sp.sin(theta3)))# #Equation4# eq4 = sp.Eq(m2*y2dd, (T2*sp.cos(theta2))-(T3*sp.cos(theta3))-(m2*9.81))# #Equation5# eq5 = sp.Eq(m3*x3dd, (-T3*sp.sin(theta3))# #Equation6# eq6 = sp.Eq(m3*y3dd, (T3*sp.cos(theta3)) - (m3*9.81))# sp.solve([eq1, eq2, eq3, eq4, eq5, eq6], (theta1.diff().diff(), theta2.diff().diff(), theta3.diff().diff(), T1, T2, T3)), dict=True