In Python: 如何绘制洛伦兹系统在变参数值下的分岔图？

Question

所以，我已经在 Mathematica 中看到了我的问题的编码解决方案，但是由于对 Mathematica 的了解很少，我还没有能够重现它。

这就是我要用 Python 做的事情：https://mathematica.stackexchange.com/questions/159211/how-to-make-a-bifurcation-diagram-of-the-lorenz-system-under-a-varying-parameter

我认为我的错误在于理解如何计算我正在寻找的内容以及如何调整我的可视化以使其看起来像 link 中的那样，但欢迎提出任何想法。

我目前的代码如下所示：

def lorenz_system(x,y,z,r,s=10,b=6):
    x_dot = s*(y-x)
    y_dot = r*x-y-x*z
    z_dot = x*z-b*z
    return x_dot, y_dot, z_dot

dr = 0.1              # parameter step size
r=np.arange(40,200,dr)  # parameter range
dt = 0.001             # time step
t = np.arange(0,10,dt) # time range

#initialize solution arrays
xs = np.empty(len(t) + 1)
ys = np.empty(len(t) + 1)
zs = np.empty(len(t) + 1)

#initial values x0,y0,z0 for the system
xs[0], ys[0], zs[0] = (1, 1, 1)

for R in r:
    for i in range(len(t)):
        #approximate numerical solutions to system
        x_dot, y_dot, z_dot = lorenz_system(xs[i], ys[i], zs[i],R)
        xs[i + 1] = xs[i] + (x_dot * dt)
        ys[i + 1] = ys[i] + (y_dot * dt)
        zs[i + 1] = zs[i] + (z_dot * dt)
    #calculate and plot the peak values of the z solution
    for i in range(0,len(zs)-1):
        #using only the positive values in the z solution
        if zs[i]>0:
            #find the local maxima
            if (zs[i-1] < zs[i] and zs[i] > zs[i+1]):
                if (zs[i]<=1000):
                    #plot the local maxima point of the z solution that used the parameter R in r
                    plt.scatter(R,zs[i], color='black')
plt.xlim(0,200)        
plt.ylim(0,400)

Answer 1

lorenz_system函数有bug，应该是z_dot = x * y - b * z.

链接的答案还“使用一个运行的最终值作为下一个的初始条件，作为一种靠近吸引子的简单方法。”，并绘制局部最小值和局部最大值。

这是使用您的代码获得类似图的方法

import numpy as np
import matplotlib.pyplot as plt


def lorenz_system(x, y, z, r, b=10, s=6):
    x_dot = b * (y - x)
    y_dot = r * x - y - x * z
    z_dot = x * y - s * z
    return x_dot, y_dot, z_dot


dr = 0.1  # parameter step size
r = np.arange(40, 200, dr)  # parameter range
dt = 0.001  # time step
t = np.arange(0, 10, dt)  # time range

# initialize solution arrays
xs = np.empty(len(t) + 1)
ys = np.empty(len(t) + 1)
zs = np.empty(len(t) + 1)

# initial values x0,y0,z0 for the system
xs[0], ys[0], zs[0] = (1, 1, 1)


# Save the plot points coordinates and plot the with a single call to plt.plot
# instead of plotting them one at a time, as it's much more efficient
r_maxes = []
z_maxes = []
r_mins = []
z_mins = []


for R in r:
    # Print something to show everything is running
    print(f"{R=:.2f}")
    for i in range(len(t)):
        # approximate numerical solutions to system
        x_dot, y_dot, z_dot = lorenz_system(xs[i], ys[i], zs[i], R)
        xs[i + 1] = xs[i] + (x_dot * dt)
        ys[i + 1] = ys[i] + (y_dot * dt)
        zs[i + 1] = zs[i] + (z_dot * dt)
    # calculate and save the peak values of the z solution
    for i in range(1, len(zs) - 1):
        # save the local maxima
        if zs[i - 1] < zs[i] and zs[i] > zs[i + 1]:
            r_maxes.append(R)
            z_maxes.append(zs[i])
        # save the local minima
        elif zs[i - 1] > zs[i] and zs[i] < zs[i + 1]:
            r_mins.append(R)
            z_mins.append(zs[i])

    # "use final values from one run as initial conditions for the next to stay near the attractor"
    xs[0], ys[0], zs[0] = xs[i], ys[i], zs[i]


plt.scatter(r_maxes, z_maxes, color="black", s=0.5, alpha=0.2)
plt.scatter(r_mins, z_mins, color="red", s=0.5, alpha=0.2)

plt.xlim(0, 200)
plt.ylim(0, 400)
plt.show()

结果：

In Python: 如何绘制洛伦兹系统在变参数值下的分岔图？

In Python: How to make a bifurcation diagram of the Lorenz system under a varying parameter value?

python

physics

chaos

lorenz-system