DWH.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Apr 17 15:53:42 2019

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

Biased Double-Well Potential
"""

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

plt.close('all')
T = 400
Amp = 3.5

def solve_flow(y0,c0,lim = [-3,3,-3,3]):

    def flow_deriv(x_y, t, c0):
        #"""Compute the time-derivative of a Medio system."""
        x, y = x_y
        #window = signal.triang(T)
        return [y,-0.5*y - x**3 + 2*x + x*(2*np.pi/T)*Amp*np.cos(2*np.pi*t/T) + Amp*np.sin(2*np.pi*t/T)]
        #return [y,-0.33*y - 2*x]
        #return [y,-0.99*y - x**3 + 2*x + (2*np.pi/T)*24*signal.triang(t/T)]



#    tt = np.zeros(shape=(tloopmax,))
#    xtt = np.zeros(shape=(tloopmax,))
#    cc = np.zeros(shape=(tloopmax,))
#    for tloop in range(0,tloopmax):
#        
#    
#        tlo = (tloop-1)*delt
#        thi = tloop*delt
#    
#        if tloop < tloopmax/2:
#            c = c0 + 3*np.abs(c0)*tloop/(tloopmax/2)
#        else:
#            c = c0 + 3*np.abs(c0) - 3*np.abs(c0)*(tloop-tloopmax/2)/(tloopmax/2)
#            
#    
#
## Solve for the trajectories
#        t = np.linspace(tlo, thi, 11)
#        x_t = integrate.odeint(flow_deriv, y0, t, args=(c,))
#        
#        szt, dum = np.shape(x_t)
#        tt[tloop], xtt[tloop] = x_t[szt-1]
#        cc[tloop] = c
    tsettle = np.linspace(0,T,101)   
    yinit = y0;
    x_tsettle = integrate.odeint(flow_deriv,yinit,tsettle,args=(T,))
    y0 = x_tsettle[100,:]
    
    t = np.linspace(0, 1.5*T, 2001)
    x_t = integrate.odeint(flow_deriv, y0, t, args=(T,))
    c  = Amp*np.sin(2*np.pi*t/T)
        
    return t, x_t, c

eps = 0.0001
xc = np.zeros(shape=(100,))
X = np.zeros(shape=(100,))
Y = np.zeros(shape=(100,))
Z = np.zeros(shape=(100,))
for loop in range(0,100):
    c = -1.2 + 2.4*loop/100 + eps
    xc[loop]=c
    
    coeff = [-1, 0, 2, c]
    y = np.roots(coeff)
    
    xtmp = np.real(y[0])
    ytmp = np.real(y[1])
    
    X[loop] = np.min([xtmp,ytmp])
    Y[loop] = np.max([xtmp,ytmp])
    Z[loop]= np.real(y[2])


plt.figure(1)
lines = plt.plot(xc,X,xc,Y,xc,Z)
plt.setp(lines, linewidth=0.5)
plt.show()
plt.title('Roots')

y0 = [1.9, 0]
c0 = -2.


t, x_t, c = solve_flow(y0,c0)
y1 = x_t[:,0]
y2 = x_t[:,1]

plt.figure(2)
lines = plt.plot(t,y1)
plt.setp(lines, linewidth=0.5)
plt.show()
plt.ylabel('X Position')
plt.xlabel('Time')

plt.figure(3)
lines = plt.plot(c,y1)
plt.setp(lines, linewidth=0.5)
plt.show()
plt.ylabel('X Position')
plt.xlabel('Control Parameter')
plt.title('Hysteresis Figure')