From 18f762d8e74534a37937c8bf72f7ae9df7631272 Mon Sep 17 00:00:00 2001 From: David Nolte Date: Fri, 8 Oct 2021 10:28:29 -0400 Subject: [PATCH] Update FlipPhone, GravSynch and StandmapHom --- FlipPhone.py | 68 ++++++++++++++++++ GravSynch.py | 190 +++++++++++++++++++++++++++++++++++++++++++++++++ StandmapHom.py | 4 +- 3 files changed, 260 insertions(+), 2 deletions(-) create mode 100644 FlipPhone.py create mode 100644 GravSynch.py diff --git a/FlipPhone.py b/FlipPhone.py new file mode 100644 index 0000000..20d6ea0 --- /dev/null +++ b/FlipPhone.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Thurs Oct 7 19:38:57 2021 + +@author: David Nolte +Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019) + +FlipPhone Example +""" +import numpy as np +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') + +I1 = 0.45 +I2 = 0.5 +I3 = 0.55 + +def solve_lorenz(param, max_time=300.0, angle=0.0): + +# Flip Phone + def flow_deriv(x_y_z, t0): + + x, y, z = x_y_z + + yp1 = ((I2-I3)/I1)*y*z; + yp2 = ((I3-I1)/I2)*z*x; + yp3 = ((I1-I2)/I3)*x*y; + + + return [yp1, yp2, yp3] + model_title = 'Flip Phone' + + t = np.linspace(0, max_time/4, int(250*max_time/4)) + + # Solve for trajectories + x0 = [[0.01,1,0.01]] + t = np.linspace(0, max_time, int(250*max_time)) + x_t = np.asarray([integrate.odeint(flow_deriv, x0i, t) + for x0i in x0]) + + x, y, z = x_t[0,:,:].T + lines = ax.plot(x, y, z, '-') + plt.setp(lines, linewidth=0.5) + + ax.view_init(30, angle) + plt.show() + plt.title(model_title) + plt.savefig('Flow3D') + + return t, x_t + +ax.set_xlim((-1.1, 1.1)) +ax.set_ylim((-1.1, 1.1)) +ax.set_zlim((-1.1, 1.1)) + +t, x_t = solve_lorenz(max_time,angle=30) + +plt.figure(2) +lines = plt.plot(t,x_t[0,:,0],t,x_t[0,:,1],t,x_t[0,:,2]) +plt.setp(lines, linewidth=1) + + diff --git a/GravSynch.py b/GravSynch.py new file mode 100644 index 0000000..172c297 --- /dev/null +++ b/GravSynch.py @@ -0,0 +1,190 @@ +#!/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 +from pathlib import Path +import time + +tstart = time.time() + +plt.close('all') + +Grav = 6.674*10**(-11) +M0 = 1.989*10**30 +M = .9*M0 +prM = 'Mo' + str(int(M/M0)) +c = 3*10**8 +Rs = 2*Grav*M/c**2 +#R = 19000 + 14959787 +R = 4*Rs +scale = 10**-7 + +Nfac = 25 # 25 +N = 20 # 50 + +# model_case 1 = complete graph (Kuramoto transition) + +facoef = .2 +nodecouple = nx.complete_graph(N) + +# function: omegout, yout = coupleN(G) +def coupleN(G): + # function: yd = flow_deriv(x_y) + def flow_deriv(y, t0): + + yp = np.zeros(shape=(N,)) + for omloop in range(N): + temp = omega[omloop] + linksz = G.nodes[omloop]['numlink'] + for cloop in range(linksz): + cindex = G.nodes[omloop]['link'][cloop] + k = G.nodes[omloop]['coupling'][cloop] + + temp = temp + k * np.sin(y[cindex] - y[omloop]) + + yp[omloop] = 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 + + for nodeloop in range(N): + omega[nodeloop] = G.nodes[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(1, 1000, 1000) + 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) + +l = [] +grav = [] +omegatemp = np.array(()) +h = 10 + +print(Rs/R/scale) +h_dil = h * np.sqrt(1-Rs/R) +lab = "_T" + "&S" +for i in range(N): + h_dil = (i*h) * np.sqrt(1 - Rs / (R+h_dil)) + #l.append(i*h_dil) + l.append(h_dil) + + grav.append(2*Grav*M/(c**2*R*scale*(1 + l[i]/R))) + omegatemp = np.append(omegatemp,(grav[i])) + +Nlink = N * (N - 1) // 2 +omega = np.zeros(shape=(N,)) + + +meanomega = np.mean(omegatemp) +omega = (omegatemp - meanomega) +print(omegatemp) +print(meanomega) +print(omega) +sto = np.std(omega) + +lnk = np.zeros(shape=(N,), dtype=int) +for loop in range(N): + nodecouple.nodes[loop]['element'] = omega[loop] + nodecouple.nodes[loop]['link'] = list(nx.neighbors(nodecouple, loop)) + nodecouple.nodes[loop]['numlink'] = np.size(list(nx.neighbors(nodecouple, loop))) + lnk[loop] = np.size(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 + for nodeloop in range(N): + nodecouple.nodes[nodeloop]['coupling'] = np.zeros(shape=(lnk[nodeloop],)) + for linkloop in range(lnk[nodeloop]): + nodecouple.nodes[nodeloop]['coupling'][linkloop] = fac + + facval[facloop] = fac * avgdegree*scale + + 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]/scale + +plt.figure(1) +lines = plt.plot(xx, yy) +plt.setp(lines, linewidth=0.5) +plt.xticks(fontsize = 8) +plt.xlabel('Couple Constant g') +plt.yticks(fontsize = 8) +plt.ylabel('Frequency') + +plt.title('Kuramoto Synchronization Transition') + +filename_out = Path("Plot_N" + str(N) + '_' + prM + lab + ".png") +plt.savefig('GravSynch.png',dpi = 300) +plt.show() + + +elapsed_time = time.time() - tstart +print('elapsed time = ', format(elapsed_time, '.2f'), 'secs') \ No newline at end of file diff --git a/StandmapHom.py b/StandmapHom.py index 5173367..dad562f 100644 --- a/StandmapHom.py +++ b/StandmapHom.py @@ -47,8 +47,6 @@ plt.ylim(ymin=-2.5,ymax=2.5) -plt.savefig('StandMap') - K = eps eps0 = 5e-7 @@ -203,4 +201,6 @@ plt.plot(x,y,linewidth =lnwid) plt.plot(x2,y2,linewidth =lnwid) +plt.savefig('StandMap.png',dpi=900) +#plt.savefig('StandMap.tiff',dpi = 500, format = "tiff")