import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

Péndulo

def oscilador(t, A, T):
    return A * np.cos(2 * np.pi * t / T)

The history saving thread hit an unexpected error (OperationalError('database is locked')).History will not be written to the database.

Dependencia con parámetros iniciales

t, y = np.loadtxt("ideal/pendulo_angulo_10_ruido_normal.txt", unpack=True)

Parametros buenos

params, _ = curve_fit(oscilador, t, y, p0=(10, 2))
plt.plot(t, y, "o")
plt.plot(t, oscilador(t, *params))

params

array([9.90817092, 2.0109836 ])

Periodo chico

params, _ = curve_fit(oscilador, t, y, p0=(10, 1))
plt.plot(t, y, "o")
plt.plot(t, oscilador(t, *params))

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

Amplitud chica

params, _ = curve_fit(oscilador, t, y, p0=(0.001, 2))
plt.plot(t, y, "o")
plt.plot(t, oscilador(t, *params))

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

Estimación de errores con los residuos

def figura(t, y, func, params):
    y_fit = func(t, *params)
    residuos = y - y_fit
    error = np.std(residuos)
    
    # Grafico
    fig = plt.figure(figsize=(10, 4))
    left, right = fig.add_gridspec(1, 2)

    ax = left.subgridspec(2, 1, height_ratios=(3, 2)).subplots(sharex=True)

    # Ajuste
    ax[0].plot(t, y, "o")
    ax[0].plot(t, y_fit)
    ax[0].set(title=f"A={params[0]:.3f}  T={params[1]:.3f}")

    # Residuos
    ax[1].axhline(0, color="k", ls="--")
    ax[1].scatter(t, residuos)

    # Histograma residuos
    ax_hist = fig.add_subplot(right)
    ax_hist.hist(residuos, bins="fd")
    ax_hist.set(title=f"$\sigma = {error:.2f}$")

Péndulo con ruido gaussiano

np.random.normal(0, 1)

0.2874378816893253

Angulo 10

t, y = np.loadtxt("ideal/pendulo_angulo_10_ruido_normal.txt", unpack=True)
params, _ = curve_fit(oscilador, t, y, p0=(10, 2))
figura(t, y, oscilador, params)

Angulo 70

t, y = np.loadtxt("ideal/pendulo_angulo_70_ruido_normal.txt", unpack=True)
params, _ = curve_fit(oscilador, t, y, p0=(70, 2))
figura(t, y, oscilador, params)

Angulo 140

t, y = np.loadtxt("ideal/pendulo_angulo_140_ruido_normal.txt", unpack=True)
params, _ = curve_fit(oscilador, t, y, p0=(140, 3))
figura(t, y, oscilador, params)

Péndulo con ruido uniforme

t, y = np.loadtxt("ideal/pendulo_angulo_10_ruido_uniforme.txt", unpack=True)
params, _ = curve_fit(oscilador, t, y, p0=(10, 2))
figura(t, y, oscilador, params)

1 / np.sqrt(12)

0.2886751345948129

Péndulo con ruido por redondeo

t, y = np.loadtxt("ideal/pendulo_angulo_10_ruido_redondeo.txt", unpack=True)
params, _ = curve_fit(oscilador, t, y, p0=(10, 2))
figura(t, y, oscilador, params)