4D 混沌系统 Lyapunov 指数

Posted 2023-02-18

【中文标题】4D 混沌系统 Lyapunov 指数

【英文标题】：4D chaotic system Lyapunov exponent

【发布时间】：2022-01-03 22:33:04

【问题描述】：

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

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

【问题讨论】：

最后一个方程的导数应该是[y,x,0,-c]。还是等式有误？

我的错。但是还是不行

更新了我的问题

【参考方案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

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

【讨论】：

哇哦.. 请解释一下 U = u[4:20].reshape([4,4]) L = u[20:24]

u 向量在0, n, n+n*n, n+n*n+n 位置（之前）处拆分，外部位置不拆分。对于n=4 值0, 4, 20, 24 结果。

如果系统有 5D 那么... 0 , n, n+nn , n+nn+n, 那么呢？

是的，模式保持不变，然后在索引0, 5, 30, 35之前使用组件分离。

您是否从f 的第一个组件中删除了u？