Add files via upload

nolte · Jan 29, 2021 · 423fcb7316fec9a0262e115340115b3d09334319 · 423fcb7
1 parent 16617a1
commit 423fcb7316fec9a0262e115340115b3d09334319
diff --git a/DWH.py b/DWH.py
@@ -0,0 +1,114 @@
+#!/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')
diff --git a/DampedDriven.py b/DampedDriven.py
@@ -0,0 +1,107 @@
+#!/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')
diff --git a/Duffing.py b/Duffing.py
@@ -0,0 +1,98 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed May 21 06:03:32 2018
+
+@author: nolte
+Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)
+
+Duffing oscillator
+"""
+
+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('Duffing.py')
+
+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
+
+px1 = np.random.rand(1)
+xp1 = np.random.rand(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]
+
+y1 = x_t[:,0]
+y2 = x_t[:,1]
+y3 = x_t[:,2]
+
+plt.figure(2)
+lines = plt.plot(y1[1:2000],y2[1:2000],'ko',ms=1)
+plt.setp(lines, linewidth=0.5)
+plt.show()
+
+for cloop in range(0,3):
+
+#phase = np.random.rand(1)*np.pi;
+    phase = np.pi*cloop/3
+
+    repnum = 5000
+    px = np.zeros(shape=(2*repnum,))
+    xvar = np.zeros(shape=(2*repnum,))
+    cnt = -1
+    testwt = np.mod(tspan-phase,T)-0.5*T;
+    last = testwt[1]
+    for loop in range(2,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]
+            xvar[cnt] = (y1[loop]-y1[loop-1])*del1 + y1[loop-1]
+            last = testwt[loop]
+        else:
+            last = testwt[loop]
+
+    plt.figure(3)
+    if cloop == 0:
+        lines = plt.plot(xvar,px,'bo',ms=1)
+    elif cloop == 1:
+        lines = plt.plot(xvar,px,'go',ms=1)
+    else:
+        lines = plt.plot(xvar,px,'ro',ms=1)
+
+    plt.show()
+
+plt.savefig('Duffing')