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

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

Disipación proporcional a la velocidad

Conocido también como la ley de Stokes, es una aproximación de la fuerza de fricción para un cuerpo (esférico) que se mueve en un fluido.

$$ |F_v| = 6\pi \mu R v $$

donde $\mu$ es el coeficiente de viscosidad del fluido, $R$ es el radio de la esfera, y $v$ la velocidad.

Agrupando todo eso en una constante $b$, tenemos:

$$ |F_v| = b \, v $$

Si resolvemos la ecuación de Newton para un cuerpo bajo efectos de esta fuerza:

$$ m \ddot{x} = -b \dot{x} $$

llegamos a la siguiente solución:

$$ x(t) = v_0 \, \exp{(-bt)} + C $$

t = np.linspace(0, 10, 100)

for b in (1, 2, 3):
    x = 5 * np.exp(-b * t)
    plt.plot(t, x, label=b)

plt.xlabel("Tiempo t")
plt.ylabel("Posicion x")
plt.legend(title="b")

<matplotlib.legend.Legend at 0x7fe8187c8b50>

Una mejor manera de graficar exponenciales es en escala logaritmica. Si

$$ x(t) = A e^{bt} $$

Aplicando logaritmo a ambos lados:

$$ \log\Big( x(t) \Big) = \log(A) + b t = a + b t $$

queda una recta.

t = np.linspace(0, 10, 100)

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

ax[0].set(ylabel="x")
ax[1].set(ylabel="log(x)")
ax[2].set(yscale="log", ylabel="x")
for a in ax:
    a.set(xlabel="Tiempo")
    a.grid()

for b in (1, 2, 3):
    x = 5 * np.exp(-b * t)
    ax[0].plot(t, x)
    ax[1].plot(t, np.log(x))
    ax[2].plot(t, x)
    
fig.tight_layout()

A partir de la pendiente y la ordenada, podemos obtener $\log{v_0}$ y $b$:

$$ x(t) = v_0 \, \exp{(-bt)} + C $$

Ojo, para este tipo de gráficos, es importante que la constante sea $0$.

t = np.linspace(0, 5, 100)

@ipywidgets.interact(c=(-1, 1, 0.01))
def _(c):
    fig, ax = plt.subplots(1, 3, figsize=(12, 4))

    ax[0].axhline(0, color="k", ls="--")
    ax[0].set(ylabel="x", ylim=(-0.5, 6))
    ax[1].set(ylabel="log(x)")
    ax[2].set(yscale="log", ylabel="x")
    for a in ax:
        a.set(xlabel="Tiempo")
        a.grid()

    for b in (1, 2, 3):
        x = 5 * np.exp(-b * t) + c
        ax[0].plot(t, x)
        ax[1].plot(t, np.log(x))
        ax[2].plot(t, x)

    fig.tight_layout()

Oscilador armónico amortiguado

Planteando la ecuación de Newton, donde agregamos un término disipativo,

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

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

$$ \ddot{x} + 2\gamma \dot{x} + \omega_0^2 x = 0 $$

donde $b = 2 \gamma$, que da lugar a un movimiento sinusoidal amortiguado:

$$ x(t) = A \exp{(-\gamma t)} \cos(\omega t + \phi) $$

donde $\omega^2 = \omega_0^2 + \gamma^2$ (a chequear)

t = np.linspace(0, 100, 1000)
y = np.exp(-0.1 * t) * np.cos(t)
plt.plot(t, y)

[<matplotlib.lines.Line2D at 0x7f293eda0940>]

Pendulo con disipación en la aproximación de angulos pequeños

Es un oscilador armonico amortiguado.

def pendulo_ode(t, y, a, b):
    """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]) - 2 * b * y[1]
    return dy


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


@ipywidgets.interact(
    angulo_inicial=(10, 170, 10),
    disipacion=(0, 0.2, 0.01),
    ruido=(0, 10, 0.1),
    T_inicial=(1.0, 3.0),
    t_max=(10, 100, 10),
)
def _(angulo_inicial=0.01, disipacion=0, ruido=0, T_inicial=2, t_max=50):
    L = 1

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

    # Ajuste
    def func(t, A, T, b):
        return A * np.exp(-b * t) * np.cos(2 * np.pi * t / T)

    popt, pcov = optimize.curve_fit(
        func, t, y, p0=(angulo_inicial, T_inicial, disipacion)
    )

    # Grafico
    fig, ax = plt.subplots(2, 1, figsize=(12, 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, T, b = popt
    fig.suptitle(f"{A=:.2f}   {T=:.2f}   {b=:.2f}")

Datasets

Pendulo para 3 angulos y 4 disipaciones:

Angulos: 10°, 70°, 140°

Disipación: 0, 0.03, 0.06, 0.09 Hz

for angulo in (10, 70, 140):
    for disp in (0, 0.03, 0.06, 0.09):
        t, y = pendulo(angulo=np.deg2rad(angulo), disipacion=disp, largo=1, t_max=50)
        y = np.rad2deg(y)

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

Guia

1. Estimar el coeficiente de disipación y el periodo para ángulo pequeño.

2. Repetir para ángulo grande.