import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def f(y, t, params):
theta, omega = y
Q, d, Omega = params
derivs = [omega,
-omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
return derivs
Q = 2.0
d = 1.5
Omega = 0.65
theta0 = 0.0
omega0 = 0.0
params = [Q, d, Omega]
y0 = [theta0, omega0]
tStop = 200.
tInc = 0.05
t = np.arange(0., tStop, tInc)
psoln = odeint(f, y0, t, args=(params,))
fig = plt.figure(1, figsize=(8,8))
ax1 = fig.add_subplot(311)
ax1.plot(t, psoln[:,0])
ax1.set_xlabel('time')
ax1.set_ylabel('theta')
ax2 = fig.add_subplot(312)
ax2.plot(t, psoln[:,1])
ax2.set_xlabel('time')
ax2.set_ylabel('omega')
ax3 = fig.add_subplot(313)
twopi = 2.0*np.pi
ax3.plot(psoln[:,0]%twopi, psoln[:,1], '.', ms=1)
ax3.set_xlabel('theta')
ax3.set_ylabel('omega')
ax3.set_xlim(0., twopi)
plt.tight_layout()
plt.show()