نمای لیاپانوف (Lyapunov Exponent) چیست؟ — از صفر تا صد

1import numpy as np
2import matplotlib.pyplot as plt
3
4#Define the Logistic Equation. This will be used for most part of this algorithm
5#Let us use this time the form of the logistic equation without any coefficient at the beginning
6#compared to the previous exercise of using 4*r*x*(1-x)
7def Log(x,r):
8    return r*x*(1-x)
9
10#1. Lyapunov Exponent
11#The Lyapunov Exponent of a Logistic Map is given by the infinite lim of 1/n multiplied by the summation of ln(df(x_i)/dx) from
12#zero to n-1.
13
14#We begin by determining the derivative of the Logistic Map and getting it's natural logarithm
15def Lyapunov (x,r):
16    return np.log(abs(r - 2*r*x))
17
18#We then define our parameter in order for us to iterate the summation function later
19n = 10000
20r = np.linspace(2.5,4.0,n)
21iteration = 1000
22lyp = np.zeros(n)
23x = 1e-5*np.ones(n)
24
25#The following code allows us to create a summation for the Lyapunov Exponent
26for i in range (iteration):
27    x = Log(x,r)
28    lyp += Lyapunov(x,r)
29fig, ax2 = plt.subplots(1,1, figsize=(8,9), sharex=True)
30#Horizontal Line at y = 0
31ax2.axhline(0,color='k',lw=0.5,alpha=0.5)
32#Negative Values
33ax2.plot(r[lyp < 0], lyp[lyp < 0] / iteration, '.k', alpha=0.5, ms=.5)
34#Positive Values
35ax2.plot(r[lyp >= 0], lyp[lyp >= 0] / iteration, '.r', alpha=0.5, ms=0.5)
36ax2.set_xlim(2.5,4)
37ax2.set_ylim(-2,1)
38ax2.set_xlabel('r')
39ax2.set_ylabel('Lambda')
40ax2.set_title('Liapunov Exponents of the Logistic Map')
41
42#2. Bifurcation Diagram
43#The bifurcation diagram is a simple iteration of the logisitc map. It's simply a plot of all the possible r values and their corresponding X values
44rs = np.linspace(2.5,4,1000) #We define the range of values our paramter r will range. In this case we want it to iterate from 0 to 4 10,000 times
45x_n = [] #We place our (x,r) values
46for r in rs: #We iterate across all the r values from 0 to 4 
47    x = 0.00001 #We Set an Initial X value
48    for i in range (900):
49        x = Log(x,r)
50    for i in range (100):
51        x = Log(x,r)
52        x_n.append([x,r])
53
54x_n = np.array(x_n)
55fig, ax1 = plt.subplots(1, 1, figsize=(8,9), sharey=True)
56ax1.plot(x_n[:,1], x_n[:,0],',k',alpha=0.25)
57ax1.set_xlim(2.5,4)
58ax1.set_ylim(0,1)
59ax1.set_xlabel('r')
60ax1.set_ylabel('x_n')
61ax1.set_title('Bifurcation Diagram of the Logistic Map')
62   
63
64#3. Time Series
65#A way to determine if a function is chaotic is by checking it's Time Series to see if a distict pattern is visible
66#The time series is simply the values of our population at different time steps
67
68x = 0.5 #We define the an initial x population. I
69r = 3.828427125 #We Set an initial r condition first
70x_n = [] #We let this be a list in which all the iterated values of x will be placed
71t = [] #We let this be the iteration step. We define our iteration step to be the "time" component of our time series
72#Iterate the Logistic Map
73for i in range(50):
74    x = Log(x,r)
75    x_n.append(x)
76    t.append(i)
77    yaxis = np.array(x_n)
78    xaxis = np.array(t)
79    
80# Plot    
81fig, ax0 = plt.subplots(1, 1, figsize=(8, 9), sharey=True)
82plt.plot(xaxis,yaxis)
83plt.xlabel('t')
84plt.ylabel('x_n')
85plt.title('Time series of Logistic map with r = 3.828427125')
86
87#4. Cycles of the Logistic Map
88#An alternative way of determining if the logistic map is chaotic is by checking it's 
89
90#The function first defines the range and the logistic map as a function if it's parameter r. We then set an initial x value then
91#iterate the logistic map again while plotting values along the x and y coordinate
92def plot_system(x0, r, n, ax=None):
93    t = np.linspace(0, 1)
94    ax.plot(t, Log(t, r), 'k', lw=2)
95    ax.plot([0, 1], [0, 1], 'k', lw=2)
96
97    x = x0
98    for i in range(n):
99        y = Log(x, r)
100        # Plot the two lines.
101        ax.plot([x, x], [x, y], 'k', lw=1)
102        ax.plot([x, y], [y, y], 'k', lw=1)
103        # Position is plot with increasing opacity every time it overlaps itself
104        ax.plot([x], [y], 'ok', ms=10,
105                alpha=(i + 1) / n)
106        x = y
107
108    ax.set_xlim(0, 1)
109    ax.set_ylim(0, 1)
110
111#Plots
112fig, ax = plt.subplots(1, 1, figsize=(8, 9), sharey=True)
113plot_system(0.5, 3.828427125, 100, ax)
114plt.xlabel ('x')
115plt.ylabel('x_n')
116plt.title('x vs x_n of the Logistic map with r = 3.828427125')