 CoCalc Public FilesAnalysing Chaotic Systems.ipynb
Author: Patrick Davies
Views : 133
Description: Python code for plotting the Bifurcation diagrams and the Lorenz Attractor System
Compute Environment: Ubuntu 20.04 (Default)
In [ ]:


In :
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#
xN = np.linspace(0.0, 1.0, 101)
r = 5
xNPlusOne = r * xN*(1 - xN)

#
plt.figure(figsize = (12, 8))
plt.title("x\u2099\u208a\u2081 as a function of x\u2099")
plt.plot(xN, xNPlusOne, color = 'r', linestyle = '-', marker = '.', markerfacecolor = 'b', markeredgecolor = 'b')
plt.xlabel("x\u2099")
plt.ylabel("x\u2099\u208a\u2081")
plt.grid(color = 'g')
plt.savefig('Initial Logistics Map Graph')
plt.show() In :
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#Dripping tap example - constants alpha and beta
Npoints = 100
Time = np.zeros(Npoints) #Create array for time between drops
Time = 0.99 #Initial time value - anything below 1, has no effect on outcome
xAxis = np.zeros(Npoints) #Create array to store iteration number for the graph
index = 0 #Set index counter to 0
alpha = 2#Growth factor - changeable
beta = 1#Decay factor - changeable

for index in range (0, Npoints - 1): #repeat for as long as index value is between 0 and 99
Time[index+1] = alpha*Time[index] * (1 - (beta*Time[index])) #Next time value using chaos theory equation and current time value
xAxis[index] = index #Record iteration number in array
index += 1 #Increase index value by 1, eventually will end the loop

print(Time[index])

#print(Time) #debugging print statement
#print(index) #debugging print statement
#Graph plotting block for time between drops after each iteration/drop
fig = plt.figure(figsize = (8, 6))
plt.title("Convergence over many iterations graph (X = 0.99)")
plt.xlabel("Number of iterations")
plt.ylabel("x\u2099")
plt.scatter(xAxis, Time, marker = '|') # Using scatter plot as it just shows the plotted points
#plt.plot(xAxis, Time, color = 'b', linestyle = '-')
plt.savefig('Convergence over many iterations graph')
plt.grid(color = 'g')


0.49999999999999994 In :
from math import isclose

Npoints = 100
Npoints_i = 10000

alpha_i = np.linspace(0, 4, Npoints_i)

index_i = 0
Population = np.zeros(Npoints) #Create array for time between drops
Population = 0.6 #Initial time value - anything below 1, has no effect on outcome
Final_pop = np.zeros(Npoints_i)
Final_pop_1 = np.zeros(Npoints_i)
Final_pop_2 = np.zeros(Npoints_i)
Final_pop_3 = np.zeros(Npoints_i)
beta = 1

for index_i in range (0, Npoints_i):
index = 0
Population = 0.6 #Initial time value - anything below 1, has no effect on outcome
for index in range (0, Npoints - 1): #repeat for as long as index value is between 0 and 99
Population[index+1] = alpha_i[index_i]*Population[index] * (1 - (beta*Population[index])) #Next population value using chaos theory equation and current population

index += 1 #Increase index value by 1, eventually will end the loop

Final_pop[index_i] = Population[index]
if isclose(Population[index-1], Population[index], abs_tol = 1e-3) == True: #Population[index-1] == Population[index]:
Final_pop_1[index_i] = - 1
#print(-1)
else:
Final_pop_1[index_i] = Population[index - 1]

if isclose(Population[index-2], Population[index], abs_tol = 1e-3) == True or isclose(Population[index-2], Population[index-1], abs_tol = 1e-3) == True:#Population[index-2] == (Population[index] or Population[index-1]):
Final_pop_2[index_i] = - 1
#print(-2)
else:
Final_pop_2[index_i] = Population[index - 2]

if (isclose(Population[index-3], Population[index], abs_tol = 1e-3) == True) or (isclose(Population[index-3], Population[index-2], abs_tol = 1e-3) == True) or (isclose(Population[index-3], Population[index-1], abs_tol = 1e-3) == True): #Population[index-3] == (Population[index] or Population[index-1] or Population[index-2]):
Final_pop_3[index_i] = - 1
#print(-3)
else:
Final_pop_3[index_i] = Population[index - 3]
index_i += 1

fig = plt.figure(figsize = (16, 9))
plt.title("How alpha affects final population")
plt.xlabel("Alpha")
plt.ylabel("Convergence Value")
plt.scatter(alpha_i, Final_pop, marker = '.', color = 'k', s = 0.25, alpha = .25) # Using scatter plot as it just shows the plotted points
plt.scatter(alpha_i, Final_pop_1, marker = '.', color = 'k', s = 0.25, alpha = .25) # Using scatter plot as it just shows the plotted points
plt.scatter(alpha_i, Final_pop_2, marker = '.', color = 'k', s = 0.25, alpha = .25) # Using scatter plot as it just shows the plotted points
plt.scatter(alpha_i, Final_pop_3, marker = '.', color = 'k', s = 0.25, alpha = .25) # Using scatter plot as it just shows the plotted points
#plt.plot(alpha_i, Final_pop, color = 'b', linestyle = '-')
#plt.plot(alpha_i, Final_pop_1, color = 'b', linestyle = '-')
#plt.plot(alpha_i, Final_pop_2, color = 'b', linestyle = '-')
#plt.plot(alpha_i, Final_pop_3, color = 'b', linestyle = '-')
plt.ylim(bottom = -0.25)
plt.grid(color = 'g')
plt.savefig('Bifurcation diagram 1', optimise = True)

<ipython-input-5-2f732b3ffc3d>:60: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "optimise" which is no longer supported as of 3.3 and will become an error two minor releases later plt.savefig('Bifurcation diagram 1', optimise = True) In [ ]:


In :
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
#
def func(r, xN):
return r * xN*(1 - xN)
#
nPoints = 10000
r = np.linspace(0, 4.0, nPoints)
xN = 0.00001
iterations = 1000
last = 100
#
fig, (ax1) = plt.subplots(1, 1, figsize=(16, 9),
sharex=True)
for i in range(iterations):
xN = func(r, xN)
if i >= (iterations - last):
ax1.plot(r, xN, ',k', alpha=.25)
ax1.set_xlim(0, 4)
ax1.set_title("How alpha affects final population")
plt.xlabel("Alpha")
plt.ylabel("Convergence Value")
plt.savefig("Bifurcation diagram 2") ## Lorenz Attractor

In :
import time
startTime = time.time()

from scipy.integrate import odeint
from matplotlib import animation

sigma = 10
rho = 28
beta = 8/3

pos_init1 = [1.0, 1.0, 1.0]
pos_init2 = [1.0, 1.0, -1.0]       #Defining the initial positions of the 2 'particles'#
t = np.linspace(0, 100, 100000)

def Traj(init, t):
x, y, z = init
return sigma * (y - x), (x*(rho-z))-y, (x*y)-(beta*z)             #Defining the L.A differential equations which are to be solved#

positions1 = odeint(Traj, pos_init1, t)                              #Solving D.Es using odeint method#
positions2 = odeint(Traj, pos_init2, t)

fig = plt.figure(figsize = [12, 12])
ax = fig.gca(projection = "3d")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.grid(False)

ax.plot(positions1[:, 0], positions1[:, 1], positions1[:, 2], alpha = 1, color = 'darkorange', linewidth = 1)
#ax.plot(positions2[:, 0], positions2[:, 1], positions2[:, 2], alpha = 1, color = 'darkblue', linewidth = 1)
ax.plot(positions1[0, 0], positions1[0, 1], positions1[0, 2], marker = '.', color = 'darkorange')
#ax.plot(positions2[0, 0], positions2[0, 1], positions2[0, 2],ed marker = '.', color = 'darkblue')
#ax.plot(positions1[3999, 0], positions1[3999, 1], positions1[3999, 2], marker = '.', color = 'red')
#ax.plot(positions2[3999, 0], positions2[3999, 1], positions2[3999, 2], marker = '.', color = 'black')
#plt.draw()
plt.show()
print(positions1[3999, 0], positions1[3999, 1], positions1[3999, 2])
print(positions2[3999, 0], positions2[3999, 1], positions2[3999, 2])
#fig.savefig('odeint.png', optimise = True)

print("Time to execute", (time.time() - startTime), "Seconds")              #Calculates run time#

fig3 = plt.figure(figsize = [12, 12])
plt.xlabel('Time (seconds)')
plt.ylabel('x position (m)')
plt.grid(True)
plt.plot(t, positions1[:, 0], alpha = 1, color = 'darkorange', linewidth = 1)
plt.plot(t, positions2[:, 0], alpha = 1, color = 'darkblue', linewidth = 1)                #Graph of x component vs time#
plt.xlim(0, 20)
plt.legend()
plt.savefig('xPositions.png', optimise = True)
plt.show()

fig4 = plt.figure(figsize = [12, 12])
plt.xlabel('Time (seconds)')
plt.ylabel('y position (m)')
plt.grid(True)
plt.plot(t, positions1[:, 1], alpha = 1, color = 'darkorange', linewidth = 1)
plt.plot(t, positions2[:, 1], alpha = 1, color = 'darkblue', linewidth = 1)                #Graph of y component vs time#
plt.xlim(0, 20)
plt.legend()
plt.savefig('yPositions.png', optimise = True)
plt.show()

fig5 = plt.figure(figsize = [12, 12])
plt.xlabel('Time (seconds)')
plt.ylabel('z position (m)')
plt.grid(True)
plt.plot(t, positions1[:, 2], alpha = 1, color = 'darkorange', linewidth = 1)
plt.plot(t, positions2[:, 2], alpha = 1, color = 'darkblue', linewidth = 1)                #Graph of y component vs time#
plt.xlim(0, 20)
plt.legend()
plt.savefig('zPositions.png', optimise = True)
plt.show()


<ipython-input-3-4cf088092ad9>:19: MatplotlibDeprecationWarning: Calling gca() with keyword arguments was deprecated in Matplotlib 3.4. Starting two minor releases later, gca() will take no keyword arguments. The gca() function should only be used to get the current axes, or if no axes exist, create new axes with default keyword arguments. To create a new axes with non-default arguments, use plt.axes() or plt.subplot(). ax = fig.gca(projection = "3d")