阻尼谐振器

无激振力

问题表述

\[\ddot{x}=-2\zeta\omega_0\点{x}-\omega_0^2 x\]

我们定义

\[\begin{split}\左\{ \BEGIN{array}{ll} x_1=x\\ x_2=\点{x} \end{数组} \对。\end{split}\]

我们是这样写的

\[\begin{split}\点{X}(T)=\BEGIN{pMatrix}\点{x}\\ -2\zeta\omega_0\点{x}-\omega_0^2 x\end{pMatrix}=f\Left(X，t\Right)\end{split}\]

%matplotlib notebook
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import integrate
import ipywidgets as ipw

阶跃响应

def ode(X, t, zeta, omega0):
    """
    Free Harmonic Oscillator ODE
    """
    x, dotx = X
    ddotx = -2*zeta*omega0*dotx - omega0**2*x
    return [dotx, ddotx]

def update(zeta = 0.05, omega0 = 2.*np.pi):
    """
    Update function.
    """
    X0 = [1., 0.]
    sol = integrate.odeint(ode, X0, t, args = (zeta, omega0))
    line0.set_ydata(sol[:, 0])

Nt = 1000
t = np.linspace(0., 10., Nt)
dummy = np.zeros_like(t)
fig = plt.figure()
line0, = plt.plot(t, dummy, label = "position")
plt.grid()
plt.ylim(-1., 1.)
plt.xlabel("Time, $t$")
plt.ylabel("Amplitude, $a$")
plt.legend()

ipw.interact(update, zeta = (0., 1., 0.01),
           omega0 = (2.*np.pi*0.05, 2.*np.pi*5, 2.*np.pi*0.01));

<IPython.core.display.Javascript object>

interactive(children=(FloatSlider(value=0.05, description='zeta', max=1.0, step=0.01), FloatSlider(value=6.283…

正弦激振力

问题表述

\[\ddot{x}=-2\zeta\omega_0\dot{x}-\omega_0^2x+F_m\sin(\omega_d t)\]

就像上面一样，我们可以把这个方程式写成适用于所有人 \(t\) ，

\[\begin{split}\点{X}(T)=\BEGIN{pMatrix}\点{x}\\ -2\zeta\omega_0\点{x}-\omega_0^2 x+F_m\sin(Omega_D T)\end{pMatrix}=f\Left(X，t\Right)\end{split}\]

def odeDrive(X, t, zeta, omega0, omegad_omega0):
    """
    Driven Harmonic Oscillator ODE
    """
    x, dotx = X
    omegad = omegad_omega0 * omega0
    ddotx = -2*zeta*omega0*dotx - omega0**2*x + F_m * np.sin(omegad * t)
    return [dotx, ddotx]

def update(zeta = 0.05, omega0 = 2.*np.pi, omegad_omega0 = 1.):
    """
    Update function.
    """
    X0 = np.zeros(2)
    sol = integrate.odeint(odeDrive, X0, t, args = (zeta, omega0, omegad_omega0))
    line0.set_ydata(sol[:, 0])

Nt = 1000
F_m = 1.
t = np.linspace(0., 10., Nt)
dummy = np.zeros_like(t)
fig = plt.figure()
line0, = plt.plot(t, dummy, label = "position")
plt.grid()
plt.ylim(-1., 1.)
plt.xlabel("Time, $t$")
plt.ylabel("Amplitude, $a$")
plt.legend()

ipw.interact(update, zeta = (0., .2, 0.01),
             omega0 = (2.*np.pi*0.5, 2.*np.pi*5, 2.*np.pi*0.01),
             omegad_omega0 = (0.1, 2., 0.05));

<IPython.core.display.Javascript object>

interactive(children=(FloatSlider(value=0.05, description='zeta', max=0.2, step=0.01), FloatSlider(value=6.283…

无激励相图

omega0 = 2. * np.pi * 1.
zeta   = .1
ad = 1.
omegad = 2. * np.pi * 1.2
Nt = 1000
tf = 10.
t = np.linspace(0., tf, Nt+1)
X0 = [1., 0.]
X = integrate.odeint(ode, X0, t, args = (zeta, omega0))

plt.figure()
plt.title("Phase plane with $\zeta =$ "+str(zeta))
plt.plot(X[:,0], X[:,1])
plt.grid()
plt.xlabel("$x$, [$m$]")
plt.ylabel("$\dot{x}$, [$m.s^{-1}$] ")
plt.show();

<IPython.core.display.Javascript object>

##带激励的相图

omegad_omega0 = 1.
X = integrate.odeint(odeDrive, X0, t, args = (zeta, omega0, omegad_omega0))

plt.figure()
plt.title("Phase plane with $\zeta =$ "+str(zeta))
plt.plot(X[:,0], X[:,1])
plt.grid()
plt.xlabel("$x$, [$m$]")
plt.ylabel("$\dot{x}$, [$m.s^{-1}$] ")
plt.show();

<IPython.core.display.Javascript object>