问题描述
我想生成陈的超混沌序列。公式如下:Chen's hyperchaotic sequence equations
附上我写的代码。
import math
a = 36
b = 3
c = 28
d = 16
k = 0.2
def chen(x0,y0,z0,q0):
xdot = a * (y0 - x0)
ydot = (-x0 * z0) + (d * x0) + (c * y0) - q0
zdot = (x0 * y0) - (b * z0)
qdot = x0 + k
return xdot,ydot,zdot,qdot
def chaotic_seq(x0,q0,length):
for i in range(length):
xdot,qdot = chen(x0,q0)
if math.isnan(xdot) or math.isnan(ydot) or math.isnan(zdot) or math.isnan(qdot):
print(i)
x0 = xdot
y0 = ydot
z0 = zdot
q0 = qdot
if __name__ == '__main__':
x0 = 0.3
y0 = -0.4
z0 = 1.2
q0 = 1
length = 2048
chaotic_seq(x0=x0,y0=y0,z0=z0,q0=q0,length=length)
我面临的问题是,在 'i=11' 之后,所有值(xdot、ydot、zdot、qdot)都是 NaN。
解决方法
我认为你有几个问题。
一个是您的函数似乎很快就会遇到 Python float
的溢出错误,它返回 nan
值,因此您需要一种支持比默认值更高的精度值的数据类型float
是 Python 内置数据类型提供的。因此,您可能会考虑使用 numpy
库的 float128
数据类型(如下所示)——或者使用 decimal
模块(未显示)进行调查。
其次,点符号表示输入变量的变化率。例如,xdot
是微分表达式 dx/dt
的简写。
您可以添加一个时间增量变量(例如,t
),它可以小幅度地更改 x0
、y0
、z0
和 q0
的值增量,模拟它们各自的差异。
这是运行 2048 次迭代的脚本的修改版本:
#!/usr/bin/env python
import sys
import numpy as np
a = np.float128(36)
b = np.float128(3)
c = np.float128(28)
d = np.float128(16)
k = np.float128(0.2)
t = np.float128(0.001)
def chen(x0,y0,z0,q0):
xdot = a * (y0 - x0)
ydot = (-x0 * z0) + (d * x0) + (c * y0) - q0
zdot = (x0 * y0) - (b * z0)
qdot = q0 + k
return xdot,ydot,zdot,qdot
def chaotic_seq(x0,q0,length):
for i in range(length):
xdot,qdot = chen(x0,q0)
if np.isnan(xdot) or np.isnan(ydot) or np.isnan(zdot) or np.isnan(qdot):
raise OverflowError("Overflow in dot variable calculation")
x0 += t * xdot
y0 += t * ydot
z0 += t * zdot
q0 += t * qdot
sys.stdout.write('after: [{}] {}\t{}\t{}\t{}\n'.format(i,x0,q0))
if __name__ == '__main__':
x0 = np.float128(0.3)
y0 = np.float128(-0.4)
z0 = np.float128(1.2)
q0 = np.float128(1)
length = 2048
chaotic_seq(x0=x0,y0=y0,z0=z0,q0=q0,length=length)
,
您的代码远未实现您的愿望:您将不得不在某个时间点求解该微分方程,而在上述示例中您并没有这样做。这就解释了为什么您的值会迅速发散到无穷大,然后开始变成 NaN。
使用 scipy 求解微分方程,我们得到以下代码,结果看似令人满意:
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
a = 36
b = 3
c = 12
d = 7
k = 0.2
def chen(_,y):
x0,q0 = y
xdot = a * (y0 - x0)
ydot = (-x0 * z0) + (d * x0) + (c * y0) - q0
zdot = (x0 * y0) - (b * z0)
qdot = x0 + k
return xdot,qdot
def chaotic_seq(x0,length):
return solve_ivp(chen,[0,length],[x0,q0])
x0 = 0.3
y0 = -0.4
z0 = 1.2
q0 = 1
length = 50
sol = chaotic_seq(x0=x0,length=length).y
# plot x,y
plt.plot(sol[0],sol[1])