Plotting the bifurcation diagram for Ikeda map

Question

I'm trying to plot the bifurcation diagram for Ikeda map. I wrote a code in Python to get the points of this diagram, but it seems that for $u > 1$ the points diverge and my code doesn't work anymore. The plot I'm getting from my code is presented below.

I would like to know whether this plot is correct or not. Furthermore, Is there any way to obtain the bifurcation diagram for $u > 1$?

My code is the following:

A = []
X = []
Y = []

steps = 100

x = np.random.rand()
y = np.random.rand()


aparameter = np.linspace(0,1,800)

for k in aparameter:
    i = 0
    while i < steps:
        t = 0.4 - 6.0/(1.0 + x*x + y*y)
        x1 = 1.0 + k*(x*np.cos(t) - y*np.sin(t))
        y1 = k*(x*np.sin(t) + y*np.cos(t))
        if i > steps - 200:
            A.append(k)
            X.append(x1)
            Y.append(y1)
        x = x1
        y = y1
        i = i + 1

# display plot
plt.subplots(figsize=(fig_sizes,fig_sizes)) 
plt.plot(A,X, ',', color='black', markersize = 0.1, label = '$x_n$')
plt.plot(A,Y, ',', color='blue', markersize = 0.1, label = '$y_n$')
plt.xlabel(r'$a$')
plt.ylabel(r'$x_n, y_n$')
plt.ylim([-4, 4])
plt.legend(loc="upper left")
plt.axis('on')
plt.show()

Thanks in advance!

Answer 1

The Ikeda Map can be expressed as,

$$ x_{n+1} = R + C_2[x_n \cos(t_n) - y_n \sin(t_n)] \\ y_{n+1} = C_2[x_n \sin(t_n) + y_n \cos(t_n)] \\ t_n = C_1 - \frac{C_3}{1 + x_n^2 +y_n^2} $$

In the plot in the original question I was considering $R = 1, C_1 = 0.4, C_3 = 6$ and $C_2$ as a free bifurcation parameter. Therefore, indeed the plot is correct, but it's not what I was looking for. The plot I need to analyze can be obtained by considering exactly the same program, just changing the free parameter to $C_3$ and setting $C_2 = 0.9$. Thefore the code becomes,

A = []
X = []
Y = []

steps = 1000

x = np.random.rand()
y = np.random.rand()


aparameter = np.linspace(6.6, 7.6, 1000)

for k in aparameter:
    i = 0
    while i < steps:
        t = 0.4 - k/(1.0 + x*x + y*y)
        x1 = 1.0 + 0.9*(x*np.cos(t) - y*np.sin(t))
        y1 = 0.9*(x*np.sin(t) + y*np.cos(t))
        if i > steps - 200:
            A.append(k)
            X.append(x1)
            Y.append(y1)
        x = x1
        y = y1
        i = i + 1

# display plot
plt.subplots(figsize=(7,7)) 
plt.plot(A,X, ',', color='black', markersize = 0.1)
plt.xlabel(r'$a$')
plt.ylabel(r'$x_n$')
plt.axis('on')
plt.show()


plt.subplots(figsize=(7,7)) 
plt.plot(A,Y, ',', color='black', markersize = 0.1)
plt.xlabel(r'$a$')
plt.ylabel(r'$y_n$')
plt.axis('on')
plt.show()

Furnishing the following plots:

From which it's possible to identify the periodic window.