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. 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() Continue reading...