import numpy as np
import matplotlib.pyplot as plt
import ipywidgets
from scipy import integrate, optimize

plt.rc("font", size=16)

Oscilador armónico

Ley de Hooke

Robert Hooke (1635-1703) fue un científico inglés, que, entre sus muchas contribuciones, formuló la ley de elasticidad o ley de Hooke.

Fig. 1: imagen de Hooke.

La ley dice que la fuerza elástica $F$ es proporcional al estiramiento $\Delta x$:

$$ |F| = k \; \Delta x $$

Experimento estático

La fuerza que hace el resorte tiene que cancelar al peso de la masa:

$$ |F_e| = k \Delta x = m g = |P| $$

Experimento dinámico

Planteando la ecuación de Newton para un cuerpo bajo el efecto de esta fuerza

$$ m \ddot{x} = -kx $$

se llega a la ecuación del oscilador armónico, generalmente escrita como:

$$ \ddot{x} + \omega^2 x = 0 $$

que da lugar a un movimiento armónico o sinusoidal:

$$ x(t) = A \cos(\omega t + \phi) $$

Apartamientos pequeños del equilibro

El movimiento armónico se observa siempre que tengamos un sistema en equilibrio (estable) y realicemos pequeños apartamientos del equilibrio.

En general, las fuerzas $F$ pueden tener dependencias complejas con la posición $x$:

$$ F = F(x) $$

x = np.linspace(-1, 1)
p = [-1, 0, -1, 0, -1, 0]
F = np.polyval(p, x)
V = -np.polyval(np.polyint(p), x)

fig, ax = plt.subplots(1, 2, figsize=(12, 4))

ax[0].plot(x, F)
ax[0].set_ylabel("Fuerza")

ax[1].plot(x, V)
ax[1].set_ylabel("Potencial")

for a in ax:
    a.set_xlabel("Posicion")
    a.axhline(0, color="k", ls="--")
    a.axvline(0, color="k", ls="--")
    a.grid()

Cerca del equilibrio (estable), cuando $F(x_0)=0$ (y $F'(x_0) < 0$), podemos aproximar la función por un desarrollo en serie de potencias (serie de Taylor):

$$ F(x) = F(x_0) + F'(x_0) \; (x - x_0) + F''(x_0) \; (x - x_0)^2 + \ldots $$

Y quedarnos hasta orden 1:

$$ F(x) \approx F'(x_0) \; (x - x_0) $$

que es la ley de Hooke, donde $F'(x_0)$ es la constante $k$, y $(x - x_0)$ es el desplazamiento $\Delta x$ respecto del equilibrio.

Péndulo

Masa puntual colgada de un hilo de largo $L$.

Se mueve a lo largo de una circunferencia parametrizada por el ángulo $\theta$.

!

Ecuacion de Newton:

$$ m L \ddot{\theta} = -m g \sin{\theta} $$

Reordenamos:

$$ \frac{d^2 \theta}{dt^2} + \frac{g}{L} \sin{\theta} = 0 $$

Aproximación ángulos pequeños ($\sin{\theta} \approx \theta$):

$$ \frac{d^2 \theta}{dt^2} + \frac{g}{L} \theta = 0 $$

que es la ecuación de un oscilador armónico:

$$ \frac{d^2 \theta}{dt^2} + \omega^2 \theta = 0 $$

Partiendo de un ángulo $\theta_0$ a velocidad $0$, la solución queda:

$$ \theta(t) = \theta_0 \, \cos(\omega t) $$

Vamos a analizar esta aproximación con datos simulados.

def pendulo_ode(t, y, a):
    """Ecuaciones diferenciales para el péndulo ideal.

    a = g / L
    """
    dy = np.empty(2)
    dy[0] = y[1]
    dy[1] = -a * np.sin(y[0])
    return dy


def pendulo(angulo, largo, t_max=10, tol=1e-6):
    sol = integrate.solve_ivp(
        pendulo_ode,
        t_span=(0, t_max),
        y0=(angulo, 0),
        args=(9.8 / largo,),
        atol=tol,
        rtol=tol,
    )
    return sol.t, sol.y[0]


@ipywidgets.interact(angulo_inicial=(10, 170, 10), T_inicial=(1.0, 5.0), ruido=(0, 10, 0.1))
def _(angulo_inicial=0.01, ruido=0, T_inicial=2):
    L = 1

    # Datos
    angulo = np.deg2rad(angulo_inicial)
    t, y = pendulo(angulo=angulo, largo=L)
    y = np.rad2deg(y)  # reconvierto a grados
    np.random.seed(42)
    y += np.random.normal(scale=ruido * np.abs(y), size=y.size)

    # Ajuste
    def func(t, A, w):
        return A * np.cos(w * t)

    popt, pcov = optimize.curve_fit(func, t, y, p0=(angulo_inicial, 2*np.pi/T_inicial))

    # Grafico
    fig, ax = plt.subplots(2, 1, figsize=(6, 6))

    ax[0].plot(t, y, "o--")  # Datos
    ax[0].plot(t, func(t, *popt))  # Ajuste

    ax[1].plot(t, y - func(t, *popt), "o--")  # Residuos
    ax[1].axhline(0, color="k", lw=2)

    A, w = popt
    T = 2 * np.pi / w
    fig.suptitle(f"{A=:.2f}\n{T=:.2f}")

Datasets

Trayectoria angular del péndulo para 3 angulos iniciales y 3 tipos de ruidos.

Angulos iniciales: 10°, 70° y 140°.

Tipos de ruido:

1. Normal o gaussiano, centrado en $\mu=0$ y con desviación estandar $\sigma=1$.
2. Uniforme entre $-0.5$ y $0.5$.
3. De redondeo: se redondeó a 0 decimales ($14.62..° \mapsto 15°$)

http://users.df.uba.ar/maurosilber/datos.zip

for angulo in (10, 70, 140):
    t, y = pendulo(angulo=np.deg2rad(angulo), largo=1)
    y = np.rad2deg(y)

    np.random.seed(42)
    np.savetxt(f"ideal/pendulo_angulo_{angulo}_ruido_normal.txt",
               np.array([t, np.random.normal(y, scale=1, size=y.size)]).T)

    np.random.seed(42)
    np.savetxt(f"ideal/pendulo_angulo_{angulo}_ruido_uniforme.txt",
               np.array([t, y + np.random.uniform(-0.5, 0.5, size=y.size)]).T)

    np.random.seed(42)
    np.savetxt(f"ideal/pendulo_angulo_{angulo}_ruido_redondeo.txt",
               np.array([t, np.round(y)]).T)

Guia

1. Analicen la dependencia del ajuste no lineal con los parametros iniciales (para angulo 10° y ruido normal)

2. Estimen los errores de los datos a partir de los residuos para los diferentes tipos de ruidos y angulos. Por ejemplo, para ruido normal les debería dar $\sigma = 1$.