python - Python 中的耦合映射点阵

Question

我尝试绘制具有边界条件的一维空间扩展系统的分岔图

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

我在获得所需输出数字时遇到问题可能是因为我使用的瞬变数。有人可以通过交叉检查我的代码来帮助我吗？我应该选择什么 nTransients 值或者我应该忽略多少瞬变？

我的Python代码如下:

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()

有人能告诉我pylab最后显示的图形有点还是线作为标记，如果是线那么如何用点绘制。

这是我的引用输出图像，使用像素后:

Answer 1

目前还不清楚您想要的输出到底是什么，但我猜您的目标是看起来像这张图片 from Wikipedia 的东西。 :

按照这个假设，我尽了最大的努力，但我猜你的方程式(带有边界条件等)给你的东西看起来并不那么漂亮。这是我的结果:

这个图本身可能看起来不是最好的东西，但是，如果你放大，你真的可以看到一些美丽的细节(这是从图的中心开始， fork 的两条臂向下延伸的地方, 相遇, 然后再次 fork ):

请注意，我使用了水平线，alpha=0.1(最初您使用的是实心垂直线，这就是结果看起来不太好的原因)。

---

代码!

我基本上对您的程序进行了一些修改以使其矢量化:我删除了 p 上的 for 循环，这使得整个程序几乎立即运行。这使我能够对 p 使用更密集的采样，并允许我绘制水平线。

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()

我不想向您解释整个程序:我使用了一些标准的 numpy 技巧，包括 broadcasting , 但除此之外，我没有做太多修改。我根本没有修改您的 LM 函数。

如果您有任何问题，请随时在评论中问我!我很乐意解释您想要解释的细节。

关于 transient 和迭代的注意事项:希望现在程序运行得更快，您可以尝试自己玩弄这些元素。对我来说， transient 的数量似乎决定了情节保持“确定性外观”的时间。迭代次数只会增加绘图线的密度，因此将其增加到一个点以上对我来说似乎没有意义。

我尝试将瞬变的数量一直增加到 10,000。这是我的实验结果，供您引用: