问题描述
我在 3d 图形上绘制了洛伦兹吸引子。我想通过在洛伦兹图上绘制代表不同增长率大小的不同颜色点来了解繁殖向量的增长率大小如何影响洛伦兹吸引子。
这是我目前的代码:
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,c = "g")
elif (3.2 < GR[k] <= 6.4):
plt.scatter(0.5*k,c = "b")
else:
plt.scatter(0.5*k,c = "r")
x1,z1 构成洛伦兹吸引子,GR 是繁殖向量的增长率,其中前 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,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,y1[50*k],z1[50*k],c = "r")
但是在 plt.scatter(0.5*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],t_eval = times)
Xpn = solve_ivp(lorenz_solveivp(),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])
dX = np.array([1,1])/(np.sqrt(3))
dt = 0.01
n = 8
ng = 1000
GR,t = growth_rate(X0,ng)
#plot lorenz function
x1,z1 = sol1.y
fig = plt.figure(figsize=(8,8)) #specify the size of plot
ax = fig.add_subplot(111,projection = '3d')
ax.plot(x1,z1)
# add growth rate markers
ax = fig.add_subplot
for k in range(100):
if (GR[k] <= 0):
plt.scatter(0.5*k,color = "y")
elif (0 < GR[k] <= 3.2):
plt.scatter(0.5*k,color = "g")
elif (3.2 < GR[k] <= 6.4):
plt.scatter(0.5*k,color = "b")
else:
plt.scatter(0.5*k,color = "r")
解决方法
我修复了你的代码。代码有两个主要问题
- 默认情况下
matplotlib
具有二维线图plt.plot
默认为二维。我通过直接在您的 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],t_eval = times)
Xpn = solve_ivp(lorenz_solveivp(),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])
dX = np.array([1,1])/(np.sqrt(3))
dt = 0.01
n = 8
ng = 1000
GR,t = growth_rate(X0,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,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()