4D 混沌系统 Lyapunov 指数

Question

我正在尝试研究 4 维混沌吸引子 Lyapunov 谱，到目前为止，下面提到的代码适用于三维系统，但在 4D 和 5D 系统中会出现错误

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

def diff_Lorenz(u):
    x,y,z,w= u
    f = [a*(y-x) , x*z+w, b-x*y, z*y-c*w]
    Df = [[-a,a,0,0], [z,0, x,1], [-y, -x, 0,0],[0,z,y,-c]]
    return np.array(f), np.array(Df)


def LEC_system(u):
    #x,y,z = u[:3]
    U = u[2:18].reshape([4,4])
    L = u[12:15]
    f,Df = diff_Lorenz(u[:4])
    A = U.T.dot(Df.dot(U))
    dL = np.diag(A).copy();
    for i in range(4):
        A[i,i] = 0
        for j in range(i+1,4): A[i,j] = -A[j,i]
    dU = U.dot(A)
    return np.concatenate([f,dU.flatten(),dL])


a=6;b=11;c=5;

u0 = np.ones(4)
U0 = np.identity(4)
L0 = np.zeros(4)
u0 = np.concatenate([u0, U0.flatten(), L0])
t = np.linspace(0,10,301)
u = odeint(lambda u,t:LEC_system(u),u0,t, hmax=0.05)
L = u[5:,12:15].T/t[5:]

# plt.plot(t[5:],L.T) 
# plt.show()
p1=L[0,:];p2=L[1,:];p3=L[2,:];p4=L[3,:]
L1 = np.mean(L[0,:]);L2=np.average(L[1,:]);L3=np.average(L[2,:]);L4=np.average(L[3,:])
t1 = np.linspace(0,100,len(p1))
plt.plot(t1,p1);plt.plot(t1,p2);plt.plot(t1,p3);plt.plot(t1,p4)

# plt.show()
print('LES= ',L1,L2,L3,L4)

输出错误为

D:\anaconda3\lib\site-packages\scipy\integrate\odepack.py:247: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
~\AppData\Local\Temp/ipykernel_7008/1971199288.py in <module>
     32 # plt.plot(t[5:],L.T)
     33 # plt.show()
---> 34 p1=L[0,:];p2=L[1,:];p3=L[2,:];p4=L[3,:]
     35 L1=np.mean(L[0,:]);L2=np.average(L[1,:]);L3=np.average(L[2,:]);L4=np.average(L[3,:])
     36 t1 = np.linspace(0,100,len(p1))

IndexError: index 3 is out of bounds for axis 0 with size 3

怎么了？

预期输出为 L1=.5162,L2=-.0001,L3=-4.9208,L4=-6.5954

Answer 1

在LEC_system(u)中，平面向量u包含序列

状态向量，n 个分量，
特征基础U，一个n x n矩阵
累积指数 L、n 个分量。

对于 n=4，这转化为分解

def LEC_system(u):
    #x,y,z,w = u[:4]
    U = u[4:20].reshape([4,4])
    L = u[20:24]
    f,Df = diff_Lorenz(u[:4])
    A = U.T.dot(Df.dot(U))
    dL = np.diag(A).copy();
    for i in range(4):
        A[i,i] = 0
        for j in range(i+1,4): A[i,j] = -A[j,i]
    dU = U.dot(A)
    return np.concatenate([f,dU.flatten(),dL])

当然，在集成后的评估中，同样必须使用正确的状态向量段

L = u[5:,20:24].T/t[5:]

然后我明白了

并且在集成到 t=60

之后仅使用后四分之一的图表

LES=  0.029214865425355396 -0.43816854013111833 -4.309199339754925 -6.28183676249535

这仍然不是预期值，因为它似乎沿着与解曲线横向的所有方向收缩。

4D 混沌系统 Lyapunov 指数

4D chaotic system Lyapunov exponent

python-3.x

matplotlib

ode

chaos