Perturbed.py

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

@author: nolte
"""
from IPython import get_ipython
get_ipython().magic('reset -f') 

import numpy as np
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from scipy import integrate
from matplotlib import pyplot as plt
from matplotlib import cm
import time
import os


plt.close('all')


# 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.02          # 0.6
    delt = 0.0       # 0.1
    w = 3/4          # 0.7
    k = 2
    phase = 0
    px1 = 1.9635
    xp1 = 0
    w1 = 0
    def flow_deriv(x_y_z,tspan):
        x, y, z = x_y_z
        a = y
        b = F*np.cos(-w*tspan + k*x + phase) - np.sin(x) - delt*y
        c = w
        return[a,b,c]
        
else:
    alpha = -1       # -1
    beta = 1         # 1
    F = 0.002       # 0.3
    delta = 0.0    # 0.15
    w = 1
    k = 1
    phase = np.random.random()
    px1 = 0
    xp1 = 0
    w1 = 0
    def flow_deriv(x_y_z,tspan):
        x, y, z = x_y_z
        a = y
        b = F*np.cos(-w*tspan + k*x + phase) - alpha*x - beta*x**3 - delta*y
        c = w
        return[a,b,c]
                

T = 2*np.pi/w

x_y_z = [xp1, px1, w1]

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

tlim = 200000  # number of points
nt = 40000   #stop time
tspan = np.linspace(1,nt,tlim)
x_t = integrate.odeint(flow_deriv, x0, tspan, rtol=1e-8)
siztmp = np.shape(x_t)
siz = siztmp[0]

y1 = np.zeros(shape=(2*tlim,))
y2 = np.zeros(shape=(2*tlim,))
if model_case == 1:
    y1tmp = np.mod(x_t[:,0]-np.pi,2*np.pi)-np.pi
    y2tmp = x_t[:,1]

    y1[0:tlim] = y1tmp
    y1[tlim:2*tlim] = y1tmp+2*np.pi
    y2[0:tlim] = y2tmp
    y2[tlim:2*tlim] = y2tmp

    y3 = x_t[:,2]    
    Energy = 0.5*x_t[:,1]**2 + 1 - np.cos(x_t[:,0])
else:      
    y1 = x_t[:,0]
    y2 = x_t[:,1]
    y3 = x_t[:,2]
    Energy = 0.5*x_t[:,1]**2 + 1 + 0.5*alpha*x_t[:,0]**2 + 0.25*beta*x_t[:,0]**4

    

plt.figure(1)
lines = plt.plot(y1,y2,'ko',ms=1)
plt.setp(lines, linewidth=0.5)
plt.title('Phase Portrait')
plt.show()


plt.figure(2)
lines = plt.plot(y3[0:3000],y2[0:3000])
plt.setp(lines, linewidth=0.5)
plt.title('Velocity')
plt.show()

plt.figure(3)
lines = plt.plot(y3[0:3000],Energy[0:3000])
plt.setp(lines, linewidth=0.5)
plt.title('Energy')
plt.show()


# First-Return Map
repnum = 5000
px = np.zeros(shape=(2*repnum,))
xvartmp = np.zeros(shape=(2*repnum,))
cnt = -1
testwt = np.mod(tspan,T)-0.5*T;
last = testwt[0]
for loop in range(1,siz):
    if (last < 0)and(testwt[loop] > 0):
        cnt = cnt+1
        del1 = -testwt[loop-1]/(testwt[loop] - testwt[loop-1])
        px[cnt] = (y2[loop]-y2[loop-1])*del1 + y2[loop-1]
        xvartmp[cnt] = (x_t[loop,0]-x_t[loop-1,0])*del1 + x_t[loop-1,0]
        #xvar[cnt] = y1[loop]
        last = testwt[loop]
    else:
        last = testwt[loop]
 


# Plot First Return Map
if model_case == 1:
    xvar = np.mod(xvartmp-np.pi,2*np.pi)-np.pi
    pxx = np.zeros(shape=(2*cnt,))
    xvarr = np.zeros(shape=(2*cnt,))
    
    xvarr[0:cnt] = xvar[0:cnt]
    xvarr[cnt:2*cnt] = xvar[0:cnt]+2*np.pi
    pxx[0:cnt] = px[0:cnt]
    pxx[cnt:2*cnt] = px[0:cnt]
      
    plt.figure(4)
    lines = plt.plot(xvarr,pxx,'ko',ms=0.5)
    plt.xlim(xmin=0, xmax=2*np.pi)
    plt.title('First Return Map')
    plt.show()
     
    plt.savefig('PPendulum')
else:
    xvar = xvartmp
    plt.figure(4)
    lines = plt.plot(xvar,px,'ko',ms=0.5)
    #mpl.pyplot.xlim(xmin=0, xmax=2*np.pi)
    plt.title('First Return Map')
    plt.show()    
    
    
    plt.savefig('PDoubleWell')