1cejhlnber6qmkhrl1 Ctha.gif

Modelado de un sistema caótico

El péndulo es un sistema de física clásica con el que todos estamos familiarizados. Ya sea un reloj de pie o un niño en un columpio, hemos visto el movimiento regular y periódico del péndulo. Un péndulo simple está bien definido en la física clásica, pero el péndulo doble (un péndulo unido al extremo de otro péndulo) es caos literal. En este artículo, nos basaremos en nuestra comprensión intuitiva de los péndulos y modelaremos el caos del péndulo doble. La física es interesante y los métodos numéricos necesarios son una herramienta esencial en el arsenal de cualquiera.

Figura 1: Ejemplo de un péndulo doble caótico

En este artículo vamos a:

  • Aprenda sobre el movimiento armónico y modele el comportamiento de un solo péndulo
  • Aprende los fundamentos de la teoría del caos.
  • Modelar numéricamente el comportamiento caótico de un péndulo doble

Movimiento armónico simple

Describimos el movimiento oscilante periódico de un péndulo como movimiento armónico. El movimiento armónico ocurre cuando hay movimiento en un sistema que es balanceado por una fuerza restauradora proporcional en la dirección opuesta a dicho movimiento. Vemos un ejemplo de esto en la figura 2 donde una masa en un resorte está siendo jalada hacia abajo debido a la gravedad, pero esto pone energía en el resorte que luego retrocede y jala la masa hacia arriba. Al lado del sistema de resortes, vemos la altura de la masa dando vueltas en un círculo llamado diagrama fasorial que ilustra aún más el movimiento regular del sistema.

Figura 2: Ejemplo de movimiento armónico simple de una masa sobre un resorte

El movimiento armónico puede ser amortiguado (disminuyendo en amplitud debido a las fuerzas de arrastre) o impulsado (aumentando en amplitud debido a la adición de una fuerza externa), pero comenzaremos con el caso más simple de movimiento armónico indefinido sin fuerzas externas actuando sobre él (movimiento no amortiguado). ). Este tipo de movimiento es una buena aproximación para modelar un solo péndulo que oscila en un ángulo pequeño/baja amplitud. En este caso, podemos modelar el movimiento con la ecuación 1 a continuación.

Ecuación 1: Movimiento armónico simple para un péndulo de ángulo pequeño

Podemos codificar fácilmente esta función y simular un péndulo simple a lo largo del tiempo.

def simple_pendulum(theta_0, omega, t, phi):
theta = theta_0*np.cos(omega*t + phi)
return theta

#parameters of our system
theta_0 = np.radians(15) #degrees to radians

g = 9.8 #m/s^2
l = 1.0 #m
omega = np.sqrt(g/l)

phi = 0 #for small angle

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals
theta = []
for t in time_span:
theta.append(simple_pendulum(theta_0, omega, t, phi))

#Convert back to cartesian coordinates
x = l*np.sin(theta)
y = -l*np.cos(theta) #negative to make sure the pendulum is facing down

Figura 3: simulación de péndulo simple

Movimiento de péndulo completo con mecánica lagrangiana

Un péndulo de ángulo pequeño simple es un buen comienzo, pero queremos ir más allá y modelar el movimiento de un péndulo completo. Como ya no podemos usar aproximaciones de ángulo pequeño es mejor modelar el péndulo usando mecánica lagrangiana. Esta es una herramienta esencial en física que nos cambia de mirar las fuerzas en un sistema a mirar la energía en un sistema. Estamos cambiando nuestro marco de referencia de fuerza impulsora vs fuerza restauradora a cinético contra potencial energía.

El Lagrangain es la diferencia entre la energía cinética y potencial dada en la ecuación 2.

Ecuación 2: El Lagrangiano

Sustituyendo en la Cinética y el Potencial de un péndulo dado en la ecuación 3 se obtiene el Lagrangain para un péndulo visto en la ecuación 4

Ecuación 3: Energía cinética y potencial de un péndulo
Ecuación 4: Lagrangiano para un péndulo

Con el Lagrangiano para un péndulo ahora describimos la energía de nuestro sistema. Hay un último paso matemático para transformar esto en algo sobre lo que podamos construir una simulación. Necesitamos volver a la referencia dinámica / orientada a la fuerza desde la referencia de energía usando el Ecuación de Euler-Lagrange. Usando esta ecuación podemos usar el Lagrangiano para obtener la aceleración angular de nuestro péndulo.

Ecuación 5: aceleración angular de la ecuación de Euler-Lagrange

Después de hacer las matemáticas, tenemos la aceleración angular que podemos usar para obtener la velocidad angular y el ángulo mismo. Esto requerirá algunos integracion numerica que se presentará en nuestra simulación de péndulo completo. Incluso para un solo péndulo, el dinámica no lineal significa que no hay una solución analítica para resolver para theta, de ahí la necesidad de una solución numérica. La integración es bastante simple (pero poderosa), usamos la aceleración angular para actualizar la velocidad angular y la velocidad angular para actualizar theta sumando la primera cantidad a la última y multiplicándola por algún paso de tiempo. Esto nos da una aproximación del área bajo la curva de aceleración/velocidad. Cuanto menor sea el paso de tiempo, más precisa será la aproximación.

def full_pendulum(g,l,theta,theta_velocity, time_step):
#Numerical Integration
theta_acceleration = -(g/l)*np.sin(theta) #Get acceleration
theta_velocity += time_step*theta_acceleration #Update velocity with acceleration
theta += time_step*theta_velocity #Update angle with angular velocity
return theta, theta_velocity

g = 9.8 #m/s^2
l = 1.0 #m

theta = [np.radians(90)] #theta_0
theta_velocity = 0 #Start with 0 velocity
time_step = 20/300 #Define a time step

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals
for t in time_span:
theta_new, theta_velocity = full_pendulum(g,l,theta[-1], theta_velocity, time_step)
theta.append(theta_new)

#Convert back to cartesian coordinates
x = l*np.sin(theta)
y = -l*np.cos(theta)

Figura 4: Simulación de un péndulo lleno

Hemos simulado un péndulo completo, pero sigue siendo un sistema bien definido. Ha llegado el momento de adentrarse en el caos del doble péndulo.

Caos, en el sentido matemático, se refiere a sistemas que son muy sensibles a sus condiciones iniciales. Incluso pequeños cambios en el inicio del sistema darán lugar a comportamientos muy diferentes a medida que el sistema evolucione. Esto describe perfectamente el movimiento del péndulo doble. A diferencia del péndulo simple, no es un sistema que se comporte bien y evolucionará de una manera muy diferente incluso con cambios leves en el ángulo de inicio.

Para modelar el movimiento del péndulo doble, usaremos el mismo enfoque lagrangiano que antes (ver derivación completa).

También usaremos el mismo esquema de integración numérica que antes al implementar esta ecuación en el código y encontrar theta.

#Get theta1 acceleration 
def theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):
mass1 = -g*(2*m1 + m2)*np.sin(theta1)
mass2 = -m2*g*np.sin(theta1 - 2*theta2)
interaction = -2*np.sin(theta1 - theta2)*m2*np.cos(theta2_velocity**2*l2 + theta1_velocity**2*l1*np.cos(theta1 - theta2))
normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta1_ddot = (mass1 + mass2 + interaction)/normalization

return theta1_ddot

#Get theta2 acceleration
def theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):
system = 2*np.sin(theta1 - theta2)*(theta1_velocity**2*l1*(m1 + m2) + g*(m1 + m2)*np.cos(theta1) + theta2_velocity**2*l2*m2*np.cos(theta1 - theta2))
normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta2_ddot = system/normalization
return theta2_ddot

#Update theta1
def theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):
#Numerical Integration
theta1_velocity += time_step*theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)
theta1 += time_step*theta1_velocity
return theta1, theta1_velocity

#Update theta2
def theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):
#Numerical Integration
theta2_velocity += time_step*theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)
theta2 += time_step*theta2_velocity
return theta2, theta2_velocity

#Run full double pendulum
def double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span):
theta1_list = [theta1]
theta2_list = [theta2]

for t in time_span:
theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)
theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)
theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list) #Pendulum 1 x
y1 = -l1*np.cos(theta1_list) #Pendulum 1 y

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list) #Pendulum 2 x
y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list) #Pendulum 2 y

return x1,y1,x2,y2

#Define system parameters
g = 9.8 #m/s^2

m1 = 1 #kg
m2 = 1 #kg

l1 = 1 #m
l2 = 1 #m

theta1 = np.radians(90)
theta2 = np.radians(45)

theta1_velocity = 0 #m/s
theta2_velocity = 0 #m/s

theta1_list = [theta1]
theta2_list = [theta2]

time_step = 20/300

time_span = np.linspace(0,20,300)
x1,y1,x2,y2 = double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span)

Figura 5: Simulación de doble péndulo

¡Por fin lo hemos hecho! Hemos modelado con éxito un péndulo doble, pero ahora es el momento de observar algo de caos. Nuestra simulación final será de dos péndulos dobles con condiciones de inicio ligeramente diferentes. Fijaremos un péndulo para que tenga un theta 1 de 90 grados y el otro para tener una theta 1 de 91 grados. Veamos qué pasa.

Figura 6: 2 péndulos dobles con condiciones de inicio ligeramente diferentes

Podemos ver que ambos péndulos comienzan con trayectorias similares pero divergen rápidamente. Esto es lo que queremos decir cuando decimos caos, incluso una diferencia de 1 grado en cascadas de ángulo en un comportamiento final muy diferente.

En este artículo aprendimos sobre el movimiento del péndulo y cómo modelarlo. Partimos del modelo de movimiento armónico más simple y construimos hasta el complejo y caótico péndulo doble. En el camino aprendimos sobre el Lagrangiano, el caos y la integración numérica.

El péndulo doble es el ejemplo más simple de un sistema caótico. Estos sistemas existen en todas partes en nuestro mundo desde dinámica poblacional, climae incluso billar. Podemos tomar las lecciones que hemos aprendido del doble péndulo y aplicarlas siempre que nos encontremos con un sistema caótico.

Conclusiones clave

  • Los sistemas caóticos son muy sensibles a las condiciones iniciales y evolucionarán de formas muy diferentes incluso con pequeños cambios en su inicio.
  • Cuando se trata de un sistema, especialmente un sistema caótico, ¿hay otro marco de referencia para mirarlo que haga que sea más fácil trabajar con él? (Como el marco de referencia de fuerza al marco de referencia de energía)
  • Cuando los sistemas se vuelven demasiado complicados, necesitamos implementar soluciones numéricas para resolverlos. Estas soluciones son simples pero poderosas y ofrecen buenas aproximaciones al comportamiento real.

Todas las figuras utilizadas en este artículo fueron creadas por el autor o son de Imágenes Matemáticas y lleno bajo el Licencia de documentación libre GNU 1.2

Mecánica Clásica, John Taylor https://neuroself.files.wordpress.com/2020/09/taylor-2005-classical-mechanics.pdf

Péndulo Simple

def makeGif(x,y,name):
!mkdir frames

counter=0
images = []
for i in range(0,len(x)):
plt.figure(figsize = (6,6))

plt.plot([0,x[i]],[0,y[i]], "o-", color = "b", markersize = 7, linewidth=.7 )
plt.title("Pendulum")
plt.xlim(-1.1,1.1)
plt.ylim(-1.1,1.1)
plt.savefig("frames/" + str(counter)+ ".png")
images.append(imageio.imread("frames/" + str(counter)+ ".png"))
counter += 1
plt.close()

imageio.mimsave(name, images)

!rm -r frames

def simple_pendulum(theta_0, omega, t, phi):
theta = theta_0*np.cos(omega*t + phi)
return theta

#parameters of our system
theta_0 = np.radians(15) #degrees to radians

g = 9.8 #m/s^2
l = 1.0 #m
omega = np.sqrt(g/l)

phi = 0 #for small angle

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals
theta = []
for t in time_span:
theta.append(simple_pendulum(theta_0, omega, t, phi))

x = l*np.sin(theta)
y = -l*np.cos(theta) #negative to make sure the pendulum is facing down

Péndulo

def full_pendulum(g,l,theta,theta_velocity, time_step):
theta_acceleration = -(g/l)*np.sin(theta)
theta_velocity += time_step*theta_acceleration
theta += time_step*theta_velocity
return theta, theta_velocity

g = 9.8 #m/s^2
l = 1.0 #m

theta = [np.radians(90)] #theta_0
theta_velocity = 0
time_step = 20/300

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals
for t in time_span:
theta_new, theta_velocity = full_pendulum(g,l,theta[-1], theta_velocity, time_step)
theta.append(theta_new)

#Convert back to cartesian coordinates
x = l*np.sin(theta)
y = -l*np.cos(theta)

#Use same function from simple pendulum
makeGif(x,y,"pendulum.gif")

Péndulo Doble

def theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):
mass1 = -g*(2*m1 + m2)*np.sin(theta1)
mass2 = -m2*g*np.sin(theta1 - 2*theta2)
interaction = -2*np.sin(theta1 - theta2)*m2*np.cos(theta2_velocity**2*l2 + theta1_velocity**2*l1*np.cos(theta1 - theta2))
normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta1_ddot = (mass1 + mass2 + interaction)/normalization

return theta1_ddot

def theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):
system = 2*np.sin(theta1 - theta2)*(theta1_velocity**2*l1*(m1 + m2) + g*(m1 + m2)*np.cos(theta1) + theta2_velocity**2*l2*m2*np.cos(theta1 - theta2))
normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta2_ddot = system/normalization
return theta2_ddot

def theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

theta1_velocity += time_step*theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)
theta1 += time_step*theta1_velocity
return theta1, theta1_velocity

def theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

theta2_velocity += time_step*theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)
theta2 += time_step*theta2_velocity
return theta2, theta2_velocity

def double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span):
theta1_list = [theta1]
theta2_list = [theta2]

for t in time_span:
theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)
theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)
theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list)
y1 = -l1*np.cos(theta1_list)

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list)
y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list)

return x1,y1,x2,y2

#Define system parameters, run double pendulum
g = 9.8 #m/s^2

m1 = 1 #kg
m2 = 1 #kg

l1 = 1 #m
l2 = 1 #m

theta1 = np.radians(90)
theta2 = np.radians(45)

theta1_velocity = 0 #m/s
theta2_velocity = 0 #m/s

theta1_list = [theta1]
theta2_list = [theta2]

time_step = 20/300

time_span = np.linspace(0,20,300)
for t in time_span:
theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)
theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)
theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list)
y1 = -l1*np.cos(theta1_list)

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list)
y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list)

#Make Gif
!mkdir frames

counter=0
images = []
for i in range(0,len(x1)):
plt.figure(figsize = (6,6))

plt.figure(figsize = (6,6))
plt.plot([0,x1[i]],[0,y1[i]], "o-", color = "b", markersize = 7, linewidth=.7 )
plt.plot([x1[i],x2[i]],[y1[i],y2[i]], "o-", color = "b", markersize = 7, linewidth=.7 )
plt.title("Double Pendulum")
plt.xlim(-2.1,2.1)
plt.ylim(-2.1,2.1)
plt.savefig("frames/" + str(counter)+ ".png")
images.append(imageio.imread("frames/" + str(counter)+ ".png"))
counter += 1
plt.close()

imageio.mimsave("double_pendulum.gif", images)

!rm -r frames

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *