0

我正在尝试研究 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

4

1 回答 1

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

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

于 2021-11-25T19:16:04.600 回答