$p$ is momentum (action variable)
$\phi$ is an angle
$t$ is time
$H(\phi,p,t)$ time dependent Hamiltonian
$b$ drift rate
$\epsilon$ primary resonance term, resonance strength. This resonance has a frequency $\omega_0 = \sqrt{\epsilon}$
$\mu$ secondary resonance strength.
$\nu$ perturbation frequency
We will be mapping at the period $T = 2 \pi/\nu$
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
state vector
y = [phi,p] dimension 2
Hamiltonian
H(phi,p,t) = p^2/2 - b*t*p - epsilon * cos(phi) - mu*cos(phi - nu t)
Change of state vector
dy/dt = [d phi/dt, dp/dt]
= [dH/dp, -dH/dphi] Hamilton's equations
= [p - b*t, -epsilon*sin(phi) - mu*sin(phi - nu*t)]
# dy/dt = func(y,t)
# here are the equations of motion
def func(y,t,epsilon,mu,nu,b):
return [y[1] - b*t,-epsilon*np.sin(y[0]) - mu*np.sin(y[0] - nu*t)]
# the jacobian Dfunc, you can pass this to integrator if you want to
def jacobian(y,t,epsilon,mu,nu,b):
return [[0,1],[-epsilon*np.cos(y[0]) - mu*np.cos(y[0]-nu*t),0]]
# order for jacobian is [[dy0dt/dy0,dy0dt/dy1],[dy1dt/dy0,dy1dt/dy1]
twopi = 2.0*np.pi # global
# integrate one period (T = 2pi/nu)
# return y (our statevector) and new time after the integration
# epsilon, mu, nu, b are parameters for the Hamiltonian model
# here y0 is initial conditions
# t0 is initial time
# calls numpy routine odeint for integration
# passes jacobian() and func() (should be defined above)
def one_step(y0,t0,epsilon,mu,nu,b):
step = 2.0*np.pi/nu # mapping at this period!
t1 = t0 + step # final time
time = [t0,t1]
# now do the integration from t0 to t1
yvec = odeint(func, y0, time, Dfun=jacobian, args=(epsilon,mu,nu,b)
,rtol=1e-8,atol=1e-8,mxstep=10000000)
# if you ask for full_output=1 it appends that to yvec!
p1=yvec[1][1]
phi1=yvec[1][0]%twopi
# print("one_step phi,p,tn = ",phi1,p1,t1/step) # testing!
y1 =np.array([phi1,p1])
return y1,t1 # return new vector y and time
# integrate nperiods
# return arrays of points (that can then be plotted!)
# y0 is initial conditions, [phi0,p0]
# t0 is initial time
# epsilon, mu, nu, b are parameters of Hamiltonian model
# calls routine one_step()
def retpoints(y0,t0,nperiods,epsilon,mu,nu,b):
phivec = np.zeros(0)
pvec = np.zeros(0)
tvec = np.zeros(0)
y=y0; t=t0;
phivec = np.append(phivec,y[0])
pvec = np.append(pvec,y[1])
tvec = np.append(tvec,t)
for i in range(nperiods): # integrating in steps
y1,t1=one_step(y,t,epsilon,mu,nu,b)
y=y1
t=t1
phivec = np.append(phivec,y[0])
pvec = np.append(pvec,y[1])
tvec = np.append(tvec,t)
return phivec,pvec,tvec
# calling above but now plotting points!
# npoints is number of points plotted
# y0 is initial conditions, [phi0,p0]
# t0 is initial time
# marker lets you chose the color to plot!
# epsilon, mu, nu, b are parameters of Hamiltonian model
# do we correct for drift?, not important if b=0
def pltpts(y0,t0,npoints,marker,epsilon,mu,nu,b):
phivec,pvec,tvec = retpoints(y0,t0,npoints,epsilon,mu,nu,b)
pdrift=pvec - b*tvec # correct for drift when plotting!
plt.plot(phivec,pdrift,marker,markersize=2) # plot it
# set up plot
plt.xlabel(r'$\phi$',fontsize=22); plt.ylabel(r'$p$',fontsize=22)
plt.xlim((0,twopi))
# parameters for hamiltonian model
nu=1.0
epsilon=-0.5 # resonance #1
mu=-0.05 # resonance #2
npoints=500 # number of points to plot
b=0.00 # no drift
t0=0.0 #initial time
# do the integrations and plot them
# here first parameter is initial conditions for integrations
y0=[np.pi,0.2]; pltpts(y0,t0,npoints,'y.',epsilon,mu,nu,b)
y0=[np.pi,0.4]; pltpts(y0,t0,npoints,'r.',epsilon,mu,nu,b)
y0=[np.pi,0.8]; pltpts(y0,t0,npoints,'g.',epsilon,mu,nu,b)
y0=[np.pi,1.0]; pltpts(y0,t0,npoints,'b.',epsilon,mu,nu,b)
y0=[np.pi,1.05]; pltpts(y0,t0,npoints,'k.',epsilon,mu,nu,b)
y0=[np.pi,1.2]; pltpts(y0,t0,npoints,'b.',epsilon,mu,nu,b)
y0=[np.pi,1.5]; pltpts(y0,t0,npoints,'m.',epsilon,mu,nu,b)
y0=[np.pi,2.0]; pltpts(y0,t0,npoints,'k.',epsilon,mu,nu,b)
y0=[np.pi,-1.7]; pltpts(y0,t0,npoints,'y.',epsilon,mu,nu,b)
y0=[np.pi,-2.0]; pltpts(y0,t0,npoints,'r.',epsilon,mu,nu,b)
y0=[np.pi,-1.5]; pltpts(y0,t0,npoints,'k.',epsilon,mu,nu,b)
y0=[np.pi,-1.3]; pltpts(y0,t0,npoints,'m.',epsilon,mu,nu,b)
y0=[np.pi,0.5]; pltpts(y0,t0,npoints,'c.',epsilon,mu,nu,b)
y0=[1.0,-0.1]; pltpts(y0,t0,npoints,'y.',epsilon,mu,nu,b)
y0=[0.0,0.8]; pltpts(y0,t0,npoints,'r.',epsilon,mu,nu,b)
y0=[0.0,-0.2]; pltpts(y0,t0,npoints,'g.',epsilon,mu,nu,b)