Oscilador Armónico Unidimensional en Python

El movimiento oscilatorio armónico es uno de los modelos físicos en el cual una partícula oscila en torno a una posición de equilibrio. En este caso se va analizar el movimiento oscilatorio armónico unidimensional en los casos simple, amortiguado y forzado. La ecuación que analiza este fenómeno físico tan particular es la siguiente:

$${\large \begin{equation} \frac{d^2 x}{dt^2} + \frac{b}{m}\frac{dx}{dt} + \frac{k}{m}x = \frac{F}{m}\sin{(w_{f}t)} \end{equation} }$$

Donde:

m es la masa de la partícula

b es la constante de amortiguamiento

k es la constante de elasticidad

F es la fuerza oscilante

$w_{f}$ es la frecuencia oscilante

Antes de empezar a simular el movimiento del oscilador armónico en Python debemos importar algunas librerías de antemano.

import matplotlib.pyplot as plt

import numpy as np

from scipy.integrate import odeint

import matplotlib.patches as mpatches

Ahora necesitamos definir una función que simule la ecuación del movimiento del oscilador armónico en Python.

def mps(p0, m, k, b, F, w):

    

    plt.figure(figsize=(14,8))

    

    def df1(f1,t):

        x1, v1 = f1[0], f1[1] 

        dx1 = v1  

        dv1 = (F/m)*np.sin(w*t) - (b/m)*v1 - (k/m)*x1 

        return [dx1,dv1]

    

    def df2(f2,t):

        x2, v2 = f2[0], f2[1] 

        dx2 = v2  

        dv2 = - (b/m)*v2 - (k/m)*x2 

        return [dx2,dv2]

    

    def df3(f3,t):

        x3, v3 = f3[0], f3[1] 

        dx3 = v3  

        dv3 = - (k/m)*x3 

        return [dx3,dv3]

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

    sol1 = odeint(df1, p0, t)

    sol2 = odeint(df2, p0, t)

    sol3 = odeint(df3, p0, t)

    f1 = sol1[:,0]

    f2 = sol2[:,0]

    f3 = sol3[:,0]

    

    plt.grid()

    plt.plot(t,f1, color = 'black', lw = 3, label = "Forzada")

    plt.plot(t,f2, color = 'blue', lw = 3, label = "Amortiguada")

    plt.plot(t,f3, color = 'red', lw = 3, label = "Simple")

    plt.title("OSCILACIONES", size = 30)

    plt.xlabel("TIEMPO", size = 20)

    plt.ylabel("POSICIÓN", size = 20)

    plt.legend(fontsize = 20, bbox_to_anchor = (1, 1))

    plt.show()

    return mps

Si aplicamos el código para las siguientes condiciones iniciales con los valores de las siguientes constantes.

mps([0, 4], 1, 40, 0.4, 10, 4)

 

A continuación una animación del modelo físico con las mismas condiciones iniciales