python - matplotlib 中的 fork 图

Question

我正在尝试获取以下方程的分岔图:

(x 是 t 的函数)

如:

这是我的片段:

import numpy as np
import matplotlib.pyplot as plt

def pitch(r, x):
    return r * x + np.power(x,3)- np.power(x,5)

n = 10000
r = np.linspace(-200, 200, n)
iterations = 1000
last = 100

x = 0

for i in range(iterations):
    x = pitch(r,x)

if i >= (iterations - last):
    plt.plot(r,x, ',k', alpha=0.02)

plt.title("Bifurcation diagram")

plt.show()

但是生成的图并不是它应该的样子:

编辑:

这是我最近的尝试:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def pitch(s,x,r):
    x = s[0]
    dxdt = r * x + np.power(x,3)- np.power(x,5)
    return [dxdt]

t = np.linspace(0,100)
s0=[-50]

r = np.linspace(-200, 200)

for i in r:
    s = odeint(pitch,s0,t, args=(i,))
    plt.plot(s,i,',k', alpha=0.02)

plt.title("Bifurcation diagram")

plt.show()

出现此错误:

raise ValueError("x and y must have same first dimension") ValueError: x and y must have same first dimension

您能给我一些建议来解决这个问题吗？!

Answer 1

我找到了这篇文章的链接，并决定发表一些评论，这些评论可能对将来偶然发现它的人有所帮助。

我没有详分割析这个方程，但从第一眼就可以清楚地看出，当r时会发生一些有趣的事情。接近于0。

因此我们可以研究 r in [-10,10] 系统的行为

您正确使用odeint而不是使用自己编写的欧拉方法求解柯西问题。

这个方程有一个吸引子，因为它很快就会“忘记”初始条件并滑向吸引子，但吸引子的选择取决于我们从 0 开始的位置。大的正初始条件将滑向负吸引子，反之亦然，如 - x^5是定义整个行为的术语 x .

我们需要做的是，对于范围内的每个 r，在每个初始条件下解滑动到的吸引子处做一个标记。

我们首先创建一个 Canvas 来放置标记:

diagram = np.zeros((200,200))

然后对于 (r,s0) 的每个组合我们在 Canvas 上放置一个点 (r,s[-1]) .

这是完整的代码

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def pitch(s,x,r):
    x = s[0]
    dxdt = r * x + np.power(x,3)- np.power(x,5)
    return [dxdt]

t = np.arange(0,100,2)
s0=[-50]

N = 200 # Number of points along each side of the diagram
diagram = np.zeros((N,N))

rmin,rmax = -10,10
rrange = np.arange(rmin, rmax,(rmax-rmin)/N)   

smin,smax = -5.0,5.0
srange = np.arange(smin,smax,2*(smax-smin)/N)    

for i in rrange:
    for s0 in srange:
        s = odeint(pitch,[s0],t, args=(i,))
        imgx = int((i-rmin)*N/(rmax-rmin))
        imgy = int((s[-1]-smin)/(smax-smin)*N)
        imgx = min(N-1,max(0,imgx)) # make sure we stay 
        imgy = min(N-1,max(0,imgy)) # within the diagram bounds
        diagram[imgy,imgx] = 1

plt.title("Bifurcation diagram")
plt.imshow(np.flipud(diagram),cmap=cm.Greys,
           extent=[rmin,rmax,smin,smax],aspect=(rmax-rmin)/(smax-smin))
plt.xlabel("r")
plt.ylabel("x")
plt.show()

以及结果图

当您通过设置 (rmin,rmax) 放大到 0 附近的区域时至(-0.5,0.5)您可以看到该图的分支并非从 0 开始

相反，如原始帖子中绘制的图表所示，分支大约从 r=-0.25 开始。