diff --git a/DWH.py b/DWH.py new file mode 100644 index 0000000..035b6ab --- /dev/null +++ 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 new file mode 100644 index 0000000..334962b --- /dev/null +++ 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 new file mode 100644 index 0000000..3f7c62e --- /dev/null +++ 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') diff --git a/Flow2D.py b/Flow2D.py new file mode 100644 index 0000000..d2c607c --- /dev/null +++ b/Flow2D.py @@ -0,0 +1,153 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 16 07:38:57 2018 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +2D Flow examples: Medio, van der Pol, Fitzhugh-Nagumo +""" +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +# model_case 1 = Medio +# model_case 2 = vdP +# model_case 3 = Fitzhugh-Nagumo +model_case = int(input('Input Model Case (1-3)')) + +def solve_flow(param,lim = [-3,3,-3,3],max_time=10.0): + + if model_case == 1: +# Medio 2D flow + def flow_deriv(x_y, t0, a,b,c,alpha): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [a*y + b*x*(c - y**2),-x+alpha] + model_title = 'Medio Economics' + + elif model_case == 2: +# van der pol 2D flow + def flow_deriv(x_y, t0, alpha,beta): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y,-alpha*x+beta*(1-x**2)*y] + model_title = 'van der Pol Oscillator' + + else: +# Fitzhugh-Nagumo + def flow_deriv(x_y, t0, alpha, beta, gamma): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y-alpha,-gamma*x+beta*(1-y**2)*y] + model_title = 'Fitzhugh-Nagumo Neuron' + + + + plt.figure() + xmin = lim[0] + xmax = lim[1] + ymin = lim[2] + ymax = lim[3] + plt.axis([xmin, xmax, ymin, ymax]) + + + N=144 + colors = plt.cm.prism(np.linspace(0, 1, N)) + + x0 = np.zeros(shape=(N,2)) + ind = -1 + for i in range(0,12): + for j in range(0,12): + ind = ind + 1; + x0[ind,0] = ymin-1 + (ymax-ymin+2)*i/11 + x0[ind,1] = xmin-1 + (xmax-xmin+2)*j/11 + + # Solve for the trajectories + t = np.linspace(0, max_time, int(250*max_time)) + x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param) + for x0i in x0]) + + for i in range(N): + x, y = x_t[i,:,:].T + lines = plt.plot(x, y, '-', c=colors[i]) + plt.setp(lines, linewidth=1) + + plt.show() + plt.title(model_title) + plt.savefig('Flow2D') + + return t, x_t + + +def solve_flow2(param,max_time=20.0): + + if model_case == 1: +# Medio 2D flow + def flow_deriv(x_y, t0, a,b,c,alpha): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [a*y + b*x*(c - y**2),-x+alpha] + model_title = 'Medio Economics' + x0 = np.zeros(shape=(2,)) + x0[0] = 1 + x0[1] = 1 + elif model_case == 2: +# van der pol 2D flow + def flow_deriv(x_y, t0, alpha,beta): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y,-alpha*x+beta*(1-x**2)*y] + model_title = 'van der Pol Oscillator' + x0 = np.zeros(shape=(2,)) + x0[0] = 0 + x0[1] = 4.5 + + else: +# Fitzhugh-Nagumo + def flow_deriv(x_y, t0, alpha, beta, gamma): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y-alpha,-gamma*x+beta*(1-y**2)*y] + model_title = 'Fitzhugh-Nagumo Neuron' + x0 = np.zeros(shape=(2,)) + x0[0] = 1 + x0[1] = 1 + + + # Solve for the trajectories + t = np.linspace(0, max_time, int(250*max_time)) + x_t = integrate.odeint(flow_deriv, x0, t, param) + + + return t, x_t + + + + +if model_case == 1: + param = (0.9,0.7,0.5,0.6) # Medio + lim = (-7,7,-5,5) +elif model_case == 2: + param = (5, 2.5) # van der Pol + lim = (-7,7,-10,10) +else: + param = (0.02,0.5,0.2) # Fitzhugh-Nagumo + lim = (-7,7,-4,4) + +t, x_t = solve_flow(param,lim) + + +t, x_t = solve_flow2(param) +plt.figure(2) +lines = plt.plot(t,x_t[:,0],t,x_t[:,1],'-') + + + + + + + diff --git a/Flow2DBorder.py b/Flow2DBorder.py new file mode 100644 index 0000000..aa2b148 --- /dev/null +++ b/Flow2DBorder.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 16 07:38:57 2018 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Border initial conditions for 2D flow +""" +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +# model_case 1 = Medio +# model_case 2 = vdP +# model_case 3 = Fitzhugh-Nagumo +model_case = int(input('Input Model Case (1-3)')) + +def solve_flow(param,lim = [-3,3,-3,3],max_time=20.0): + + if model_case == 1: +# Medio 2D flow + def flow_deriv(x_y, t0, a,b,c,alpha): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [a*y + b*x*(c - y**2),-x+alpha] + model_title = 'Medio Economics' + + elif model_case == 2: +# van der pol 2D flow + def flow_deriv(x_y, t0, alpha,beta): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y,-alpha*x+beta*(1-x**2)*y] + model_title = 'van der Pol Oscillator' + + else: +# Fitzhugh-Nagumo + def flow_deriv(x_y, t0, alpha, beta, gamma): + #"""Compute the time-derivative of a Medio system.""" + x, y = x_y + return [y-alpha,-gamma*x+beta*(1-y**2)*y] + model_title = 'Fitzhugh-Nagumo Neuron' + + + + plt.figure() + xmin = lim[0] + xmax = lim[1] + ymin = lim[2] + ymax = lim[3] + plt.axis([xmin, xmax, ymin, ymax]) + + + N = 24*4 + 1 + x0 = np.zeros(shape=(N,2)) + ind = -1 + for i in range(0,24): + ind = ind + 1 + x0[ind,0] = xmin + (xmax-xmin)*i/23 + x0[ind,1] = ymin + ind = ind + 1 + x0[ind,0] = xmin + (xmax-xmin)*i/23 + x0[ind,1] = ymax + ind = ind + 1 + x0[ind,0] = xmin + x0[ind,1] = ymin + (ymax-ymin)*i/23 + ind = ind + 1 + x0[ind,0] = xmax + x0[ind,1] = ymin + (ymax-ymin)*i/23 + ind = ind + 1 + x0[ind,0] = 0.05 + x0[ind,1] = 0.05 + + colors = plt.cm.prism(np.linspace(0, 1, N)) + + + # Solve for the trajectories + t = np.linspace(0, max_time, int(250*max_time)) + x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param) + for x0i in x0]) + + for i in range(N): + x, y = x_t[i,:,:].T + lines = plt.plot(x, y, '-', c=colors[i]) + plt.setp(lines, linewidth=1) + + plt.show() + plt.title(model_title) + + return t, x_t + + + +if model_case == 1: + param = (0.9,0.7,0.5,0.6) # Medio + lim = (-7,7,-5,5) +elif model_case == 2: + param = (5, 0.5) # van der Pol + lim = (-7,7,-10,10) +else: + param = (0.02,0.5,0.2) # Fitzhugh-Nagumo + lim = (-7,7,-4,4) + +t, x_t = solve_flow(param,lim) + +if model_case == 1: + plt.savefig('Medio') +elif model_case == 2: + plt.savefig('VdP') +else: + plt.savefig('Fitzhugh-Nagumo') diff --git a/Flow3D.py b/Flow3D.py new file mode 100644 index 0000000..27d3662 --- /dev/null +++ b/Flow3D.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 16 07:38:57 2018 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +3D Flow examples: Lorenz, Rossler, Chua +""" +import numpy as np +import matplotlib as mpl +from mpl_toolkits.mplot3d import Axes3D +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') +fig = plt.figure() +ax = fig.add_axes([0, 0, 1, 1], projection='3d') +ax.axis('on') + +# model_case 1 = Lorenz +# model_case 2 = Rossler +# model_case 3 = Chua +model_case = int(input('Enter Model Case (1-3)')); + +def solve_lorenz(param, max_time=8.0, angle=0.0): + + if model_case == 1: +# Lorenz 3D flow + def flow_deriv(x_y_z, t0, sigma, beta, rho): + #"""Compute the time-derivative of a Lorenz system.""" + x, y, z = x_y_z + return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z] + model_title = 'Lorenz Attractor' + + elif model_case == 2: +# Rossler 3D flow + def flow_deriv(x_y_z, t0, sigma, beta, rho): + #"""Compute the time-derivative of a Medio system.""" + x, y, z = x_y_z + return [-y-z, x + sigma*y, beta + z*(x - rho)] + model_title = 'Rossler Attractor' + + else: +# Chua 3D flow + def flow_deriv(x_y_z, t0, alpha, beta, c, d): + #"""Compute the time-derivative of a Medio system.""" + x, y, z = x_y_z + f = c*x + 0.5*(d-c)*(abs(x+1)-abs(x-1)) + return [alpha*(y-x-f), x-y+z, -beta*y] + model_title = 'Chua Attractor' + + + + N=12 + colors = plt.cm.prism(np.linspace(0, 1, N)) + + + # Choose random starting points, uniformly distributed from -15 to 15 + np.random.seed(1) + x0 = init1 + init2*np.random.random((N, 3)) + + # Settle-down Solve for the trajectories + t = np.linspace(0, max_time/4, int(250*max_time/4)) + x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param) + for x0i in x0]) + + # Solve for trajectories + x0 = x_t[0:N,int(250*max_time/4)-1,0:3] + t = np.linspace(0, max_time, int(250*max_time)) + x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t, param) + for x0i in x0]) + + + # choose a different color for each trajectory + # colors = plt.cm.viridis(np.linspace(0, 1, N)) + # colors = plt.cm.rainbow(np.linspace(0, 1, N)) + # colors = plt.cm.spectral(np.linspace(0, 1, N)) + colors = plt.cm.prism(np.linspace(0, 1, N)) + + for i in range(N): + x, y, z = x_t[i,:,:].T + lines = ax.plot(x, y, z, '-', c=colors[i]) + plt.setp(lines, linewidth=0.5) + + ax.view_init(30, angle) + plt.show() + plt.title(model_title) + plt.savefig('Flow3D') + + return t, x_t + +if model_case == 1: + param = (10, 8/3, 28) # Lorenz + ax.set_xlim((-25, 25)) + ax.set_ylim((-35, 35)) + ax.set_zlim((5, 55)) + max_time = 50.0 + init1 = -15 + init2 = 30 +elif model_case == 2: + param = (0.2, 0.2, 5.7) # Rossler + ax.set_xlim((-15, 15)) + ax.set_ylim((-15, 15)) + ax.set_zlim((0, 20)) + max_time = 200 + init1 = -15 + init2 = 30 +else: + param = (15.6, 28.0, -0.7, -1.14 ) # Chua + ax.set_xlim((-3, 3)) + ax.set_ylim((-1, 1)) + ax.set_zlim((-3, 3)) + max_time = 100 + init1 = 0 + init2 = 0.1 + + +t, x_t = solve_lorenz(param, max_time,angle=30) + +plt.figure(2) +lines = plt.plot(t,x_t[1,:,0],t,x_t[1,:,1],t,x_t[1,:,2]) +plt.setp(lines, linewidth=1) + +for i in range(4): + plt.figure(3) + lines = plt.plot(x_t[i,:,0],x_t[i,:,1]) + plt.setp(lines,linewidth=0.5) + + plt.figure(4) + lines = plt.plot(x_t[i,:,1],x_t[i,:,2]) + plt.setp(lines,linewidth=0.5) + + plt.figure(5) + lines = plt.plot(x_t[i,:,0],x_t[i,:,2]) + plt.setp(lines,linewidth=0.5) + + diff --git a/Hamilton4D.py b/Hamilton4D.py new file mode 100644 index 0000000..20279bd --- /dev/null +++ b/Hamilton4D.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 18 06:03:32 2018 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +4D Hamiltonian flows +""" + +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 = Heiles +# model_case 2 = Crescent +print(' ') +print('Hamilton4D.py') +print('Case: 1 = Heiles') +print('Case: 2 = Crescent') +model_case = int(input('Enter the Model Case (1-2)')) + +if model_case == 1: + E = 1 # Heiles: 1, 0.3411 Crescent: 0.05, 1 + epsE = 0.3411 # 3411 + def flow_deriv(x_y_z_w,tspan): + x, y, z, w = x_y_z_w + a = z + b = w + c = -x - epsE*(2*x*y) + d = -y - epsE*(x**2 - y**2) + return[a,b,c,d] +else: + E = .095 # Crescent: 0.1, 1 + epsE = 1 + def flow_deriv(x_y_z_w,tspan): + x, y, z, w = x_y_z_w + a = z + b = w + c = -(epsE*(y-2*x**2)*(-4*x) + x) + d = -(y-epsE*2*x**2) + return[a,b,c,d] + +prms = np.sqrt(E) +pmax = np.sqrt(2*E) + +# Potential Function +if model_case == 1: + V = np.zeros(shape=(100,100)) + for xloop in range(100): + x = -2 + 4*xloop/100 + for yloop in range(100): + y = -2 + 4*yloop/100 + V[yloop,xloop] = 0.5*x**2 + 0.5*y**2 + epsE*(x**2*y - 0.33333*y**3) +else: + V = np.zeros(shape=(100,100)) + for xloop in range(100): + x = -2 + 4*xloop/100 + for yloop in range(100): + y = -2 + 4*yloop/100 + V[yloop,xloop] = 0.5*x**2 + 0.5*y**2 + epsE*(2*x**4 - 2*x**2*y) + +fig = plt.figure(1) +contr = plt.contourf(V,100, cmap=cm.coolwarm, vmin = 0, vmax = 10) +fig.colorbar(contr, shrink=0.5, aspect=5) +fig = plt.show() + +repnum = 250 +mulnum = 64/repnum + +np.random.seed(1) +for reploop in range(repnum): + px1 = 2*(np.random.random((1))-0.499)*pmax + py1 = np.sign(np.random.random((1))-0.499)*np.real(np.sqrt(2*(E-px1**2/2))) + xp1 = 0 + yp1 = 0 + + x_y_z_w0 = [xp1, yp1, px1, py1] + + tspan = np.linspace(1,1000,10000) + x_t = integrate.odeint(flow_deriv, x_y_z_w0, tspan) + siztmp = np.shape(x_t) + siz = siztmp[0] + + if reploop % 50 == 0: + plt.figure(2) + lines = plt.plot(x_t[:,0],x_t[:,1]) + plt.setp(lines, linewidth=0.5) + plt.show() + time.sleep(0.1) + #os.system("pause") + + y1 = x_t[:,0] + y2 = x_t[:,1] + y3 = x_t[:,2] + y4 = x_t[:,3] + + py = np.zeros(shape=(2*repnum,)) + yvar = np.zeros(shape=(2*repnum,)) + cnt = -1 + last = y1[1] + for loop in range(2,siz): + if (last < 0)and(y1[loop] > 0): + cnt = cnt+1 + del1 = -y1[loop-1]/(y1[loop] - y1[loop-1]) + py[cnt] = y4[loop-1] + del1*(y4[loop]-y4[loop-1]) + yvar[cnt] = y2[loop-1] + del1*(y2[loop]-y2[loop-1]) + last = y1[loop] + else: + last = y1[loop] + + plt.figure(3) + lines = plt.plot(yvar,py,'o',ms=1) + plt.show() + +if model_case == 1: + plt.savefig('Heiles') +else: + plt.savefig('Crescent') + + diff --git a/Heiles.py b/Heiles.py new file mode 100644 index 0000000..14cf838 --- /dev/null +++ b/Heiles.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Apr 18 06:03:32 2018 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Hamiltonian models +""" + +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') +print('Heiles.py') + +# model_case 1 = Heiles +# model_case 2 = Crescent +model_case = int(input('Enter the Model Case (1-2)')) + +if model_case == 1: + E = 1 # Heiles: 1, 0.3411 Crescent: 0.05, 1 + epsE = 0.3411 # 3411 + def flow_deriv(x_y_z_w,tspan): + x, y, z, w = x_y_z_w + a = z + b = w + c = -x - epsE*(2*x*y) + d = -y - epsE*(x**2 - y**2) + return[a,b,c,d] +else: + E = 1 # Heiles: 1, 0.3411 Crescent: 0.05, 1 + epsE = 0.2 # 0.2 + def flow_deriv(x_y_z_w,tspan): + x, y, z, w = x_y_z_w + a = z + b = w + c = -(epsE*(y-2*x**2)*(-4*x) + x) + d = -(y-epsE*2*x**2) + return[a,b,c,d] + +prms = np.sqrt(E) +pmax = np.sqrt(2*E) + + +# Potential Function +if model_case == 1: + V = np.zeros(shape=(100,100)) + for xloop in range(100): + x = -4 + 8*xloop/100 + for yloop in range(100): + y = -4 + 8*yloop/100 + V[yloop,xloop] = 0.5*x**2 + 0.5*y**2 + epsE*(x**2*y - 0.33333*y**3) +else: + V = np.zeros(shape=(100,100)) + for xloop in range(100): + x = -4 + 8*xloop/100 + for yloop in range(100): + y = -4 + 8*yloop/100 + V[yloop,xloop] = 0.5*x**2 + 0.5*y**2 + epsE*(2*x**4 - 2*x**2*y) + +mxV = np.int(np.max(V)) +fig = plt.figure(1) +contr = plt.contourf(V,mxV, cmap=cm.coolwarm, vmin = 0, vmax = 30) +fig.colorbar(contr, shrink=0.5, aspect=5) +fig = plt.show() + +repnum = 250 +mulnum = 64/repnum + +np.random.seed(1) +for reploop in range(repnum): + px1 = 2*(np.random.random((1))-0.499)*pmax + py1 = np.sign(np.random.random((1))-0.499)*np.real(np.sqrt(2*(E-px1**2/2))) + xp1 = 0 + yp1 = 0 + + x_y_z_w0 = [xp1, yp1, px1, py1] + + tspan = np.linspace(1,1000,10000) + x_t = integrate.odeint(flow_deriv, x_y_z_w0, tspan) + siztmp = np.shape(x_t) + siz = siztmp[0] + + if reploop % 50 == 0: + plt.figure(2) + lines = plt.plot(x_t[:,0],x_t[:,1]) + plt.setp(lines, linewidth=0.5) + plt.pause(0.05) + plt.show() + + y1 = x_t[:,0] + y2 = x_t[:,1] + y3 = x_t[:,2] + y4 = x_t[:,3] + + py = np.zeros(shape=(2*repnum,)) + yvar = np.zeros(shape=(2*repnum,)) + cnt = -1 + last = y1[1] + for loop in range(2,siz): + if (last < 0)and(y1[loop] > 0): # Test for x going positive + cnt = cnt+1 + del1 = -y1[loop-1]/(y1[loop] - y1[loop-1]) # Interpolate to the Poincare plane + py[cnt] = y4[loop-1] + del1*(y4[loop]-y4[loop-1]) + yvar[cnt] = y2[loop-1] + del1*(y2[loop]-y2[loop-1]) + last = y1[loop] + else: + last = y1[loop] + + plt.figure(3) + lines = plt.plot(yvar,py,'o',ms=1) + plt.pause(0.05) + plt.show() + +plt.savefig('Heiles') diff --git a/HenonHeiles.py b/HenonHeiles.py new file mode 100644 index 0000000..757859f --- /dev/null +++ b/HenonHeiles.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Jan 30 10:45:00 2019 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Henon-Heiles model of galactic dynamics +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') +print(' ') +print('HenonHeiles.py') + +E = 1 +epsE = 0.35 +pmax = np.sqrt(2*E) + +def flow_deriv(x_y_z_w, tspan): + x, y, z, w = x_y_z_w + a = z + b = w + c = -x - epsE*(2*x*y) + d = -y - epsE*(x**2 - y**2) + return [a,b,c,d] + +px1 = 2*(np.random.random(1)-0.499)*pmax; +py1 = np.sign(np.random.random(1)-0.5)*np.sqrt(2*(E-px1**2/2)); +xp1 = 0; +yp1 = 0; + +x_y_z_w = [xp1, yp1, px1, py1] + +tspan = np.linspace(1,1000,10000) +x_t = integrate.odeint(flow_deriv, x_y_z_w, tspan) + +y1 = x_t[:,0] +y2 = x_t[:,1] +y3 = x_t[:,2] +y4 = x_t[:,3] + +plt.figure(1) +lines = plt.plot(y1,y2) +plt.setp(lines,linewidth=0.5) + + +repnum = 256 # 32 +np.random.seed(1) +for reploop in range(2,repnum): + + px1 = 2*(np.random.random(1)-0.499)*pmax; + py1 = np.sign(np.random.random(1)-0.5)*np.sqrt(2*(E-px1**2/2)); + xp1 = 0; + yp1 = 0; + + x_y_z_w = [xp1, yp1, px1, py1] + + tspan = np.linspace(1,1000,10000) + x_t = integrate.odeint(flow_deriv, x_y_z_w, tspan) + siztmp = np.shape(x_t) + siz = siztmp[0] + + y1 = x_t[:,0] + y2 = x_t[:,1] + y3 = x_t[:,2] + y4 = x_t[:,3] + + py = np.zeros(shape=(2*repnum,)) + yvar = np.zeros(shape=(2*repnum,)) + cnt = -1 + last = y1[1] + for loop in range(2,siz): + if (last < 0)and(y1[loop] > 0): + cnt = cnt+1 + del1 = -y1[loop-1]/(y1[loop] - y1[loop-1]) + py[cnt] = y4[loop-1] + del1*(y4[loop]-y4[loop-1]) + yvar[cnt] = y2[loop-1] + del1*(y2[loop]-y2[loop-1]) + last = y1[loop] + else: + last = y1[loop] + + + plt.figure(2) + lines = plt.plot(yvar,py,'o',ms=1) + plt.show() + + + + + diff --git a/Hill.py b/Hill.py new file mode 100644 index 0000000..804e1e0 --- /dev/null +++ b/Hill.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue May 28 11:50:24 2019 +@author: nolte +Blog site: https://galileo-unbound.blog/2019/07/19/getting-armstrong-aldrin-and-collins-home-from-the-moon-apollo-11-and-the-three-body-problem/ +""" + +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') + +womega = 1 +R = 1 +eps = 1e-6 + +M1 = 1 +M2 = 1/10 +chsi = M2/M1 + +x1 = -M2*R/(M1+M2) +x2 = x1 + R + + +def poten(y,c): + + rp0 = np.sqrt(y**2 + c**2); + thetap0 = np.arctan(y/c); + + rp1 = np.sqrt(x1**2 + rp0**2 - 2*np.abs(rp0*x1)*np.cos(np.pi-thetap0)); + rp2 = np.sqrt(x2**2 + rp0**2 - 2*np.abs(rp0*x2)*np.cos(thetap0)); + V = -M1/rp1 -M2/rp2 - E; + + return [V] + +def flow_deriv(x_y_z,tspan): + x, y, z, w = x_y_z + + r1 = np.sqrt(x1**2 + x**2 - 2*np.abs(x*x1)*np.cos(np.pi-z)); + r2 = np.sqrt(x2**2 + x**2 - 2*np.abs(x*x2)*np.cos(z)); + + + yp = np.zeros(shape=(4,)) + yp[0] = y + yp[1] = -womega**2*R**3*(np.abs(x)-np.abs(x1)*np.cos(np.pi-z))/(r1**3+eps) - womega**2*R**3*chsi*(np.abs(x)-abs(x2)*np.cos(z))/(r2**3+eps) + x*(w-womega)**2 + yp[2] = w + yp[3] = 2*y*(womega-w)/x - womega**2*R**3*chsi*abs(x2)*np.sin(z)/(x*(r2**3+eps)) + womega**2*R**3*np.abs(x1)*np.sin(np.pi-z)/(x*(r1**3+eps)) + + return yp + + +r0 = 0.64 +v0 = 0.3 +theta0 = 0 +vrfrac = 1 + +rp1 = np.sqrt(x1**2 + r0**2 - 2*np.abs(r0*x1)*np.cos(np.pi-theta0)) +rp2 = np.sqrt(x2**2 + r0**2 - 2*np.abs(r0*x2)*np.cos(theta0)) +V = -M1/rp1 - M2/rp2 +T = 0.5*v0**2 +E = T + V + +vr = vrfrac*v0 +W = (2*T - v0**2)/r0 + + +y0 = [r0, vr, theta0, W] + +tspan = np.linspace(1,2000,20000) + +y = integrate.odeint(flow_deriv, y0, tspan) + +xx = y[1:20000,0]*np.cos(y[1:20000,2]); +yy = y[1:20000,0]*np.sin(y[1:20000,2]); + + +plt.figure(1) +lines = plt.plot(xx,yy) +plt.setp(lines, linewidth=0.5) +plt.show() diff --git a/Lozi.py b/Lozi.py new file mode 100644 index 0000000..9b3be73 --- /dev/null +++ b/Lozi.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed May 2 16:17:27 2018 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Lozi discrete map +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +#plt.close('all') + +B = -1 +C = 0.5 + +np.random.seed(2) + +plt.figure(1) + +for eloop in range(0,100): + + xlast = np.random.normal(0,1,1) + ylast = np.random.normal(0,1,1) + + xnew = np.zeros(shape=(500,)) + ynew = np.zeros(shape=(500,)) + for loop in range(0,500): + xnew[loop] = 1 + ylast - C*abs(xlast) + ynew[loop] = B*xlast + xlast = xnew[loop] + ylast = ynew[loop] + + plt.plot(np.real(xnew),np.real(ynew),'o',ms=1) + plt.xlim(xmin=-1.25,xmax=2) + plt.ylim(ymin=-2,ymax=1.25) + + +plt.savefig('Lozi') + + diff --git a/NetDynamics.py b/NetDynamics.py new file mode 100644 index 0000000..d4d3153 --- /dev/null +++ b/NetDynamics.py @@ -0,0 +1,238 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 11 08:56:41 2019 + +@author: nolte + +D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019) +""" + +# https://www.python-course.eu/networkx.php +# https://networkx.github.io/documentation/stable/tutorial.html +# https://networkx.github.io/documentation/stable/reference/functions.html + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt +import networkx as nx +from UserFunction import linfit +import time + +tstart = time.time() + +plt.close('all') + + + +Nfac = 25 # 25 +N = 36 # 50 +width = 0.2 + +# model_case 1 = Complete Graph +# model_case 2 = Barabasi +# model_case 3 = Power Law +# model_case 4 = D-dimensional Hypercube +# model_case 5 = Erdos Renyi +# model_case 6 = Random Regular +# model_case 7 = Strogatz +# model_case 8 = Hexagonal lattice +# model_case 9 = Tree +# model_case 10 == 2D square lattice + +model_case = int(input('Input Model Case (1-10)')) +if model_case == 1: # Complete Graph + facoef = 0.2 + nodecouple = nx.complete_graph(N) +elif model_case == 2: # Barabasi + facoef = 2 + k = 3 + nodecouple = nx.barabasi_albert_graph(N, k, seed=None) +elif model_case == 3: # Power law + facoef = 3 + k = 3 + triangle_prob = 0.3 + nodecouple = nx.powerlaw_cluster_graph(N, k, triangle_prob) +elif model_case == 4: + Dim = 6 + facoef = 3 + nodecouple = nx.hypercube_graph(Dim) + N = nodecouple.number_of_nodes() +elif model_case == 5: # Erdos-Renyi + facoef = 5 + prob = 0.1 + nodecouple = nx.erdos_renyi_graph(N, prob, seed=None, directed=False) +elif model_case == 6: # Random + facoef = 5 + nodecouple = nx.random_regular_graph(3, N, seed=None) +elif model_case == 7: # Watts + facoef = 7 + k = 4; # nearest neighbors + rewire_prob = 0.2 # rewiring probability + nodecouple = nx.watts_strogatz_graph(N, k, rewire_prob, seed=None) +elif model_case == 8: + facoef = 8 + rows = 4 + colm = 8 + nodecouple = nx.hexagonal_lattice_graph(rows, colm, periodic=True, with_positions=False) + N = nodecouple.number_of_nodes() +elif model_case == 9: # k-fold tree + facoef = 16 + k = 3 # degree + h = 3 # height + sm = 0 + for lp in range(h+1): + sm = sm + k**lp + N = sm + nodecouple = nx.balanced_tree(k, h) +elif model_case == 10: # square lattice + facoef = 3 + m = 6 + n = 6 + nodecouple = nx.grid_2d_graph(m, n, periodic=False) + N = nodecouple.number_of_nodes() + +plt.figure(1) +nx.draw(nodecouple) +#nx.draw_circular(nodecouple) +#nx.draw_spring(nodecouple) +#nx.draw_spectral(nodecouple) +print('Number of nodes = ',nodecouple.number_of_nodes()) +print('Number of edges = ',nodecouple.number_of_edges()) +#print('Average degree = ',nx.k_nearest_neighbors(nodecouple)) + +# function: omegout, yout = coupleN(G) +def coupleN(G): + + # function: yd = flow_deriv(x_y) + def flow_deriv(y,t0): + + yp = np.zeros(shape=(N,)) + ind = -1 + for omloop in G.node: + ind = ind + 1 + temp = omega[ind] + linksz = G.node[omloop]['numlink'] + for cloop in range(linksz): + cindex = G.node[omloop]['link'][cloop] + indx = G.node[cindex]['index'] + g = G.node[omloop]['coupling'][cloop] + + + temp = temp + g*np.sin(y[indx]-y[ind]) + + yp[ind] = temp + + yd = np.zeros(shape=(N,)) + for omloop in range(N): + yd[omloop] = yp[omloop] + + return yd + # end of function flow_deriv(x_y) + + mnomega = 1.0 + + ind = -1 + for nodeloop in G.node: + ind = ind + 1 + omega[ind] = G.node[nodeloop]['element'] + + x_y_z = omega + + # Settle-down Solve for the trajectories + tsettle = 100 + t = np.linspace(0, tsettle, tsettle) + x_t = integrate.odeint(flow_deriv, x_y_z, t) + x0 = x_t[tsettle-1,0:N] + + t = np.linspace(0,2000,2000) + y = integrate.odeint(flow_deriv, x0, t) + siztmp = np.shape(y) + sy = siztmp[0] + + # Fit the frequency + m = np.zeros(shape = (N,)) + w = np.zeros(shape = (N,)) + mtmp = np.zeros(shape=(4,)) + btmp = np.zeros(shape=(4,)) + for omloop in range(N): + + if np.remainder(sy,4) == 0: + mtmp[0],btmp[0] = linfit(t[0:sy//2],y[0:sy//2,omloop]); + mtmp[1],btmp[1] = linfit(t[sy//2+1:sy],y[sy//2+1:sy,omloop]); + mtmp[2],btmp[2] = linfit(t[sy//4+1:3*sy//4],y[sy//4+1:3*sy//4,omloop]); + mtmp[3],btmp[3] = linfit(t,y[:,omloop]); + else: + sytmp = 4*np.floor(sy/4); + mtmp[0],btmp[0] = linfit(t[0:sytmp//2],y[0:sytmp//2,omloop]); + mtmp[1],btmp[1] = linfit(t[sytmp//2+1:sytmp],y[sytmp//2+1:sytmp,omloop]); + mtmp[2],btmp[2] = linfit(t[sytmp//4+1:3*sytmp/4],y[sytmp//4+1:3*sytmp//4,omloop]); + mtmp[3],btmp[3] = linfit(t[0:sytmp],y[0:sytmp,omloop]); + + + #m[omloop] = np.median(mtmp) + m[omloop] = np.mean(mtmp) + + w[omloop] = mnomega + m[omloop] + + omegout = m + yout = y + + return omegout, yout + # end of function: omegout, yout = coupleN(G) + + + +Nlink = N*(N-1)//2 +omega = np.zeros(shape=(N,)) +omegatemp = width*(np.random.rand(N)-1) +meanomega = np.mean(omegatemp) +omega = omegatemp - meanomega +sto = np.std(omega) + + + +lnk = np.zeros(shape = (N,), dtype=int) +ind = -1 +for loop in nodecouple.node: + ind = ind + 1 + nodecouple.node[loop]['index'] = ind + nodecouple.node[loop]['element'] = omega[ind] + nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop)) + nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop))) + lnk[ind] = len(list(nx.neighbors(nodecouple,loop))) + + +avgdegree = np.mean(lnk) +mnomega = 1 + +facval = np.zeros(shape = (Nfac,)) +yy = np.zeros(shape=(Nfac,N)) +xx = np.zeros(shape=(Nfac,)) +for facloop in range(Nfac): + print(facloop) + + fac = facoef*(16*facloop/(Nfac))*(1/(N-1))*sto/mnomega + ind = -1 + for nodeloop in nodecouple.node: + ind = ind + 1 + nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],)) + for linkloop in range (lnk[ind]): + nodecouple.node[nodeloop]['coupling'][linkloop] = fac + + facval[facloop] = fac*avgdegree + + omegout, yout = coupleN(nodecouple) # Here is the subfunction call for the flow + + for omloop in range(N): + yy[facloop,omloop] = omegout[omloop] + + xx[facloop] = facval[facloop] + +plt.figure(2) +lines = plt.plot(xx,yy) +plt.setp(lines, linewidth=0.5) +plt.show() + +elapsed_time = time.time() - tstart +print('elapsed time = ',format(elapsed_time,'.2f'),'secs') diff --git a/NetSIR.py b/NetSIR.py new file mode 100644 index 0000000..babfac8 --- /dev/null +++ b/NetSIR.py @@ -0,0 +1,293 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 11 08:56:41 2019 + +@author: nolte + +D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019) +""" + +# https://www.python-course.eu/networkx.php +# https://networkx.github.io/documentation/stable/tutorial.html +# https://networkx.github.io/documentation/stable/reference/functions.html + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt +import networkx as nx +import time +from random import random + +tstart = time.time() + +plt.close('all') + +r = 0.0001 # 0.00015 +k = 10 # connections 50 +dill = 14 # days ill +dpq = 14 # days shelter in place +fnq = 0.5 # fraction NOT sheltering in place +mr0 = 0.01 # mortality rate +mr1 = 0.03 # extra mortality rate if exceeding hospital capacity +P = 330 # population of US in Millions +HC = 0.003 # hospital capacity + +dinf = fnq*dill + (1-fnq)*np.exp(-dpq/dill)*dill; + +#betap = r*k*dinf; +#mu = 1/dill; +betap = 0.014; +mu = 0.13; + +print('beta = ',betap) +print('dinf = ',dinf) +print('betap/mu = ',betap/mu) + + +N = 128 # 50 + +# model_case 1 = Complete Graph +# model_case 2 = Barabasi +# model_case 3 = Power Law +# model_case 4 = D-dimensional Hypercube +# model_case 5 = Erdos Renyi +# model_case 6 = Random Regular +# model_case 7 = Strogatz +# model_case 8 = Hexagonal lattice +# model_case 9 = Tree + +model_case = int(input('Input Model Case (1-9)')) +if model_case == 1: # Complete Graph + facoef = 0.5 + nodecouple = nx.complete_graph(N) +elif model_case == 2: # Barabasi + facoef = 2 + k = 1 + nodecouple = nx.barabasi_albert_graph(N, k, seed=None) +elif model_case == 3: # Power law + facoef = 3 + k = 2 + triangle_prob = 0.05 + nodecouple = nx.powerlaw_cluster_graph(N, k, triangle_prob) +elif model_case == 4: + Dim = 6 + facoef = 3 + nodecouple = nx.hypercube_graph(Dim) + N = nodecouple.number_of_nodes() +elif model_case == 5: # Erdos-Renyi + facoef = 5 + prob = 0.03 + nodecouple = nx.erdos_renyi_graph(N, prob, seed=None, directed=False) +elif model_case == 6: # Random + facoef = 10 + nodecouple = nx.random_regular_graph(3, N, seed=None) +elif model_case == 7: # Watts + facoef = 7 + k = 4; # nearest neighbors + rewire_prob = 0.2 # rewiring probability + nodecouple = nx.watts_strogatz_graph(N, k, rewire_prob, seed=None) +elif model_case == 8: # Hexagonal lattice + facoef = 8 + rows = 4 + colm = 8 + nodecouple = nx.hexagonal_lattice_graph(rows, colm, periodic=True, with_positions=False) + N = nodecouple.number_of_nodes() +elif model_case == 9: # k-fold tree + facoef = 16 + k = 3 # degree + h = 4 # height + sm = 0 + for lp in range(h+1): + sm = sm + k**lp + N = sm + nodecouple = nx.balanced_tree(k, h) + +indhi = 0 +deg = 0 +for omloop in nodecouple.node: + degtmp = nodecouple.degree(omloop) + if degtmp > deg: + deg = degtmp + indhi = omloop +print('highest degree node = ',indhi) +print('highest degree = ',deg) + +plt.figure(1) +colors = [(random(), random(), random()) for _i in range(10)] +#nx.draw(nodecouple) +nx.draw_circular(nodecouple,node_size=75, node_color=colors) +#nx.draw_spring(nodecouple,node_size=75, node_color=colors) +#nx.draw_spectral(nodecouple) + +#print('Number of nodes = ',nodecouple.number_of_nodes()) +#print('Number of edges = ',nodecouple.number_of_edges()) +#print('Average degree = ',nx.k_nearest_neighbors(nodecouple)) +print(nx.info(nodecouple)) + +# function: omegout, yout = coupleN(G) +def coupleN(G,tlim): + + # function: yd = flow_deriv(x_y) + def flow_deriv(x_y,t0): + + N = np.int(x_y.size/2) + yd = np.zeros(shape=(2*N,)) + ind = -1 + for omloop in G.node: + ind = ind + 1 + temp1 = -mu*x_y[ind] + betap*x_y[ind]*x_y[N+ind] + temp2 = -betap*x_y[ind]*x_y[N+ind] + linksz = G.node[omloop]['numlink'] + for cloop in range(linksz): + cindex = G.node[omloop]['link'][cloop] + indx = G.node[cindex]['index'] + g = G.node[omloop]['coupling'][cloop] + + temp1 = temp1 + g*betap*x_y[indx]*x_y[N+ind] + temp2 = temp2 - g*betap*x_y[indx]*x_y[N+ind] + + yd[ind] = temp1 + yd[N+ind] = temp2 + + return yd + # end of function flow_deriv(x_y) + x0 = x_y + t = np.linspace(0,tlim,tlim) # 600 300 + y = integrate.odeint(flow_deriv, x0, t) + + return t,y + # end of function: omegout, yout = coupleN(G) + +lnk = np.zeros(shape = (N,), dtype=int) +ind = -1 +for loop in nodecouple.node: + ind = ind + 1 + nodecouple.node[loop]['index'] = ind + nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop)) + nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop))) + lnk[ind] = len(list(nx.neighbors(nodecouple,loop))) + +gfac = 0.1 + +ind = -1 +for nodeloop in nodecouple.node: + ind = ind + 1 + nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],)) + for linkloop in range (lnk[ind]): + nodecouple.node[nodeloop]['coupling'][linkloop] = gfac*facoef + +x_y = np.zeros(shape=(2*N,)) +for loop in nodecouple.node: + x_y[loop]=0 + x_y[N+loop]=nodecouple.degree(loop) + #x_y[N+loop]=1 +x_y[N-1]=1 +x_y[2*N-1] = x_y[2*N-1] - 0.1 +N0 = np.sum(x_y[N:2*N]) - x_y[indhi] - x_y[N+indhi] + +tlim0 = 600 +t0,yout0 = coupleN(nodecouple,tlim0) # Here is the subfunction call for the flow + + +plt.figure(2) +plt.yscale('log') +plt.gca().set_ylim(1e-3, 1) +for loop in range(N): + lines1 = plt.plot(t0,yout0[:,loop]) + lines2 = plt.plot(t0,yout0[:,N+loop]) + lines3 = plt.plot(t0,N0-yout0[:,loop]-yout0[:,N+loop]) + + plt.setp(lines1, linewidth=0.5) + plt.setp(lines2, linewidth=0.5) + plt.setp(lines3, linewidth=0.5) + + +Itot = np.sum(yout0[:,0:127],axis = 1) - yout0[:,indhi] +Stot = np.sum(yout0[:,128:255],axis = 1) - yout0[:,N+indhi] +Rtot = N0 - Itot - Stot +plt.figure(3) +plt.plot(t0,Itot,t0,Stot,t0,Rtot) +plt.legend(('Infected','Susceptible','Removed')) +plt.hold + + + + +# Repeat but innoculate +x_y = np.zeros(shape=(2*N,)) +for loop in nodecouple.node: + x_y[loop]=0 + x_y[N+loop]=nodecouple.degree(loop) + #x_y[N+loop]=1 +x_y[N-1]=1 +x_y[2*N-1] = x_y[2*N-1] - 0.1 +N0 = np.sum(x_y[N:2*N]) - x_y[indhi] - x_y[N+indhi] + +tlim0 = 50 +t0,yout0 = coupleN(nodecouple,tlim0) # Here is the subfunction call for the flow + + +# remove all edges from highest-degree node +ee = list(nodecouple.edges(indhi)) +nodecouple.remove_edges_from(ee) +print(nx.info(nodecouple)) + +#nodecouple.remove_node(indhi) +lnk = np.zeros(shape = (N,), dtype=int) +ind = -1 +for loop in nodecouple.node: + ind = ind + 1 + nodecouple.node[loop]['index'] = ind + nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop)) + nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop))) + lnk[ind] = len(list(nx.neighbors(nodecouple,loop))) + +ind = -1 +x_y = np.zeros(shape=(2*N,)) +for nodeloop in nodecouple.node: + ind = ind + 1 + nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],)) + x_y[ind] = yout0[tlim0-1,nodeloop] + x_y[N+ind] = yout0[tlim0-1,N+nodeloop] + for linkloop in range (lnk[ind]): + nodecouple.node[nodeloop]['coupling'][linkloop] = gfac*facoef + +#N1 = np.sum(yout0[tlim0-1,:]) +#x_y = yout0[tlim0-1,:] + +tlim1 = 500 +t1,yout1 = coupleN(nodecouple,tlim1) + +t = np.zeros(shape=(tlim0+tlim1,)) +yout = np.zeros(shape=(tlim0+tlim1,2*N)) +t[0:tlim0] = t0 +t[tlim0:tlim1+tlim0] = tlim0+t1 +yout[0:tlim0,:] = yout0 +yout[tlim0:tlim1+tlim0,:] = yout1 + + +plt.figure(4) +plt.yscale('log') +plt.gca().set_ylim(1e-3, 1) +for loop in range(N): + lines1 = plt.plot(t,yout[:,loop]) + lines2 = plt.plot(t,yout[:,N+loop]) + lines3 = plt.plot(t,N0-yout[:,loop]-yout[:,N+loop]) + + plt.setp(lines1, linewidth=0.5) + plt.setp(lines2, linewidth=0.5) + plt.setp(lines3, linewidth=0.5) + + +Itot = np.sum(yout[:,0:127],axis = 1) - yout[:,indhi] +Stot = np.sum(yout[:,128:255],axis = 1) - yout[:,N+indhi] +Rtot = N0 - Itot - Stot +plt.figure(3) +plt.plot(t,Itot,t,Stot,t,Rtot,linestyle='dashed') +plt.legend(('Infected','Susceptible','Removed')) +plt.hold() + + +elapsed_time = time.time() - tstart +print('elapsed time = ',format(elapsed_time,'.2f'),'secs') diff --git a/NetSIRSF.py b/NetSIRSF.py new file mode 100644 index 0000000..8e7a9fa --- /dev/null +++ b/NetSIRSF.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat May 11 08:56:41 2019 + +@author: nolte + +D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019) +""" + +# https://www.python-course.eu/networkx.php +# https://networkx.github.io/documentation/stable/tutorial.html +# https://networkx.github.io/documentation/stable/reference/functions.html + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt +import networkx as nx +import time +from random import random + +tstart = time.time() + +plt.close('all') + +betap = 0.014; +mu = 0.13; + +print('beta = ',betap) +print('betap/mu = ',betap/mu) + + +N = 128 # 50 + + +facoef = 2 +k = 1 +nodecouple = nx.barabasi_albert_graph(N, k, seed=None) + +indhi = 0 +deg = 0 +for omloop in nodecouple.node: + degtmp = nodecouple.degree(omloop) + if degtmp > deg: + deg = degtmp + indhi = omloop +print('highest degree node = ',indhi) +print('highest degree = ',deg) + +plt.figure(1) +colors = [(random(), random(), random()) for _i in range(10)] +nx.draw_circular(nodecouple,node_size=75, node_color=colors) +print(nx.info(nodecouple)) + +# function: omegout, yout = coupleN(G) +def coupleN(G,tlim): + + # function: yd = flow_deriv(x_y) + def flow_deriv(x_y,t0): + + N = np.int(x_y.size/2) + yd = np.zeros(shape=(2*N,)) + ind = -1 + for omloop in G.node: + ind = ind + 1 + temp1 = -mu*x_y[ind] + betap*x_y[ind]*x_y[N+ind] + temp2 = -betap*x_y[ind]*x_y[N+ind] + linksz = G.node[omloop]['numlink'] + for cloop in range(linksz): + cindex = G.node[omloop]['link'][cloop] + indx = G.node[cindex]['index'] + g = G.node[omloop]['coupling'][cloop] + + temp1 = temp1 + g*betap*x_y[indx]*x_y[N+ind] + temp2 = temp2 - g*betap*x_y[indx]*x_y[N+ind] + + yd[ind] = temp1 + yd[N+ind] = temp2 + + return yd + # end of function flow_deriv(x_y) + x0 = x_y + t = np.linspace(0,tlim,tlim) # 600 300 + y = integrate.odeint(flow_deriv, x0, t) + + return t,y + # end of function: omegout, yout = coupleN(G) + +lnk = np.zeros(shape = (N,), dtype=int) +ind = -1 +for loop in nodecouple.node: + ind = ind + 1 + nodecouple.node[loop]['index'] = ind + nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop)) + nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop))) + lnk[ind] = len(list(nx.neighbors(nodecouple,loop))) + +gfac = 0.1 + +ind = -1 +for nodeloop in nodecouple.node: + ind = ind + 1 + nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],)) + for linkloop in range (lnk[ind]): + nodecouple.node[nodeloop]['coupling'][linkloop] = gfac*facoef + +x_y = np.zeros(shape=(2*N,)) +for loop in nodecouple.node: + x_y[loop]=0 + x_y[N+loop]=nodecouple.degree(loop) + #x_y[N+loop]=1 +x_y[N-1 ]= 0.01 +x_y[2*N-1] = x_y[2*N-1] - 0.01 +N0 = np.sum(x_y[N:2*N]) - x_y[indhi] - x_y[N+indhi] +print('N0 = ',N0) + +tlim0 = 600 +t0,yout0 = coupleN(nodecouple,tlim0) # Here is the subfunction call for the flow + + +plt.figure(2) +plt.yscale('log') +plt.gca().set_ylim(1e-3, 1) +for loop in range(N): + lines1 = plt.plot(t0,yout0[:,loop]) + lines2 = plt.plot(t0,yout0[:,N+loop]) + lines3 = plt.plot(t0,N0-yout0[:,loop]-yout0[:,N+loop]) + + plt.setp(lines1, linewidth=0.5) + plt.setp(lines2, linewidth=0.5) + plt.setp(lines3, linewidth=0.5) + + +Itot = np.sum(yout0[:,0:127],axis = 1) - yout0[:,indhi] +Stot = np.sum(yout0[:,128:255],axis = 1) - yout0[:,N+indhi] +Rtot = N0 - Itot - Stot +plt.figure(3) +#plt.plot(t0,Itot,'r',t0,Stot,'g',t0,Rtot,'b') +plt.plot(t0,Itot/N0,'r',t0,Rtot/N0,'b') +#plt.legend(('Infected','Susceptible','Removed')) +plt.legend(('Infected','Removed')) +plt.hold + +# Repeat but innoculate highest-degree node +x_y = np.zeros(shape=(2*N,)) +for loop in nodecouple.node: + x_y[loop]=0 + x_y[N+loop]=nodecouple.degree(loop) + #x_y[N+loop]=1 +x_y[N-1] = 0.01 +x_y[2*N-1] = x_y[2*N-1] - 0.01 +N0 = np.sum(x_y[N:2*N]) - x_y[indhi] - x_y[N+indhi] + +tlim0 = 50 +t0,yout0 = coupleN(nodecouple,tlim0) + + +# remove all edges from highest-degree node +ee = list(nodecouple.edges(indhi)) +nodecouple.remove_edges_from(ee) +print(nx.info(nodecouple)) + +#nodecouple.remove_node(indhi) +lnk = np.zeros(shape = (N,), dtype=int) +ind = -1 +for loop in nodecouple.node: + ind = ind + 1 + nodecouple.node[loop]['index'] = ind + nodecouple.node[loop]['link'] = list(nx.neighbors(nodecouple,loop)) + nodecouple.node[loop]['numlink'] = len(list(nx.neighbors(nodecouple,loop))) + lnk[ind] = len(list(nx.neighbors(nodecouple,loop))) + +ind = -1 +x_y = np.zeros(shape=(2*N,)) +for nodeloop in nodecouple.node: + ind = ind + 1 + nodecouple.node[nodeloop]['coupling'] = np.zeros(shape=(lnk[ind],)) + x_y[ind] = yout0[tlim0-1,nodeloop] + x_y[N+ind] = yout0[tlim0-1,N+nodeloop] + for linkloop in range (lnk[ind]): + nodecouple.node[nodeloop]['coupling'][linkloop] = gfac*facoef + + +tlim1 = 500 +t1,yout1 = coupleN(nodecouple,tlim1) + +t = np.zeros(shape=(tlim0+tlim1,)) +yout = np.zeros(shape=(tlim0+tlim1,2*N)) +t[0:tlim0] = t0 +t[tlim0:tlim1+tlim0] = tlim0+t1 +yout[0:tlim0,:] = yout0 +yout[tlim0:tlim1+tlim0,:] = yout1 + + +plt.figure(4) +plt.yscale('log') +plt.gca().set_ylim(1e-3, 1) +for loop in range(N): + lines1 = plt.plot(t,yout[:,loop]) + lines2 = plt.plot(t,yout[:,N+loop]) + lines3 = plt.plot(t,N0-yout[:,loop]-yout[:,N+loop]) + + plt.setp(lines1, linewidth=0.5) + plt.setp(lines2, linewidth=0.5) + plt.setp(lines3, linewidth=0.5) + + +Itot = np.sum(yout[:,0:127],axis = 1) - yout[:,indhi] +Stot = np.sum(yout[:,128:255],axis = 1) - yout[:,N+indhi] +Rtot = N0 - Itot - Stot +plt.figure(3) +#plt.plot(t,Itot,'r',t,Stot,'g',t,Rtot,'b',linestyle='dashed') +plt.plot(t,Itot/N0,'r',t,Rtot/N0,'b',linestyle='dashed') +#plt.legend(('Infected','Susceptible','Removed')) +plt.legend(('Infected','Removed')) +plt.xlabel('Days') +plt.ylabel('Fraction of Sub-Population') +plt.title('Network Dynamics for COVID-19') +plt.show() +plt.hold() + +elapsed_time = time.time() - tstart +print('elapsed time = ',format(elapsed_time,'.2f'),'secs') diff --git a/PenInverted.py b/PenInverted.py new file mode 100644 index 0000000..77ae592 --- /dev/null +++ b/PenInverted.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed May 21 06:03:32 2018 + +@author: nolte +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +# model_case 1 = Pendulum +# model_case 2 = Double Well +print(' ') +print('PenInverted.py') + +F = 133.5 # 30 to 140 (133.5) +delt = 0.000 # 0.000 to 0.01 +w = 20 # 10 +def flow_deriv(x_y_z,tspan): + x, y, z = x_y_z + a = y + b = -(1 + F*np.cos(z))*np.sin(x) - delt*y + c = w + return[a,b,c] + + +T = 2*np.pi/w + +x0 = np.pi+0.3 +v0 = 0.00 +z0 = 0 + +x_y_z = [x0, v0, z0] + +# Solve for the trajectories +t = np.linspace(0, 2000, 200000) +x_t = integrate.odeint(flow_deriv, x_y_z, t) +siztmp = np.shape(x_t) +siz = siztmp[0] + +#y1 = np.mod(x_t[:,0]-np.pi,2*np.pi)-np.pi +y1 = x_t[:,0] +y2 = x_t[:,1] +y3 = x_t[:,2] + +plt.figure(1) +lines = plt.plot(t[0:2000],x_t[0:2000,0]/np.pi) +plt.setp(lines, linewidth=0.5) +plt.show() +plt.title('Angular Position') + +plt.figure(2) +lines = plt.plot(t[0:1000],y2[0:1000]) +plt.setp(lines, linewidth=0.5) +plt.show() +plt.title('Speed') + +repnum = 5000 +px = np.zeros(shape=(2*repnum,)) +xvar = np.zeros(shape=(2*repnum,)) +cnt = -1 +testwt = np.mod(t,T)-0.5*T; +last = testwt[1] +for loop in range(2,siz-1): + 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) +lines = plt.plot(xvar[0:5000],px[0:5000],'ko',ms=1) +plt.show() +plt.title('First Return Map') + +plt.figure(4) +lines = plt.plot(x_t[0:1000,0]/np.pi,y2[0:1000]) +plt.setp(lines, linewidth=0.5) +plt.show() +plt.title('Phase Space') diff --git a/Perturbed.py b/Perturbed.py new file mode 100644 index 0000000..fd78709 --- /dev/null +++ b/Perturbed.py @@ -0,0 +1,172 @@ +#!/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') + + diff --git a/SIR.py b/SIR.py new file mode 100644 index 0000000..a1f105b --- /dev/null +++ b/SIR.py @@ -0,0 +1,148 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat March 21 2020 + +@author: nolte +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +print(' ') +print('SIR.py') + +def solve_flow(param,max_time=1000.0): + + def flow_deriv(x_y,tspan,mu,betap): + x, y = x_y + + return [-mu*x + betap*x*y,-betap*x*y] + + x0 = [del1, del2] + + # Solve for the trajectories + t = np.linspace(0, int(tlim), int(250*tlim)) + x_t = integrate.odeint(flow_deriv, x0, t, param) + + + return t, x_t + + +r = 0.0002 # 0.00015 +k = 50 # connections 50 +dill = 14 # days ill +dpq = 14 # days shelter in place +fnq = 0.6 # fraction NOT sheltering in place +mr0 = 0.01 # mortality rate +mr1 = 0.03 # extra mortality rate if exceeding hospital capacity +P = 330 # population of US in Millions +HC = 0.003 # hospital capacity + +dinf = fnq*dill + (1-fnq)*np.exp(-dpq/dill)*dill; + +betap = r*k*dinf; +mu = 1/dill; + +print('beta = ',betap) +print('dinf = ',dinf) +print('betap/mu = ',betap/mu) + +del1 = .001 # infected +del2 = 1-del1 # susceptible + +tlim = np.log(P*1e6/del1)/betap + 50/betap + +param = (mu, betap) # flow parameters + +t, y = solve_flow(param) +I = y[:,0] +S = y[:,1] +R = 1 - I - S + +plt.figure(1) +lines = plt.semilogy(t,I,t,S,t,R) +plt.ylim([0.001,1]) +plt.xlim([0,tlim]) +plt.legend(('Infected','Susceptible','Removed')) +plt.setp(lines, linewidth=0.5) +plt.xlabel('Days') +plt.ylabel('Fraction of Population') +plt.title('Population Dynamics for COVID-19 in US') +plt.show() + +mr = mr0 + mr1*(0.2*np.max(I)-HC)*np.heaviside(0.2*np.max(I),HC) +Dead = mr*P*R[R.size-1] +print('Dead = ',Dead) + + +D = np.zeros(shape=(100,)) +x = np.zeros(shape=(100,)) +for kloop in range(0,5): + for floop in range(0,100): + + fnq = floop/100 + + dinf = fnq*dill + (1-fnq)*np.exp(-dpq/dill)*dill; + + k = 20 + kloop*10 + betap = r*k*dinf + + tlim = np.log(P*1e6/del1)/betap + 50/betap + + param = (mu, betap) # flow parameters + + t, y = solve_flow(param) + I = y[:,0] + S = y[:,1] + R = 1 - I - S + + mr = mr0 + mr1*(0.2*np.max(I)-HC)*np.heaviside(0.2*np.max(I),HC) + + D[floop] = mr*P*R[R.size-1] + x[floop] = fnq + + plt.figure(2) + lines2 = plt.plot(x,D) + plt.setp(lines2, linewidth=0.5) + +plt.ylabel('US Million Deaths') +plt.xlabel('Fraction NOT Quarantining') +plt.title('Quarantine and Distancing') +plt.legend(('20','30','40','50','60','70')) +plt.show() + + +label = np.zeros(shape=(9,)) +for floop in range(0,8): + + fq = floop/10.0 + + dinf = (1-fq)*dill + fq*np.exp(-dpq/dill)*dill; + + k = 50 + betap = r*k*dinf + + tlim = np.log(P*1e6/del1)/betap + 50/betap + + param = (mu, betap) # flow parameters + + t, y = solve_flow(param) + I = y[:,0] + S = y[:,1] + R = 1 - I - S + + plt.figure(3) + lines2 = plt.plot(t,I*P) + plt.setp(lines2, linewidth=0.5) + label[floop]=fq + +plt.legend(label) +plt.ylabel('US Millions Infected') +plt.xlabel('Days') +plt.title('Flattening the Curve') + + diff --git a/SIRS.py b/SIRS.py new file mode 100644 index 0000000..e0f6146 --- /dev/null +++ b/SIRS.py @@ -0,0 +1,171 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Fri July 17 2020 + +D. D. Nolte, "Introduction to Modern Dynamics: + Chaos, Networks, Space and Time, 2nd Edition (Oxford University Press, 2019) + +@author: nolte +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + + +def tripartite(x,y,z): + + sm = x + y + z + xp = x/sm + yp = y/sm + + f = np.sqrt(3)/2 + + y0 = f*xp + x0 = -0.5*xp - yp + 1; + + lines = plt.plot(x0,y0) + plt.setp(lines, linewidth=0.5) + plt.plot([0, 1],[0, 0],'k',linewidth=1) + plt.plot([0, 0.5],[0, f],'k',linewidth=1) + plt.plot([1, 0.5],[0, f],'k',linewidth=1) + plt.show() + + + +print(' ') +print('SIRS.py') + +def solve_flow(param,max_time=1000.0): + + def flow_deriv(x_y,tspan,mu,betap,nu): + x, y = x_y + + return [-mu*x + betap*x*(1-x-y),mu*x-nu*y] + + x0 = [del1, del2] + + # Solve for the trajectories + t = np.linspace(0, int(tlim), int(250*tlim)) + x_t = integrate.odeint(flow_deriv, x0, t, param) + + + return t, x_t + + + # rates per week +betap = 0.3; # infection rate +mu = 0.2; # recovery rate +nu = 0.02 # immunity decay rate + +print('beta = ',betap) +print('mu = ',mu) +print('nu =',nu) +print('betap/mu = ',betap/mu) + +del1 = 0.005 # initial infected +del2 = 0.005 # recovered + +tlim = 600 # weeks (about 12 years) + +param = (mu, betap, nu) # flow parameters + +t, y = solve_flow(param) +I = y[:,0] +R = y[:,1] +S = 1 - I - R + +plt.figure(1) +lines = plt.semilogy(t,I,t,S,t,R) +plt.ylim([0.001,1]) +plt.xlim([0,tlim]) +plt.legend(('Infected','Susceptible','Recovered')) +plt.setp(lines, linewidth=0.5) +plt.xlabel('Days') +plt.ylabel('Fraction of Population') +plt.title('Population Dynamics for COVID-19') +plt.show() + +plt.figure(2) +plt.hold(True) +for xloop in range(0,10): + del1 = xloop/10.1 + 0.001 + del2 = 0.01 + + tlim = 300 + param = (mu, betap, nu) # flow parameters + t, y = solve_flow(param) + I = y[:,0] + R = y[:,1] + S = 1 - I - R + + tripartite(I,R,S); + +for yloop in range(1,6): + del1 = 0.001; + del2 = yloop/10.1 + t, y = solve_flow(param) + I = y[:,0] + R = y[:,1] + S = 1 - I - R + + tripartite(I,R,S); + +for loop in range(2,10): + del1 = loop/10.1 + del2 = 1 - del1 - 0.01 + t, y = solve_flow(param) + I = y[:,0] + R = y[:,1] + S = 1 - I - R + + tripartite(I,R,S); + +plt.hold(False) +plt.title('Simplex Plot of COVID-19 Pop Dynamics') + +vac = [1, 0.8, 0.6] +for loop in vac: + + + # Run the epidemic to the first peak + del1 = 0.005 + del2 = 0.005 + tlim = 52 + param = (mu, betap, nu) + t1, y1 = solve_flow(param) + + # Now vaccinate a fraction of the population + st = np.size(t1) + del1 = y1[st-1,0] + del2 = y1[st-1,1] + tlim = 400 + + param = (mu, loop*betap, nu) + t2, y2 = solve_flow(param) + + t2 = t2 + t1[st-1] + + tc = np.concatenate((t1,t2)) + yc = np.concatenate((y1,y2)) + + I = yc[:,0] + R = yc[:,1] + S = 1 - I - R + + plt.figure(3) + plt.hold(True) + lines = plt.semilogy(tc,I,tc,S,tc,R) + plt.ylim([0.001,1]) + plt.xlim([0,tlim]) + plt.legend(('Infected','Susceptible','Recovered')) + plt.setp(lines, linewidth=0.5) + plt.xlabel('Days') + plt.ylabel('Fraction of Population') + plt.title('Vaccination at 1 Year') + plt.show() + +plt.hold(False) \ No newline at end of file diff --git a/SIRWave.py b/SIRWave.py new file mode 100644 index 0000000..0bf018a --- /dev/null +++ b/SIRWave.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sat March 21 2020 + +@author: nolte + +D. D. Nolte, Introduction to Modern Dynamics: Chaos, Networks, Space and Time, 2nd ed. (Oxford,2019) + +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +print(' ') +print('SIR.py') + +def solve_flow(param,max_time=1000.0): + + def flow_deriv(x_y_z_w,tspan): + In, Sn, Iq, Sq = x_y_z_w + + Inp = -mu*In + beta*knn*In*Sn + beta*knq*Iq*Sn + Snp = -beta*knn*In*Sn - beta*knq*Iq*Sn + + Iqp = -mu*Iq + beta*kqn*In*Sq + beta*kqq*Iq*Sq + Sqp = -beta*kqn*In*Sq - beta*kqq*Iq*Sq + + return [Inp, Snp, Iqp, Sqp] + + x0 = [In0, Sn0, Iq0, Sq0] + + # Solve for the trajectories + t = np.linspace(tlo, thi, thi-tlo) + x_t = integrate.odeint(flow_deriv, x0, t) + + + return t, x_t + +beta = 0.03 # infection rate +dill = 6 # mean days infectious +mu = 1/dill # decay rate +fnq = 0.2 # fraction not quarantining +fq = 1-fnq # fraction quarantining +P = 330 # Population of the US in millions +mr = 0.001 # Mortality rate +dq = 90 # Days of lock-down (this is the key parameter) + +# During quarantine +knn = 50 # Average connections per day for non-compliant group among themselves +kqq = 0 # Connections among compliant group +knq = 0 # Effect of compliaht group on non-compliant +kqn = 5 # Effect of non-clmpliant group on compliant + +initfrac = 0.0001 # Initial conditions: +In0 = initfrac*fnq # infected non-compliant +Sn0 = (1-initfrac)*fnq # susceptible non-compliant +Iq0 = initfrac*fq # infected compliant +Sq0 = (1-initfrac)*fq # susceptivle compliant + +tlo = 0 +thi = dq + +param = (mu, beta, knn, knq, kqn, kqq) # flow parameters + +t1, y1 = solve_flow(param) + +In1 = y1[:,0] +Sn1 = y1[:,1] +Rn1 = fnq - In1 - Sn1 +Iq1 = y1[:,2] +Sq1 = y1[:,3] +Rq1 = fq - Iq1 - Sq1 + +# Lift the quarantine +knn = 50 +kqq = 5 +knq = 20 +kqn = 15 + +fin1 = len(t1) +In0 = In1[fin1-1] +Sn0 = Sn1[fin1-1] +Iq0 = Iq1[fin1-1] +Sq0 = Sq1[fin1-1] + +tlo = fin1 +thi = fin1 + 365-dq + +param = (mu, beta, knn, knq, kqn, kqq) + +t2, y2 = solve_flow(param) + +In2 = y2[:,0] +Sn2 = y2[:,1] +Rn2 = fnq - In2 - Sn2 +Iq2 = y2[:,2] +Sq2 = y2[:,3] +Rq2 = fq - Iq2 - Sq2 + +fin2 = len(t2) +t = np.zeros(shape=(fin1+fin2,)) +In = np.zeros(shape=(fin1+fin2,)) +Sn = np.zeros(shape=(fin1+fin2,)) +Rn = np.zeros(shape=(fin1+fin2,)) +Iq = np.zeros(shape=(fin1+fin2,)) +Sq = np.zeros(shape=(fin1+fin2,)) +Rq = np.zeros(shape=(fin1+fin2,)) + +t[0:fin1] = t1 +In[0:fin1] = In1 +Sn[0:fin1] = Sn1 +Rn[0:fin1] = Rn1 +Iq[0:fin1] = Iq1 +Sq[0:fin1] = Sq1 +Rq[0:fin1] = Rq1 + + +t[fin1:fin1+fin2] = t2 +In[fin1:fin1+fin2] = In2 +Sn[fin1:fin1+fin2] = Sn2 +Rn[fin1:fin1+fin2] = Rn2 +Iq[fin1:fin1+fin2] = Iq2 +Sq[fin1:fin1+fin2] = Sq2 +Rq[fin1:fin1+fin2] = Rq2 + +plt.figure(1) +lines = plt.semilogy(t,In,t,Iq,t,(In+Iq)) +plt.ylim([0.0001,.1]) +plt.xlim([0,thi]) +plt.legend(('Non-compliant','Compliant','Total')) +#plt.setp(lines, linewidth=0.5) +plt.xlabel('Days') +plt.ylabel('Infected (Millions)') +plt.title('Infection Dynamics for COVID-19 in US') +plt.show() + +plt.figure(2) +lines = plt.semilogy(t,Rn*P*mr,t,Rq*P*mr) +plt.ylim([0.001,1]) +plt.xlim([0,thi]) +plt.legend(('Non-compliant','Compliant')) +plt.setp(lines, linewidth=0.5) +plt.xlabel('Days') +plt.ylabel('Deaths (Millions)') +plt.title('Total Deaths for COVID-19 in US') +plt.show() +D = P*mr*(Rn[fin1+fin2-1] + Rq[fin1+fin2-1]) +print('Deaths = ',D) + +plt.figure(3) +lines = plt.semilogy(t,In/fnq,t,Iq/fq) +plt.ylim([0.0001,.1]) +plt.xlim([0,thi]) +plt.legend(('Non-compliant','Compliant')) +plt.setp(lines, linewidth=0.5) +plt.xlabel('Days') +plt.ylabel('Fraction of Sub-Population') +plt.title('Population Dynamics for COVID-19 in US') +plt.show() + diff --git a/Standmap.py b/Standmap.py new file mode 100644 index 0000000..5c9618a --- /dev/null +++ b/Standmap.py @@ -0,0 +1,60 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed May 2 16:17:27 2018 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Standard map +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +eps = 0.97 + +np.random.seed(2) + +plt.figure(1) + +for eloop in range(0,200): + + rlast = 2*np.pi*(0.5-np.random.random()) + thlast = 2*np.pi*np.random.random() + + # rold = 2.0*pi*(0.5-rand); + # thetold = 2.0*pi*rand; + + + rplot = np.zeros(shape=(200,)) + thetaplot = np.zeros(shape=(200,)) + for loop in range(0,200): + rnew = rlast + eps*np.sin(thlast) + thnew = np.mod(thlast+rnew,2*np.pi) + + thetaplot[loop] = np.mod(thnew-np.pi,2*np.pi) - np.pi + if rnew > np.pi: + rtemp = rnew-2*np.pi + elif rnew < -np.pi: + rtemp = 2*np.pi + rnew + else: + rtemp = rnew + + rplot[loop] = np.mod(0.5 + (rtemp + np.pi)/2/np.pi,1) + + + rlast = rnew + thlast = thnew + + plt.plot(np.real(thetaplot),np.real(rplot),'o',ms=1) +# plt.xlim(xmin=-np.pi,xmax=np.pi) +# plt.ylim(ymin=-np.pi,ymax=2*np.pi) + + +plt.savefig('StandMap') + + diff --git a/StandmapHom.py b/StandmapHom.py new file mode 100644 index 0000000..5173367 --- /dev/null +++ b/StandmapHom.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun Aug 2 2020 + +"Introduction to Modern Dynamics" 2nd Edition (Oxford, 2019) + +@author: nolte +""" + +import numpy as np +from matplotlib import pyplot as plt +from numpy import linalg as LA + +plt.close('all') + +eps = 0.97 + +np.random.seed(2) + +plt.figure(1) + +for eloop in range(0,100): + + rlast = 2*np.pi*(0.5-np.random.random()) + thlast = 4*np.pi*np.random.random() + + + rplot = np.zeros(shape=(200,)) + thetaplot = np.zeros(shape=(200,)) + for loop in range(0,200): + rnew = rlast + eps*np.sin(thlast) + thnew = np.mod(thlast+rnew,4*np.pi) + + thetaplot[loop] = np.mod(thnew-np.pi,4*np.pi) + + rtemp = np.mod(rnew + np.pi,2*np.pi) + + rplot[loop] = rtemp - np.pi + + + rlast = rnew + thlast = thnew + + plt.plot(np.real(thetaplot),np.real(rplot),'o',ms=.2) + plt.xlim(xmin=np.pi,xmax=4*np.pi) + plt.ylim(ymin=-2.5,ymax=2.5) + + +plt.savefig('StandMap') + +K = eps +eps0 = 5e-7 + +J = [[1,1+K],[1,1]] +w, v = LA.eig(J) + +My = w[0] +Vu = v[:,0] # unstable manifold +Vs = v[:,1] # stable manifold + +# Plot the unstable manifold +Hr = np.zeros(shape=(100,150)) +Ht = np.zeros(shape=(100,150)) +for eloop in range(0,100): + + eps = eps0*eloop + + roldu1 = eps*Vu[0] + thetoldu1 = eps*Vu[1] + + Nloop = np.ceil(-6*np.log(eps0)/np.log(eloop+2)) + flag = 1 + cnt = 0 + + while flag==1 and cnt < Nloop: + + + ru1 = roldu1 + K*np.sin(thetoldu1) + thetau1 = thetoldu1 + ru1 + + roldu1 = ru1 + thetoldu1 = thetau1 + + if thetau1 > 4*np.pi: + flag = 0 + + Hr[eloop,cnt] = roldu1 + Ht[eloop,cnt] = thetoldu1 + 3*np.pi + cnt = cnt+1 + + +lnwid = 0.75 + +x = Ht[0:99,12] - 2*np.pi +x2 = 6*np.pi - x +y = Hr[0:99,12] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[5:39,15] - 2*np.pi +x2 = 6*np.pi - x +y = Hr[5:39,15] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[12:69,16] - 2*np.pi +x2 = 6*np.pi - x +y = Hr[12:69,16] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[15:89,17] - 2*np.pi +x2 = 6*np.pi - x +y = Hr[15:89,17] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[30:99,18] - 2*np.pi +x2 = 6*np.pi - x +y = Hr[30:99,18] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +# Plot the stable manifold +del Hr, Ht +Hr = np.zeros(shape=(100,150)) +Ht = np.zeros(shape=(100,150)) +#eps0 = 0.03 +for eloop in range(0,100): + + eps = eps0*eloop + + roldu1 = eps*Vs[0] + thetoldu1 = eps*Vs[1] + + Nloop = np.ceil(-6*np.log(eps0)/np.log(eloop+2)) + flag = 1 + cnt = 0 + + while flag==1 and cnt < Nloop: + + + thetau1 = thetoldu1 - roldu1 + ru1 = roldu1 - K*np.sin(thetau1) + + + roldu1 = ru1 + thetoldu1 = thetau1 + + if thetau1 > 4*np.pi: + flag = 0 + + Hr[eloop,cnt] = roldu1 + Ht[eloop,cnt] = thetoldu1 + cnt = cnt+1 + +x = Ht[0:79,12] + np.pi +x2 = 6*np.pi - x +y = Hr[0:79,12] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[4:39,15] + np.pi +x2 = 6*np.pi - x +y = Hr[4:39,15] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[12:69,16] + np.pi +x2 = 6*np.pi - x +y = Hr[12:69,16] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[15:89,17] + np.pi +x2 = 6*np.pi - x +y = Hr[15:89,17] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + +del x,y +x = Ht[30:99,18] + np.pi +x2 = 6*np.pi - x +y = Hr[30:99,18] +y2 = -y +plt.plot(x,y,linewidth =lnwid) +plt.plot(x2,y2,linewidth =lnwid) + + diff --git a/Standmaptwist.py b/Standmaptwist.py new file mode 100644 index 0000000..50707dc --- /dev/null +++ b/Standmaptwist.py @@ -0,0 +1,49 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed Oct. 2, 2019 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Twist map +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +eps = 0.9 + +np.random.seed(2) + +plt.figure(1) + +for eloop in range(0,50): + + rlast = np.pi*(1.5*np.random.random()-0.5) + thlast = 2*np.pi*np.random.random() + + orbit = np.int(200*(rlast+np.pi/2)) + rplot = np.zeros(shape=(orbit,)) + thetaplot = np.zeros(shape=(orbit,)) + x = np.zeros(shape=(orbit,)) + y = np.zeros(shape=(orbit,)) + for loop in range(0,orbit): + rnew = rlast + eps*np.sin(thlast) + thnew = np.mod(thlast+rnew,2*np.pi) + + rplot[loop] = rnew + thetaplot[loop] = np.mod(thnew-np.pi,2*np.pi) - np.pi + + rlast = rnew + thlast = thnew + + x[loop] = (rnew+np.pi+0.25)*np.cos(thnew) + y[loop] = (rnew+np.pi+0.25)*np.sin(thnew) + + plt.plot(x,y,'o',ms=1) + +plt.savefig('StandMap') diff --git a/UserFunction.py b/UserFunction.py new file mode 100644 index 0000000..036acda --- /dev/null +++ b/UserFunction.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Sun May 12 12:33:27 2019 + +@author: nolte +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt +import networkx as nx + +def linfit(x,y): + + meanx = np.mean(x) + meany = np.mean(y) + meanxy = np.mean(x*y) + meanx2 = np.mean(x**2) + + m = (meanxy - meanx*meany)/(meanx2 - meanx**2) + b = meany - m*meanx + + return m, b + + + + \ No newline at end of file diff --git a/WebMap.py b/WebMap.py new file mode 100644 index 0000000..d8b909e --- /dev/null +++ b/WebMap.py @@ -0,0 +1,47 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Wed May 2 16:17:27 2018 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Web maps +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') +phi = (1+np.sqrt(5))/2 +K = 1-phi # (0.618, 4) (0.618,5) (0.618,7) (1.2, 4) +q = 4 # 4, 5, 6, 7 +alpha = 2*np.pi/q + +np.random.seed(2) + +plt.figure(1) + +for eloop in range(0,1000): + + xlast = 50*np.random.random() + ylast = 50*np.random.random() + + xnew = np.zeros(shape=(300,)) + ynew = np.zeros(shape=(300,)) + for loop in range(0,300): + xnew[loop] = (xlast + K*np.sin(ylast))*np.cos(alpha) + ylast*np.sin(alpha) + ynew[loop] = -(xlast + K*np.sin(ylast))*np.sin(alpha) + ylast*np.cos(alpha) + xlast = xnew[loop] + ylast = ynew[loop] + + plt.plot(np.real(xnew),np.real(ynew),'o',ms=1) + plt.xlim(xmin=-60,xmax=60) + plt.ylim(ymin=-60,ymax=60) + + +plt.title('WebMap') +plt.savefig('WebMap') + + diff --git a/gravlens.py b/gravlens.py new file mode 100644 index 0000000..a167fa8 --- /dev/null +++ b/gravlens.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue May 28 11:50:24 2019 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Gravitational lensing and photon orbits +""" + +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') + +def create_circle(): + circle = plt.Circle((0,0), radius= 10, color = 'black') + return circle + +def show_shape(patch): + ax=plt.gca() + ax.add_patch(patch) + plt.axis('scaled') + plt.show() + +def refindex(x,y): + + A = 10 + eps = 1e-6 + + rp0 = np.sqrt(x**2 + y**2); + + n = 1/(1 - A/(rp0+eps)) + fac = np.abs((1-9*(A/rp0)**2/8)) # approx correction to Eikonal + nx = -fac*n**2*A*x/(rp0+eps)**3 + ny = -fac*n**2*A*y/(rp0+eps)**3 + + return [n,nx,ny] + +def flow_deriv(x_y_z,tspan): + x, y, z, w = x_y_z + + [n,nx,ny] = refindex(x,y) + + yp = np.zeros(shape=(4,)) + yp[0] = z/n + yp[1] = w/n + yp[2] = nx + yp[3] = ny + + return yp + +for loop in range(-5,30): + + xstart = -100 + ystart = -2.245 + 4*loop + print(ystart) + + [n,nx,ny] = refindex(xstart,ystart) + + + y0 = [xstart, ystart, n, 0] + + tspan = np.linspace(1,400,2000) + + y = integrate.odeint(flow_deriv, y0, tspan) + + xx = y[1:2000,0] + yy = y[1:2000,1] + + + plt.figure(1) + lines = plt.plot(xx,yy) + plt.setp(lines, linewidth=1) + plt.show() + plt.title('Photon Orbits') + +c = create_circle() +show_shape(c) +axes = plt.gca() +axes.set_xlim([-100,100]) +axes.set_ylim([-100,100]) + +xstart = 0 +ystart = 15 + +[n,nx,ny] = refindex(xstart,ystart) + +y0 = [xstart, ystart, n, 0] + +tspan = np.linspace(1,94,1000) + +y = integrate.odeint(flow_deriv, y0, tspan) + +xx = y[1:1000,0] +yy = y[1:1000,1] + +plt.figure(1) +lines = plt.plot(xx,yy) +plt.setp(lines, linewidth=2, color = 'black') +plt.show() diff --git a/logistic.py b/logistic.py new file mode 100644 index 0000000..3649a55 --- /dev/null +++ b/logistic.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Nov 17 15:51:35 2020 + +@author: Nolte +""" + +import numpy as np +from matplotlib import pyplot as plt + +plt.close('all') +print('logistic.py') + +# model_case 1 = logistic +# model_case 2 = cosine +# model_case 3 = cubed +model_case = int(input('Enter the Model Case (1-3)')) + +if model_case == 1: + ming = 2.9 + maxg = 4.0 + grep = 900 + ymin = 0 + ymax = 1 + xold = 0.231 + def logis(g,xin): + xout = g*xin*(1.0-xin) + return xout +elif model_case == 2: + ming = 1 + maxg = 10 + grep = 900 + ymin = -1.2 + ymax = 1.2 + xold = 0.231 + def logis(g,xin): + xout = np.cos(g*xin) + return xout +else: + ming = 1 + maxg = 3 + grep = 900 + xold = 0.231 + def logis(g,xin): + xout = g*xin*(1-xin**2) + return xout + + +rep = 64; + +plt.figure(1) +plt.axes(xlim=(ming, maxg), ylim=(ymin, ymax)) + +for gainloop in range(grep): + + g=(maxg-ming)*gainloop/grep + ming; + + for sloop in range(100): # settle-down loop + xnew = logis(g,xold) + xold = xnew + + x = np.zeros(shape = (rep)) + y = np.zeros(shape = (rep)) + ind = 0 + for rloop in range(rep-1): + xnew = logis(g,xold) + xold = xnew + + ind = ind + 1 + x[ind] = g + (maxg-ming)*(rloop-1)/(grep*rep); + y[ind] = xnew; + + + + lines = plt.plot(x,y,'o',ms=1) + plt.pause(0.05) + plt.show() + + + diff --git a/raysimple.py b/raysimple.py new file mode 100644 index 0000000..e58c507 --- /dev/null +++ b/raysimple.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue May 28 11:50:24 2019 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Eikonal equation simulations +""" + +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') + +# selection 1 = Gaussian +# selection 2 = Donut +selection = 1 + +print(' ') +print('raysimple.py') + +def refindex(x,y): + + if selection == 1: + + sig = 10 + + n = 1 + np.exp(-(x**2 + y**2)/2/sig**2) + nx = (-2*x/2/sig**2)*np.exp(-(x**2 + y**2)/2/sig**2) + ny = (-2*y/2/sig**2)*np.exp(-(x**2 + y**2)/2/sig**2) + + elif selection == 2: + + sig = 10; + r2 = (x**2 + y**2) + r1 = np.sqrt(r2) + np.expon = np.exp(-r2/2/sig**2) + + n = 1+0.3*r1*np.expon; + nx = 0.3*r1*(-2*x/2/sig**2)*np.expon + 0.3*np.expon*2*x/r1 + ny = 0.3*r1*(-2*y/2/sig**2)*np.expon + 0.3*np.expon*2*y/r1 + + + return [n,nx,ny] + + +def flow_deriv(x_y_z,tspan): + x, y, z, w = x_y_z + + n, nx, ny = refindex(x,y) + + yp = np.zeros(shape=(4,)) + yp[0] = z/n + yp[1] = w/n + yp[2] = nx + yp[3] = ny + + return yp + +V = np.zeros(shape=(100,100)) +for xloop in range(100): + xx = -20 + 40*xloop/100 + for yloop in range(100): + yy = -20 + 40*yloop/100 + n,nx,ny = refindex(xx,yy) + V[yloop,xloop] = n + +fig = plt.figure(1) +if selection == 1: + contr = plt.contourf(V,100, cmap=cm.coolwarm, vmin = 1, vmax = 2) +else: + contr = plt.contourf(V,100, cmap=cm.coolwarm, vmin = 1, vmax = 3) +fig.colorbar(contr, shrink=0.5, aspect=5) +fig = plt.show() + + +v1 = 0.707 # Change this initial condition +v2 = np.sqrt(1-v1**2) +y0 = [12, v1, 0, v2] # Change these initial conditions + +tspan = np.linspace(1,1700,1700) + +y = integrate.odeint(flow_deriv, y0, tspan) + +plt.figure(2) +lines = plt.plot(y[1:1550,0],y[1:1550,1]) +plt.setp(lines, linewidth=0.5) +plt.show() + diff --git a/trirep.py b/trirep.py new file mode 100644 index 0000000..119590d --- /dev/null +++ b/trirep.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +trirep.py +Created on Thu May 9 16:23:30 2019 + +@author: nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +Population dynamics +""" + +import numpy as np +from scipy import integrate +from matplotlib import pyplot as plt + +plt.close('all') + +def tripartite(x,y,z): + + sm = x + y + z + xp = x/sm + yp = y/sm + + f = np.sqrt(3)/2 + + y0 = f*xp + x0 = -0.5*xp - yp + 1; + + plt.figure(2) + lines = plt.plot(x0,y0) + plt.setp(lines, linewidth=0.5) + plt.plot([0, 1],[0, 0],'k',linewidth=1) + plt.plot([0, 0.5],[0, f],'k',linewidth=1) + plt.plot([1, 0.5],[0, f],'k',linewidth=1) + plt.show() + + +def solve_flow(y,tspan): + def flow_deriv(y, t0): + #"""Compute the time-derivative .""" + + f = np.zeros(shape=(N,)) + for iloop in range(N): + ftemp = 0 + for jloop in range(N): + ftemp = ftemp + A[iloop,jloop]*y[jloop] + f[iloop] = ftemp + + phitemp = phi0 # Can adjust this from 0 to 1 to stabilize (but Nth population is no longer independent) + for loop in range(N): + phitemp = phitemp + f[loop]*y[loop] + phi = phitemp + + yd = np.zeros(shape=(N,)) + for loop in range(N-1): + yd[loop] = y[loop]*(f[loop] - phi); + + if np.abs(phi0) < 0.01: # average fitness maintained at zero + yd[N-1] = y[N-1]*(f[N-1]-phi); + else: # non-zero average fitness + ydtemp = 0 + for loop in range(N-1): + ydtemp = ydtemp - yd[loop] + yd[N-1] = ydtemp + + return yd + + # Solve for the trajectories + t = np.linspace(0, tspan, 701) + x_t = integrate.odeint(flow_deriv,y,t) + return t, x_t + +# model_case 1 = zero diagonal +# model_case 2 = zero trace +# model_case 3 = asymmetric (zero trace) +print(' ') +print('trirep.py') +print('Case: 1 = antisymm zero diagonal') +print('Case: 2 = antisymm zero trace') +print('Case: 3 = random') +model_case = int(input('Enter the Model Case (1-3)')) + +N = 3 +asymm = 3 # 1 = zero diag (replicator eqn) 2 = zero trace (autocatylitic model) 3 = random (but zero trace) +phi0 = 0.001 # average fitness (positive number) damps oscillations +T = 100; + + +if model_case == 1: + Atemp = np.zeros(shape=(N,N)) + for yloop in range(N): + for xloop in range(yloop+1,N): + Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1)) + Atemp[xloop,yloop] = -Atemp[yloop,xloop] + +if model_case == 2: + Atemp = np.zeros(shape=(N,N)) + for yloop in range(N): + for xloop in range(yloop+1,N): + Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1)) + Atemp[xloop,yloop] = -Atemp[yloop,xloop] + Atemp[yloop,yloop] = 2*(0.5 - np.random.random(1)) + + tr = np.trace(Atemp) + A = Atemp + for yloop in range(N): + A[yloop,yloop] = Atemp[yloop,yloop] - tr/N + +else: + Atemp = np.zeros(shape=(N,N)) + for yloop in range(N): + for xloop in range(N): + Atemp[yloop,xloop] = 2*(0.5 - np.random.random(1)) + + tr = np.trace(Atemp) + A = Atemp + for yloop in range(N): + A[yloop,yloop] = Atemp[yloop,yloop] - tr/N + +plt.figure(3) +im = plt.matshow(A,3,cmap=plt.cm.get_cmap('seismic')) # hsv, seismic, bwr +cbar = im.figure.colorbar(im) + +M = 20 +delt = 1/M +ep = 0.01; + +tempx = np.zeros(shape = (3,)) +for xloop in range(M): + tempx[0] = delt*(xloop)+ep; + for yloop in range(M-xloop): + tempx[1] = delt*yloop+ep + tempx[2] = 1 - tempx[0] - tempx[1] + + x0 = tempx/np.sum(tempx); # initial populations + + tspan = 70 + t, x_t = solve_flow(x0,tspan) + + y1 = x_t[:,0] + y2 = x_t[:,1] + y3 = x_t[:,2] + + plt.figure(1) + lines = plt.plot(t,y1,t,y2,t,y3) + plt.setp(lines, linewidth=0.5) + plt.show() + plt.ylabel('X Position') + plt.xlabel('Time') + + tripartite(y1,y2,y3)