¿ Cómo se mueve un péndulo ?

Se dice que un sistema cualquiera, mecánico, eléctrico, neumático, etc., es un oscilador armónico si, cuando se deja en libertad fuera de su posición de equilibrio, vuelve hacia ella describiendo oscilaciones sinusoidales, o sinusoidales amortiguadas en torno a dicha posición estable.

https://es.wikipedia.org/wiki/Oscilador_armónico

Referencias:

En realidad esto es el estudio de oscilaciones.

Los sistemas mas sencillos a estudiar en oscilaciones son el sistema masa-resorte y el péndulo simple.

\begin{align} \ddot{x} + \omega_{0}^2 x &= 0, \quad \omega_{0} = \sqrt{\frac{k}{m}}\notag\\ \ddot{\theta} + \omega_{0}^{2}\, \theta &= 0, \quad\mbox{donde}\quad \omega_{0}^2 = \frac{g}{l} \end{align}

Sistema `masa-resorte`

La solución a este sistema masa-resortese explica en términos de la segunda ley de Newton. Para este caso, si la masa permanece constante y solo consideramos la dirección en $x$. Entonces,

\begin{equation} F = m\ddot{x} = m \frac{d\dot{x}}{dt} = m \frac{d^2x}{dt^2} \end{equation}

¿ Cuál es la fuerza? Ley de Hooke! \begin{equation} F = -k x, \quad k > 0 \end{equation} Vemos que la fuerza se opone al desplazamiento y su intensidad es proporcional al mismo. Y $k$ es la constante elástica o recuperadora del resorte.

Cuya solución se escribe como \begin{equation} x(t) = A \cos(\omega_{o} t) + B \sin(\omega_{o} t) \end{equation} Y su primera derivada (velocidad) sería \begin{equation} \dot{x}(t) = \omega_{0}[- A \sin(\omega_{0} t) + B\cos(\omega_{0}t)] \end{equation}

¿Cómo se ven las gráficas de x vs t, tanto para la posición y velocidad?

Esta instrucción es para que las gráficas aparezcan dentro de este entorno.



In [1]:

    
%matplotlib inline

Esta es la librería con todas las instrucciones para realizar gráficos.



In [2]:

    
import matplotlib.pyplot as plt



In [3]:

    
import matplotlib as mpl
label_size = 14
mpl.rcParams['xtick.labelsize'] = label_size 
mpl.rcParams['ytick.labelsize'] = label_size

Y esta es la librería con todas las funciones matemáticas necesarias.



In [4]:

    
import numpy as np



In [5]:

    
t = np.linspace(0, 50, 100)
plt.figure(figsize = (7, 4))
plt.plot(t, .5*np.cos(.5 *t) + .1 * np.sin(.5 *t), '-', lw = 1, ms = 4,
         label = '$x(t)$')
plt.plot(t, .5*.1*np.cos(.5*t) - .5*.5 * np.sin(.5*t), 'ro-', lw = 1, ms = 4,
         label = r'$\dot{x(t)}$')
plt.xlabel('$t$', fontsize = 20)
plt.show()



In [6]:

    
plt.figure(figsize = (7, 4))
plt.scatter(t, .5*np.cos(.5 *t) + .1 * np.sin(.5 *t), lw = 0, c = 'red',
            label = '$x(t)$')
plt.plot(t, .5*np.cos(.5 *t) + .1 * np.sin(.5 *t), 'r-', lw = 1)
plt.scatter(t, .5*.1*np.cos(.5*t) - .5*.5 * np.sin(.5*t), lw = 0, c = 'b',
            label = r'$\dot{x(t)}$')
plt.plot(t, .5*.1*np.cos(.5*t) - .5*.5 * np.sin(.5*t), 'b-', lw = 1)
plt.xlabel('$t$', fontsize = 20)
plt.legend(loc = 'best')
plt.show()

Y si consideramos un conjunto de frecuencias de oscilación, entonces



In [49]:

    
frecuencias = np.array([.1, .2 , .5, .6])
plt.figure(figsize = (7, 4))
for f in frecuencias:
    plt.plot(t, .5 * np.cos(f * t) + .1 * np.sin(f * t), '*-')
plt.xlabel('$t$', fontsize = 16)
plt.ylabel('$x(t)$', fontsize = 16)
plt.title('Oscilaciones', fontsize = 16)
plt.show()



In [9]:

    
t=np.linspace(0,50)

Estos colores, son el default de matplotlib, sin embargo existe otra librería dedicada, entre otras cosas, a la presentación de gráficos.



In [10]:

    
import seaborn as sns
sns.set(style='ticks', palette='Set2')



In [11]:

    
frecuencias = np.array([.1, .2 , .5, .6])
plt.figure(figsize = (7, 4))
for f in frecuencias:
    plt.plot(t, .5 * np.cos(f * t) + .1 * np.sin(f * t), 'o-', 
             label = '$\omega_0 = %s$'%f)
plt.xlabel('$t$', fontsize = 16)
plt.ylabel('$x(t)$', fontsize = 16)
plt.title('Oscilaciones', fontsize = 16)
plt.legend(loc='center left', bbox_to_anchor=(-.5, 0.9), prop={'size': 20})
plt.show()



In [11]:

    
frecuencias = np.array([.1, .2 , .5, .6])
simbolos = np.array(['o', 'D', '*', '.'])
plt.figure(figsize = (7, 4))
for indx, f in enumerate (frecuencias):
    plt.plot(t, .5 * np.cos(f * t) + .1 * np.sin(f * t), simbolos[indx], 
             label = '$\omega_0 = %s$'%f)
plt.xlabel('$t$', fontsize = 16)
plt.ylabel('$x(t)$', fontsize = 16)
plt.title('Oscilaciones', fontsize = 16)
plt.legend(loc='center left', bbox_to_anchor=(-.5, 0.9), prop={'size': 20})
plt.show()



In [12]:

    
frecuencias = np.array([.1, .2 , .5, .6])
    lineas = np.array(['-', '-.', ':', '.'])
    colores = np.array(['blue', 'black', 'magenta', 'orange'])
    plt.figure(figsize = (7, 4))
    for indx, f in enumerate (frecuencias):
        plt.plot(t, .5 * np.cos(f * t) + .1 * np.sin(f * t), lineas[indx], lw = 3.4,
                 color = colores[indx],
                 label = '$\omega_0 = %s$'%f)
    plt.xlabel('$t$', fontsize = 16)
    plt.ylabel('$x(t)$', fontsize = 16)
    plt.title('Oscilaciones', fontsize = 16)
    plt.legend(loc='center left', bbox_to_anchor=(-.5, 0.9), prop={'size': 20})
    plt.show()

Si queremos tener manipular un poco mas las cosas, hacemos uso de lo siguiente:



In [13]:

    
from ipywidgets import *



In [14]:

    
def masa_resorte(t = 0):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(.5 * np.cos(.5 * t) + .1 * np.sin(.5 * t), [0],  'ko', ms = 10)
    ax.set_xlim(xmin = -0.6, xmax = .6)
    ax.axvline(x=0, color = 'r')
    ax.axhline(y=0, color = 'grey', lw = 1)
    fig.canvas.draw()

interact(masa_resorte, t = (0, 50,.01));

La opción de arriba generalmente será lenta, así que lo recomendable es usar interact_manual.



In [15]:

    
def masa_resorte(t = 0):
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(.5 * np.cos(.5 * t) + .1 * np.sin(.5 * t), [0],  'ko', ms = 10)
    ax.set_xlim(xmin = -0.6, xmax = .6)
    ax.axvline(x=0, color = 'r')
    ax.axhline(y=0, color = 'grey', lw = 1)
    fig.canvas.draw()
interact_manual(masa_resorte, t = (0, 50,.01));

Péndulo simple

Ahora, si fijamos nuestra atención al movimiento de un péndulo simple (oscilaciones pequeñas). La ecuación diferencial a resolver tiene la misma forma. La diferencia más evidente es como hemos definido a $\omega_{0}$. Esto quiere decir que,

\begin{equation} \theta(t) = A\cos(\omega_{0} t) + B\sin(\omega_{0}t) \end{equation}

Si graficamos la ecuación de arriba vamos a encontrar un comportamiento muy similar al ya discutido anteriormente. Es por ello que ahora veremos el movimiento en el plano $xy$. Es decir,

\begin{align} x &= l \sin(\theta), \quad y = l \cos(\theta) \end{align}



In [16]:

    
def theta_t(a, b, g, l, t):
    omega_0 = np.sqrt(g/l)
    return a * np.cos(omega_0 * t) + b * np.sin(omega_0 * t)



In [32]:

    
def pendulo_simple(t = 0):
        fig = plt.figure(figsize = (5,5))
        ax = fig.add_subplot(1, 1, 1)
        x = 2 * np.sin(theta_t(.4, .6, 9.8, 2, t))
        y =  - 2 * np.cos(theta_t(.4, .6, 9.8, 2, t))
        ax.plot(x, y, 'ko', ms = 10)
        ax.plot([0], [0], 'rD')
        ax.plot([0, x ], [0, y], 'k-', lw = 1)
        ax.set_xlim(xmin = -2.2, xmax = 2.2)
        ax.set_ylim(ymin = -2.2, ymax = .2)
        fig.canvas.draw()
    interact_manual(pendulo_simple, t = (0, 10,.01));

Condiciones iniciales

Realmente lo que se tiene que resolver es,

\begin{equation} \theta(t) = \theta(0) \cos(\omega_{0} t) + \frac{\dot{\theta}(0)}{\omega_{0}} \sin(\omega_{0} t) \end{equation}

Actividad. Modificar el programa anterior para incorporar las condiciones iniciales.



In [41]:

    
# Solución: 
def theta_t(a,t):
    omega_0 = np.sqrt(g/l)
    return theta(a) * np.cos(omega_0 * t) + ((theta(a)/ omega_0)* np.sin(omega_0 * t) )



In [42]:

    
def pendulo_simple(t = 0):
    fig = plt.figure(figsize = (5,5))
    ax = fig.add_subplot(1, 1, 1)
    x = 2 * np.sin(theta_t( 1, t))
    y =  - 2 * np.cos(theta_t(1, t))
    ax.plot(x, y, 'ko', ms = 10)
    ax.plot([0], [0], 'rD')
    ax.plot([0, x ], [0, y], 'k-', lw = 1)
    ax.set_xlim(xmin = -2.2, xmax = 2.2)
    ax.set_ylim(ymin = -2.2, ymax = .2)
    fig.canvas.draw()
interact_manual(pendulo_simple, t = (0, 10,.01));

Espacio fase $(q, p)$

Recuerden el momento lineal se define como $ p = m v, $ donde $m$ es la masa y $v$ la velocidad. Y la posición y velocidad para el sistema masa-resorte se escriben como:

\begin{align} x(t) &= x(0) \cos(\omega_{o} t) + \frac{\dot{x}(0)}{\omega_{0}} \sin(\omega_{o} t)\\ \dot{x}(t) &= -\omega_{0}x(0) \sin(\omega_{0} t) + \dot{x}(0)\cos(\omega_{0}t)] \end{align}



In [22]:

    
k = 3  #constante elástica [N]/[m] 
m = 1 # [kg]
omega_0 = np.sqrt(k/m)
x_0 = .5
x_0_dot = .1



In [23]:

    
t = np.linspace(0, 50, 300)



In [24]:

    
x_t = x_0 *np.cos(omega_0 *t) + (x_0_dot/omega_0) * np.sin(omega_0 *t)
x_t_dot = -omega_0 * x_0 * np.sin(omega_0 * t) + x_0_dot * np.cos(omega_0 * t)



In [47]:

    
plt.figure(figsize = (7, 4))
plt.plot(t, x_t, label = '$x(t)$', lw = 1)
plt.plot(t, x_t_dot, label = '$\dot{x}(t)$', lw = 1)
plt.legend(loc='center left', bbox_to_anchor=(1.01, 0.5), prop={'size': 14})
plt.xlabel('$t$', fontsize = 18)
plt.show()



In [26]:

    
plt.figure(figsize = (5, 5))
plt.plot(x_t, x_t_dot/omega_0, 'ro', ms = 2)
plt.xlabel('$x(t)$', fontsize = 18)
plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
plt.show()



In [27]:

    
plt.figure(figsize = (5, 5))
plt.scatter(x_t,  x_t_dot/omega_0, cmap = 'viridis', c = x_t_dot, s = 8, lw = 0)
plt.xlabel('$x(t)$', fontsize = 18)
plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
plt.show()

Multiples condiciones iniciales



In [28]:

    
k = 3  #constante elástica [N]/[m] 
m = 1 # [kg]
omega_0 = np.sqrt(k/m)



In [29]:

    
t = np.linspace(0, 50, 50)



In [30]:

    
x_0s = np.array([.7, .5, .25, .1])
x_0s_dot = np.array([.2, .1, .05, .01])
cmaps = np.array(['viridis', 'inferno', 'magma', 'plasma'])



In [31]:

    
plt.figure(figsize = (6, 6))
for indx, x_0 in enumerate(x_0s):
    x_t = x_0 *np.cos(omega_0 *t) + (x_0s_dot[indx]/omega_0) * np.sin(omega_0 *t)
    x_t_dot = -omega_0 * x_0 * np.sin(omega_0 * t) + x_0s_dot[indx] * np.cos(omega_0 * t)
    plt.scatter(x_t,  x_t_dot/omega_0, cmap = cmaps[indx], 
                c = x_t_dot, s = 10, 
                lw = 0)
    plt.xlabel('$x(t)$', fontsize = 18)
    plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
    #plt.legend(loc='center left', bbox_to_anchor=(1.05, 0.5))

Trayectorias del oscilador armónico simple en el espacio fase $(x,\, \dot{x}\,/\omega_0)$ para diferentes valores de la energía.



In [ ]:

Tarea, considerar un pendulo simle em diferentes planetas, y graficar su espacio fase



In [50]:

    
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint 
%matplotlib inline



In [51]:

    
from ipywidgets import *



In [52]:

    
def theta_t(a, b, g, l, t):
    omega_0 = np.sqrt(g/l)
    return a * np.cos(omega_0 * t) + b * np.sin(omega_0 * t)



In [53]:

    
def theta_t(a, b, g, l, t):
    omega_0 = np.sqrt(g/l)
    return a * np.cos(omega_0 * t) + b * np.sin(omega_0 * t)



In [54]:

    
def pendulo_simple(t = 0):
    fig = plt.figure(figsize = (5,5))
    ax = fig.add_subplot(1, 1, 1)
    x = 2 * np.sin(theta_t(.4, .6, 9.8, 2, t))
    y =  - 2 * np.cos(theta_t(.4, .6, 9.8, 2, t))
    ax.plot(x, y, 'ko', ms = 10)
    ax.plot([0], [0], 'rD')
    ax.plot([0, x ], [0, y], 'k-', lw = 1)
    ax.set_xlim(xmin = -2.2, xmax = 2.2)
    ax.set_ylim(ymin = -2.2, ymax = .2)
    fig.canvas.draw()
interact_manual(pendulo_simple, t = (0, 10,.01));



In [55]:

    
colores = np.array(['viridis','inferno','Blues','Greys','Reds','RdPu'])



In [56]:

    
gravedades = np.array([3.7, 8.87, 9.8, 3.7, 24.8, 10.44, 8.7, 8.87, 11.15])
plt.figure(figsize = (20, 30))
for gv in gravedades:
    plt.plot(2 * np.sin(theta_t(.4, .6, gv, 2, t)), - 2 * np.cos(theta_t(.4, .6, gv, 2, t)), '-', lw = 3.4,)
plt.title('Pendulos con diferentes gravedades', fontsize = 16)
plt.show()



In [57]:

    
k = 3  #constante elástica [N]/[m] 
m = 1 # [kg]
omega_0 = np.sqrt(k/m)



In [58]:

    
t = np.linspace(0, 50, 50)



In [59]:

    
x_0s = np.array([.7, .5, .25, .1])
x_0s_dot = np.array([.2, .1, .05, .01])
cmaps = np.array(['viridis', 'inferno', 'magma', 'plasma'])



In [60]:

    
plt.figure(figsize = (6, 6))
for indx, x_0 in enumerate(x_0s):
    x_t = x_0 *np.cos(omega_0 *t) + (x_0s_dot[indx]/omega_0) * np.sin(omega_0 *t)
    x_t_dot = -omega_0 * x_0 * np.sin(omega_0 * t) + x_0s_dot[indx] * np.cos(omega_0 * t)
    plt.scatter(x_t,  x_t_dot/omega_0, cmap = cmaps[indx], 
                c = x_t_dot, s = 10, 
                lw = 0)
    plt.xlabel('$x(t)$', fontsize = 18)
    plt.ylabel('$\dot{x}(t)/\omega_0$', fontsize = 18)
    #plt.legend(loc='center left', bbox_to_anchor=(1.05, 0.5))



In [ ]: