Update

FlipPhone, GravSynch and StandmapHom

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)
+
+

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')

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")