Find the Mathematica code here
Find the Mathematica code here
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def lorenz_rhs(t,state,s,r,b):
# Break up state variable into 3:
x,y,z = state
dxdt = s*(y-x)
dydt = r*x - y - x*z
dzdt = x*y - b*z
# Assemble
return [dxdt,dydt,dzdt]
# Initial condition
y0 = [0.9, 0.0, 0.0]
# Time span for the solution (start time, end time)
t_span = (0, 50)
# Plot the solution at the following times:
tvals = np.linspace(0,t_span[1],10000)
# Solve the ODE
solution1 = solve_ivp(lorenz_rhs, # system of equations
t_span, # time span of integration
y0, # initial conditions
args=(10,28,8/3), # parameters
t_eval=tvals, # times at which to store solution
method='RK45' # method of integration
)
fig1 = plt.figure(dpi=100,figsize=(12,5))
# 3-d plot
ax1_1 = fig1.add_subplot(1,2,1, projection='3d')
ax1_1.plot(solution1.y[0], solution1.y[1], solution1.y[2], lw=0.5, color='black')
ax1_1.set_title("Solution in [x,y,z] space ")
# time-trace
ax1_2 = fig1.add_subplot(1, 2, 2)
ax1_2.set_aspect(1/2)
ax1_2.plot(solution1.t, solution1.y[2], lw=0.5, color='blue')
ax1_2.set_title("Variable 'z' against time")
plt.show()| Using the function $$z_{n+1} = 1 - | 2z_n - 1 | ^{1/2}$$ to model the data from Lorenz’s map. Find the Mathematica code here. |
Successive iteration of \(x_{n+1} = x_n r (1-x_n)\) for 50,0000 steps. Find the Mathematica code here.
Probability of finding a particular value of $x$; histogram or Probability Density Function
Run the code in this NextJournal notebook or download the Python file here
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
import matplotlib as mpl
def logisticMap(x,r):
# Applies the "map" function f in logistic map
return r*x*(1-x)
def iterateLogisticMap(x0,nsteps,r):
# Applies the logistic map recursively to initial
# point x0 for nsteps using parameter r, and
# returns an array of (r,x) values.
xlist = np.zeros([nsteps,2])
xlist[0,:] = [r,x0]
x = x0
for i in range(nsteps):
x = logisticMap(x,r)
xlist[i,:] = [r,x]
return xlist
def orbitDiagramPoints(rmin,rmax,N_r,N_x):
# Takes as input a range of r values over which to
# iterate the map forward by 2*N_x steps. The variable
# N_r controls the density of points in the r-direction.
# Then, it keeps the last N_x points. This ensures
# that transients die out and we are left with a "grid" of
# N_x points in the vertical direction and N_r points
# in the horizontal direction.
r_vals = np.linspace(rmin,rmax,N_r)
points_grid = np.zeros([len(r_vals),N_x,2])
for i,r in enumerate(r_vals):
points = iterateLogisticMap(rand(),2*N_x,r)
keep_points = points[-N_x:]
points_grid[i,:,:] = keep_points
return points_grid
def makeOrbitDiagram(rmin,rmax,N_r,xmin,xmax,Nsteps):
# Does not output anything; generates an orbit diagram.
points_list = orbitDiagramPoints(rmin,rmax,N_r,Nsteps).reshape(-1,2)
plt.figure(figsize=(6,6),dpi=300)
plt.scatter(points_list[:,0],points_list[:,1],s=0.1)
plt.xlim(rmin,rmax)
plt.ylim(xmin,xmax)
# plt.gca().set_aspect(0.8)
plt.xlabel("Parameter r")
plt.ylabel("Variable x")
plt.title("Orbit Diagram for Logistic Map")
plt.savefig(f"OrbitDiagram_x{xmin}-{xmax}r{rmin}-{rmax}.png")
plt.show()
print("Use this module to run the function: ")
print("makeOrbitDiagram(rmin,rmax,N_r,xmin,xmax,Nsteps)")