I'm trying to produce a bifurcation diagram such as this one.
Typical bifurcation diagrams only include stable fixed points, which are attractive (nearby points get pulled towards the fixed point upon iteration). These are easy enough to compute: just iterate a bunch of times and it will converge to the fixed point(s).
I am, however, trying to include the unstable (repelling) fixed points. By their nature, they will not be found using the algorithm described above.
My idea was to loosen the condition for a point to be fixed. Namely, if your value of x is within "epsilon" of f(x), then simply return x. This doesn't produce the desired output.
Note that symbolic computation of the fixed point is possible for period 2 and even period 4 fixed points, however, not beyond that (to my capabilities). I’m thus looking for a numerical approximation, not a symbolic solution.
I have included my code below. Any ideas or feedback is appreciated! Thanks.
import numpy as np
from matplotlib import pyplot as plt
# Define the logistic map
epsilon = 1e-4
def lmap(x, r):
# "Beefing up" the fixed point:
# if x is in the epsilon neighbourhood of f(x), let it be a fixed point (return x)
if abs(r * x * (1 - x) - x) < epsilon:
return x
else:
return r * x * (1 - x)
# Iterate over the logistic map to find fixed points
def get_fixed_pts(f, x0, r, iterations = 100, keep = 10, decimals = 5):
x = x0
for i in range(iterations):
x = f(x, r)
if i >= (iterations - keep):
yield np.round(x, decimals)
# Setting up some params
x0 = 0.5
n_pts = 1000
r_lims = [2, 4]
fixed_points = []
corresponding_r_values = []
# Looping over r values to find the fixed point for every r value
for r in np.linspace(*r_lims, n_pts):
# Get the fixed points and remove duplicates.
# Then, find the cycle size (number of fixed points)
pts = set(get_fixed_pts(lmap, x0, r))
cycle_size = len(pts)
# Store values
fixed_points.extend(pts)
corresponding_r_values.extend([r] * cycle_size)
# Plot all found points
plt.figure(figsize=(16, 9))
plt.scatter(corresponding_r_values, fixed_points, s = 0.01, c = 'black', label = 'Stable fixed point')
plt.title('Bifurcation diagram of logistic map')
plt.xlabel('Parameter Value')
plt.ylabel('Fixed Points')
plt.xlim(r_lims)
plt.legend()
plt.show()
I'm trying to produce a bifurcation diagram such as this one.
Typical bifurcation diagrams only include stable fixed points, which are attractive (nearby points get pulled towards the fixed point upon iteration). These are easy enough to compute: just iterate a bunch of times and it will converge to the fixed point(s).
I am, however, trying to include the unstable (repelling) fixed points. By their nature, they will not be found using the algorithm described above.
My idea was to loosen the condition for a point to be fixed. Namely, if your value of x is within "epsilon" of f(x), then simply return x. This doesn't produce the desired output.
Note that symbolic computation of the fixed point is possible for period 2 and even period 4 fixed points, however, not beyond that (to my capabilities). I’m thus looking for a numerical approximation, not a symbolic solution.
I have included my code below. Any ideas or feedback is appreciated! Thanks.
import numpy as np
from matplotlib import pyplot as plt
# Define the logistic map
epsilon = 1e-4
def lmap(x, r):
# "Beefing up" the fixed point:
# if x is in the epsilon neighbourhood of f(x), let it be a fixed point (return x)
if abs(r * x * (1 - x) - x) < epsilon:
return x
else:
return r * x * (1 - x)
# Iterate over the logistic map to find fixed points
def get_fixed_pts(f, x0, r, iterations = 100, keep = 10, decimals = 5):
x = x0
for i in range(iterations):
x = f(x, r)
if i >= (iterations - keep):
yield np.round(x, decimals)
# Setting up some params
x0 = 0.5
n_pts = 1000
r_lims = [2, 4]
fixed_points = []
corresponding_r_values = []
# Looping over r values to find the fixed point for every r value
for r in np.linspace(*r_lims, n_pts):
# Get the fixed points and remove duplicates.
# Then, find the cycle size (number of fixed points)
pts = set(get_fixed_pts(lmap, x0, r))
cycle_size = len(pts)
# Store values
fixed_points.extend(pts)
corresponding_r_values.extend([r] * cycle_size)
# Plot all found points
plt.figure(figsize=(16, 9))
plt.scatter(corresponding_r_values, fixed_points, s = 0.01, c = 'black', label = 'Stable fixed point')
plt.title('Bifurcation diagram of logistic map')
plt.xlabel('Parameter Value')
plt.ylabel('Fixed Points')
plt.xlim(r_lims)
plt.legend()
plt.show()
1 Answer
Reset to default 0To produce a bifurcation diagram of the logistic map with unstable fixed points, you can compute fixed points analytically, since unstable points repel trajectories.
Use f(x) = r * x * (1 - x) to find fixed points x = 0 and x = (r-1)/r,
and classify them as stable |f'(x)| < 1 or unstable |f'(x)| > 1 using the derivative f'(x) = r * (1 - 2x)
EDIT
The new code finds unstable fixed and periodic points by solving f^n(x) - x = 0 using Newton's method with different starting guesses. It then checks the stability of each found point by calculating the derivative and classifying it based on whether the derivative's magnitude is less than or greater than 1.
The code below can detect fixed points for higher order periodic orbits (2, 4, etc.)
import numpy as np
from matplotlib import pyplot as plt
from scipy import optimize
def lmap(x, r):
return r * x * (1 - x)
def lmap_deriv(x, r):
return r * (1 - 2 * x)
def lmap_iterate(x, r, n):
result = x
for _ in range(n):
result = lmap(result, r)
return result
def fixed_point_eq(x, r, n):
return lmap_iterate(x, r, n) - x
def get_stable_points(r, x0=0.5, iterations=500, keep=100):
x = x0
points = []
for i in range(iterations):
x = lmap(x, r)
if i >= (iterations - keep):
points.append(x)
return np.unique(np.round(points, decimals=5))
def get_all_fixed_points(r, period, guesses=[0.1, 0.5, 0.9]):
fixed_points = []
for x0 in guesses:
# Use Newton’s method to find roots of f^n(x) - x = 0
try:
root = optimize.newton(lambda x: fixed_point_eq(x, r, period), x0)
if 0 <= root <= 1: # Ensure physical bounds
fixed_points.append(round(root, 5))
except RuntimeError:
pass # Skip if Newton fails to converge
return np.unique(fixed_points)
def classify_stability(x, r, period):
h = 1e-6
x_next = lmap_iterate(x, r, period)
deriv = (lmap_iterate(x + h, r, period) - x_next) / h
return 'stable' if abs(deriv) < 1 else 'unstable'
# Parameters
r_lims = [2, 4]
n_pts = 1000
max_period = 4 # Check periods up to 4
stable_points = []
unstable_points = []
r_stable = []
r_unstable = []
for r in np.linspace(*r_lims, n_pts):
stable_pts = get_stable_points(r)
r_stable.extend([r] * len(stable_pts))
stable_points.extend(stable_pts)
# Get all fixed points for periods 1, 2, 4
for period in [1, 2, 4]:
all_pts = get_all_fixed_points(r, period)
for pt in all_pts:
if pt not in stable_pts: # Exclude points already identified as stable
stability = classify_stability(pt, r, period)
if stability == 'unstable':
unstable_points.append(pt)
r_unstable.append(r)
plt.figure(figsize=(16, 9))
plt.scatter(r_stable, stable_points, s=0.5, c='black', label='Stable Points')
plt.scatter(r_unstable, unstable_points, s=1, c='red', label='Unstable Points')
plt.title('Bifurcation Diagram of Logistic Map (Stable and Unstable Points)')
plt.xlabel('Parameter r')
plt.ylabel('x')
plt.xlim(r_lims)
plt.ylim([0, 1])
plt.legend(loc='upper left')
plt.show()