DampedDriven.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 21 06:03:32 2018

@author: David Nolte
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)

Damped-driven pendulum
"""

import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt


plt.close('all')
print('DampedDriven.py')

# model_case 1 = Pendulum
# model_case 2 = Double Well
print(' ')
print('DampedDriven.py')
print('Case: 1 = Pendulum 2 = Double Well')
model_case = int(input('Enter the Model Case (1-2)'))

if model_case == 1:
    F = 0.6          # 0.6
    delt = 0.25      # 0.1
    w = 0.7          # 0.7
    def flow_deriv(x_y_z,tspan):
        x, y, z = x_y_z
        a = y
        b = F*np.cos(w*tspan) - np.sin(x) - delt*y
        c = w
        return[a,b,c]
else:
    alpha = -1       # -1
    beta = 1         # 1
    delta = 0.3       # 0.3
    gam = 0.15    # 0.15
    w = 1
    def flow_deriv(x_y_z,tspan):
        x, y, z = x_y_z
        a = y
        b = delta*np.cos(w*tspan) - alpha*x - beta*x**3 - gam*y
        c = w
        return[a,b,c]
                
T = 2*np.pi/w

# initial conditions
px1 = .1
xp1 = .1
w1 = 0

x_y_z = [xp1, px1, w1]

# Settle-down Solve for the trajectories
t = np.linspace(0, 2000, 40000)
x_t = integrate.odeint(flow_deriv, x_y_z, t)
x0 = x_t[39999,0:3]

tspan = np.linspace(1,20000,400000)
x_t = integrate.odeint(flow_deriv, x0, tspan)
siztmp = np.shape(x_t)
siz = siztmp[0]

if model_case == 1:
    y1 = np.mod(x_t[:,0]-np.pi,2*np.pi)-np.pi
    y2 = x_t[:,1]
    y3 = x_t[:,2]    
else:      
    y1 = x_t[:,0]
    y2 = x_t[:,1]
    y3 = x_t[:,2]
    
plt.figure(1)
lines = plt.plot(y1,y2,'ko',ms=1)
plt.setp(lines, linewidth=0.5)
plt.show()

repnum = 5000
px = np.zeros(shape=(2*repnum,))
xvar = np.zeros(shape=(2*repnum,))
cnt = -1
testwt = np.mod(tspan,T)-0.5*T;
last = testwt[1]
for loop in range(2,siz):
    if (last < 0)and(testwt[loop] > 0): # check for trajectory intersection with Poincare section
        cnt = cnt+1
        del1 = -testwt[loop-1]/(testwt[loop] - testwt[loop-1])
        px[cnt] = (y2[loop]-y2[loop-1])*del1 + y2[loop-1]
        xvar[cnt] = (y1[loop]-y1[loop-1])*del1 + y1[loop-1]
        last = testwt[loop]
    else:
        last = testwt[loop]
 
plt.figure(2)
lines = plt.plot(xvar,px,'ko',ms=1)
plt.show()


if model_case == 1:
    plt.savefig('DrivenPendulum')
else:
    plt.savefig('DrivenDoubleWell')