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:

按照这个假设，我尽了最大努力，但我猜你的方程（带有边界条件等）会给你一些看起来不太漂亮的东西。这是我的结果：

这个情节本身可能看起来不是最好的，但是，如果你放大，你真的可以看到一些美丽的细节（这是从情节的中心，分叉的两条臂下来的地方，相遇，然后再次分支）：

请注意，我使用了 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 函数。

如果您有任何问题，请随时在 cmets 中向我提问！我很乐意解释您想要解释的细节。

关于瞬态和迭代的说明：希望现在程序运行得更快，您可以自己尝试使用这些元素。对我来说，瞬变的数量似乎决定了情节保持“确定性”多长时间。迭代次数只会增加情节线的密度，因此将其增加超过一个点对我来说似乎没有意义。

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