import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Our Hamiltonian $$H(\theta,p) = \frac{1}{2}p^2 - \epsilon \cos(\theta) $$
We plot $X=\theta$ and $Y=p$.
Hamilton's equations are $$ \dot \theta = \frac{\partial H}{\partial p}=p \qquad \dot p = -\frac{\partial H}{\partial \theta} = -\epsilon \sin\theta $$
# range of x and y grid
xmax = np.pi*2.0
ymax = 2
fac=1.01
# make a grid of x and y values, Y = dot X
X,Y = np.meshgrid(np.arange(-xmax,xmax,.1),np.arange(-ymax,ymax,.1) )
epsilon=0.4
H = 0.5*Y*Y - epsilon*np.cos(X) #here is the Hamiltonian
# Hamilton's equations give us a vector field U,V
U = Y
V = -epsilon*np.sin(X)
fig, ax = plt.subplots(1,1,figsize=(4,4))
ax.set_xlabel(r'$ \theta$')
ax.set_ylabel(r'd$\theta$/dt')
ax.set_xlim([-xmax*fac, xmax*fac])
ax.set_ylim([-ymax*fac, ymax*fac])
# plot the vector field here with either of the two commands below
#Q = plt.quiver(X,Y,U, V)
Q = ax.streamplot(X,Y,U, V,density=2)
