1

我在 3d 图上绘制了 Lorenz 吸引子。我想通过在洛伦兹图上绘制代表不同大小的增长率的不同颜色的点来了解繁殖向量的增长率大小如何影响洛伦兹吸引子。

这是我目前的代码:

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection = '3d')

情节洛伦兹

ax.plot(x1, y1, z1)

添加增长率标记

ax = fig.add_subplot

for k in range(100):
    if (GR[k] <= 0):
      plt.scatter(0.5*k, x1[50*k], c = "y")
    elif (0 < GR[k] <= 3.2):
      plt.scatter(0.5*k, x1[50*k], c = "g")
    elif (3.2 < GR[k] <= 6.4):
      plt.scatter(0.5*k, x1[50*k], c = "b")
    else:
      plt.scatter(0.5*k, x1[50*k], c = "r")

x1, y1, z1 组成了 lorenz 吸引子,GR 是 bred 向量的增长率,其中前 50 个是:

  [0.          10.282047    10.8977816    9.94731134   5.09550477
  -2.90817325  -8.55789949 -10.22519406  -7.08646881  -4.03173251
   0.32302345   2.48287221  -0.64753007  -1.22328369   1.14720494
   0.50083297  -1.24334573  -1.97221857   1.48577796   2.20605109
  -1.09659768  -0.82320336   1.23992983   0.32689335  -1.35888724
  -1.8668327    1.79410769   1.84711434  -1.38602027  -0.44126068
   1.28436189   0.27735059  -1.35896733  -1.81959438   1.87149091
   1.53278532  -1.54682835  -0.15104558   1.35899661   0.39353056
  -1.21200428  -1.86788144   1.69061062   1.31533289  -1.6250634
   0.01201846   1.5258175    0.71428205  -0.86708544  -1.95685686]

在该行plt.scatter(0.5*k, x1[50*k], c = "y"),出现此错误:

TypeError: loop of ufunc does not support argument 0 of type NoneType which has no callable sqrt method

我也尝试以这种方式绘制散点图:

ax = fig.add_subplot
for k in range(100):
    if (GR[k] <= 0):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], c = "y")
    elif (0 < GR[k] <= 3.2):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], c = "g")
    elif (3.2 < GR[k] <= 6.4):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], c = "b")
    else:
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], c = "r")

但是在这行plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], c = "y")出现以下错误:

TypeError: scatter() got multiple values for argument 'c'

有人能帮忙吗?

下面是一个应该可以在本地运行的 MIV:

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

#Initial Conditions
X0 = np.array([1, 1, 1])

#define lorenz function
def lorenz_solveivp(sigma=10, r=28, b=8/3):
    def rhs(t, X):
        return np.array([sigma*(X[1] - X[0]), -X[0]*X[2] + r*X[0] - X[1], X[0]*X[1] - b*X[2]])
    return rhs

#define time and apply solve_ivp
t = np.linspace(0, 50, 5001)  #equispaced points from 0 to 50, timestep of 0.01
rhs_function = lorenz_solveivp()
sol1 = solve_ivp(rhs_function, (0, 50), X0, t_eval=t)

#find the growth rates
def growth_rate(X0, dX, n, dt, ng):
    Xp0 = X0 + dX
    times = np.linspace(0, n*dt, n + 1)
    Xn_save = np.zeros((X0.size, n*(ng-1)+1))
    Xpn_save = np.zeros((Xp0.size, n*(ng-1)+1))
    t = np.zeros((n*(ng - 1) + 1))
    i = 0
    g = np.zeros(ng)

    while i < ng - 1:

        Xn = solve_ivp(lorenz_solveivp(), [0, n*dt], X0, t_eval = times)
        Xpn = solve_ivp(lorenz_solveivp(), [0, n*dt], Xp0, t_eval = times)

        Xn_save[:, n*i: n + 1 + n*i] = Xn.y
        Xpn_save[:, n*i: n + 1 + n*i] = Xpn.y
        t[n*i: n + 1 + n*i] = Xn.t + i*n*dt

        dXb = Xpn.y[:,n] - Xn.y[:,n]

        g[i+1] = np.log(np.linalg.norm(dXb)/np.linalg.norm(dX))/(n*dt)

        P_new = dXb*(np.linalg.norm(dX)/np.linalg.norm(dXb))
        Xp0 = Xn.y[:,n] + P_new
        X0 = Xn.y[:,n]
        i += 1

    return g, Xn_save, t

X0 = np.array([1, 1, 1])
dX = np.array([1, 1, 1])/(np.sqrt(3))
dt = 0.01
n = 8
ng = 1000
GR, Xn_save, t = growth_rate(X0, dX, n, dt, ng)

#plot lorenz function
x1, y1, z1 = sol1.y

fig = plt.figure(figsize=(8, 8)) #specify the size of plot
ax = fig.add_subplot(111, projection = '3d')
ax.plot(x1, y1, z1)

# add growth rate markers
ax = fig.add_subplot
for k in range(100):
    if (GR[k] <= 0):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], color = "y")
    elif (0 < GR[k] <= 3.2):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], color = "g")
    elif (3.2 < GR[k] <= 6.4):
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], color = "b")
    else:
      plt.scatter(0.5*k, x1[50*k], y1[50*k], z1[50*k], color = "r")
4

1 回答 1

0

我修复了你的代码。代码有两个主要问题

  • 默认情况下matplotlib,二维线图plt.plot默认为 2d。我通过直接在你的 3d 轴上绘制来解决这个问题
  • 删除0.5 * k,因为我不确定这是在做什么,但不是我认为的 3d 坐标系的一部分。

在此处输入图像描述

import numpy as np
from scipy.integrate import solve_ivp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

#Initial Conditions
X0 = np.array([1, 1, 1])

#define lorenz function
def lorenz_solveivp(sigma=10, r=28, b=8/3):
    def rhs(t, X):
        return np.array([sigma*(X[1] - X[0]), -X[0]*X[2] + r*X[0] - X[1], X[0]*X[1] - b*X[2]])
    return rhs

#define time and apply solve_ivp
t = np.linspace(0, 50, 5001)  #equispaced points from 0 to 50, timestep of 0.01
rhs_function = lorenz_solveivp()
sol1 = solve_ivp(rhs_function, (0, 50), X0, t_eval=t)

#find the growth rates
def growth_rate(X0, dX, n, dt, ng):
    Xp0 = X0 + dX
    times = np.linspace(0, n*dt, n + 1)
    Xn_save = np.zeros((X0.size, n*(ng-1)+1))
    Xpn_save = np.zeros((Xp0.size, n*(ng-1)+1))
    t = np.zeros((n*(ng - 1) + 1))
    i = 0
    g = np.zeros(ng)

    while i < ng - 1:

        Xn = solve_ivp(lorenz_solveivp(), [0, n*dt], X0, t_eval = times)
        Xpn = solve_ivp(lorenz_solveivp(), [0, n*dt], Xp0, t_eval = times)

        Xn_save[:, n*i: n + 1 + n*i] = Xn.y
        Xpn_save[:, n*i: n + 1 + n*i] = Xpn.y
        t[n*i: n + 1 + n*i] = Xn.t + i*n*dt

        dXb = Xpn.y[:,n] - Xn.y[:,n]

        g[i+1] = np.log(np.linalg.norm(dXb)/np.linalg.norm(dX))/(n*dt)

        P_new = dXb*(np.linalg.norm(dX)/np.linalg.norm(dXb))
        Xp0 = Xn.y[:,n] + P_new
        X0 = Xn.y[:,n]
        i += 1

    return g, Xn_save, t

X0 = np.array([1, 1, 1])
dX = np.array([1, 1, 1])/(np.sqrt(3))
dt = 0.01
n = 8
ng = 1000
GR, Xn_save, t = growth_rate(X0, dX, n, dt, ng)

#plot lorenz function
x1, y1, z1 = sol1.y

fig = plt.figure(figsize=(8, 8)) #specify the size of plot
ax = fig.add_subplot(111, projection = '3d')
ax.plot(x1, y1, z1)

# add growth rate markers
for k in range(100):
    # simplified expression to set color
    if (GR[k] <= 0):
        c = 'y'
    elif (0 < GR[k] <= 3.2):
        c = 'g'
    elif (3.2 < GR[k] <= 6.4):
        c = 'b'
    else:
        c = 'r'
    ax.scatter(x1[50*k], y1[50*k], z1[50*k], color = c)
fig.show()
于 2020-12-23T09:29:07.433 回答