Coupled map lattice in Python

Question

I attempt to plot bifurcation diagram for following one-dimensional spatially extended system with boundary conditions

x[i,n+1] = (1-eps)*(r*x[i,n]*(1-x[i,n])) + 0.5*eps*( r*x[i-1,n]*(1-x[i-1,n]) + r*x[i+1,n]*(1-x[i+1,n])) + p

I am facing problem in getting desired output figure may be because of number of transients I am using. Can someone help me out by cross-checking my code, what values of nTransients should I choose or how many transients should I ignore ?

My Python code is as follows:

import numpy as np
from numpy import *
from pylab import *

L = 60      # no. of lattice sites
eps = 0.6   # diffusive coupling strength
r = 4.0     # control parameter r

np.random.seed(1010)
ic = np.random.uniform(0.1, 0.9, L) # random initial condition betn. (0,1)

nTransients = 900 # The iterates we'll throw away
nIterates = 1000 # This sets how much the attractor is filled in
nSteps = 400 # This sets how dense the bifurcation diagram will be

pLow = -0.4
pHigh = 0.0
pInc = (pHigh-pLow)/float(nSteps)

def LM(p, x):
    x_new = []
    for i in range(L):
        if i==0:
            x_new.append((1-eps)*(r*x[i]*(1-x[i])) + 0.5*eps*(r*x[L-1]*(1-x[L-1]) + r*x[i+1]*(1-x[i+1])) + p)
        elif i==L-1:
            x_new.append((1-eps)*(r*x[i]*(1-x[i])) + 0.5*eps*(r*x[i-1]*(1-x[i-1]) + r*x[0]*(1-x[0])) + p)
        elif i>0 and i<L-1:
            x_new.append((1-eps)*(r*x[i]*(1-x[i])) + 0.5*eps*(r*x[i-1]*(1-x[i-1]) + r*x[i+1]*(1-x[i+1])) + p)
    return x_new

for p in arange(pLow, pHigh, pInc):
    # set initial conditions
    state = ic
    # throw away the transient iterations
    for i in range(nTransients):
        state = LM(p, state)
    # now stote the next batch of iterates
    psweep = []     # store p values
    x = []          # store iterates
    for i in range(nIterates):
        state = LM(p, state)
        psweep.append(p)
        x.append(state[L/2-1])
    plot(psweep, x, 'k,')   # Plot the list of (r,x) pairs as pixels

xlabel('Pinning Strength p')
ylabel('X(L/2)')

# Display plot in window
show()

Can someone also tell me figure displayed by pylab in the end has either dots or lines as a marker, if it is lines then how to get plot with dots.

This is my output image for reference, after using pixels:

Answer 1

It still isn't clear exactly what your desired output is, but I'm guessing you're aiming for something that looks like this image from Wikipedia:

Going with that assumption, I gave it my best shot, but I'm guessing your equations (with the boundary conditions and so on) give you something that simply doesn't look quite that pretty. Here's my result:

This plot by itself may not look like the best thing ever, however, if you zoom in, you can really see some beautiful detail (this is right from the center of the plot, where the two arms of the bifurcation come down, meet, and then branch away again):

Note that I have used horizontal lines, with alpha=0.1 (originally you were using solid, vertical lines, which was why the result didn't look good).

---

The code!

I essentially modified your program a little to make it vectorized: I removed the for loop over p, which made the whole thing run almost instantaneously. This enabled me to use a much denser sampling for p, and allowed me to plot horizontal lines.

from __future__ import print_function, division

import numpy as np
import matplotlib.pyplot as plt

L = 60      # no. of lattice sites
eps = 0.6   # diffusive coupling strength
r = 4.0     # control parameter r

np.random.seed(1010)
ic = np.random.uniform(0.1, 0.9, L)  # random initial condition betn. (0,1)

nTransients = 100   # The iterates we'll throw away
nIterates = 100     # This sets how much the attractor is filled in
nSteps = 4000       # This sets how dense the bifurcation diagram will be

pLow = -0.4
pHigh = 0.0
pInc = (pHigh - pLow) / nSteps

def LM(p, x):
    x_new = np.empty(x.shape)
    for i in range(L):
        if i == 0:
            x_new[i] = ((1 - eps) * (r * x[i] * (1 - x[i])) + 0.5 * eps * (r * x[L - 1] * (1 - x[L - 1]) + r * x[i + 1] * (1 - x[i + 1])) + p)
        elif i == L - 1:
            x_new[i] = ((1 - eps) * (r * x[i] * (1 - x[i])) + 0.5 * eps * (r * x[i - 1] * (1 - x[i - 1]) + r * x[0] * (1 - x[0])) + p)
        elif i > 0 and i < L - 1:
            x_new[i] = ((1 - eps) * (r * x[i] * (1 - x[i])) + 0.5 * eps * (r * x[i - 1] * (1 - x[i - 1]) + r * x[i + 1] * (1 - x[i + 1])) + p)
    return x_new

p = np.arange(pLow, pHigh, pInc)
state = np.tile(ic[:, np.newaxis], (1, p.size))

# set initial conditions
# throw away the transient iterations
for i in range(nTransients):
    state = LM(p, state)

# now store the next batch of iterates
x = np.empty((p.size, nIterates))           # store iterates
for i in range(nIterates):
    state = LM(p, state)
    x[:, i] = state[L // 2 - 1]

# Plot the list of (r,x) pairs as pixels
plt.plot(p, x, c=(0, 0, 0, 0.1))
plt.xlabel('Pinning Strength p')
plt.ylabel('X(L/2)')

# Display plot in window
plt.show()

I don't want to try explaining the whole program to you: I've used a few standard numpy tricks, including broadcasting, but otherwise, I have not modified much. I've not modified your LM function at all.

Please don't hesitate to ask me in the comments if you have any questions! I'm happy to explain specifics that you want explained.

A note on transients and iterates: Hopefully, now that the program runs much faster, you can try playing with these elements yourself. To me, the number of transients seemed to decide for how long the plot remained "deterministic-looking". The number of iterates just increases the density of plot lines, so increasing this beyond a point didn't seem to make sense to me.

I tried increasing the number of transients all the way up till 10,000. Here's my result from that experiment, for your reference: