Associated Map Grid in Python

Question

I am trying to build a bifurcation diagram for a 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 ran into the problem of getting the desired output result due to the number of transients that I use. Can someone help me by cross-checking my code, which values ​​of nTransients should be selected or how many transients should be ignored?

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 tell me that the figure displayed by pylab at the end has either dots or lines as a marker, if these are lines, then how to get a graph with dots.

This is my output image for reference after using pixels:

Answer 1

It is not yet clear what your desired result is, but I assume that you are aiming for something similar to this image from Wikipedia :

Following this assumption, I gave him the best shot, but I assume that your equations (with boundary conditions, etc.) give you something that just doesn't look pretty pretty. Here is my result:

This plot alone may not look best, however, if you zoom in, you can see some beautiful details (this is right from the center of the plot, where the two shoulders of the bifurcation descend, meet, and then turn away again):

Note that I used horizontal lines with alpha = 0.1 (initially you used solid vertical lines, so the result looked pretty good).

---

The code!

I substantially modified your program to make it vectorized: I deleted the for loop through p , which made it all work almost instantly. This allowed me to use a much denser sample for p and allow me to draw 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 to explain the whole program to you: I used several standard numpy tricks, including broadcasting , but otherwise I did not change much. I did not change your LM function at all.

Please feel free to ask me in the comments if you have any questions! I am happy to explain the specifics you want to explain.

Note on transients and iterations: I hope now that the program runs much faster, you can try to play with these elements yourself. In my opinion, the number of transients seemed to decide how long the plot remained “deterministic." The number of iterations only increases the density of storylines, so increasing this outside the point did not seem clear to me.

I tried to increase the number of transients to 10,000. Here is my result from this experiment, for your reference: