Perturbed pendulum model

$$H(\phi,p,t) = p^2/2 - bt p - \epsilon \cos(\phi) - \mu \cos(\phi - \nu t)$$

$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$

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline

Equations of motion

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)]
In [2]:
# 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
In [3]:
# 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
In [4]:
# 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
In [5]:
# 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
In [6]:
# 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)
In [ ]: