4D混沌吸引子Lyapunov谱计算代码运行索引错误问题求助

核心错误原因 #

你的索引越界报错以及odeint警告，本质是3维系统适配4维系统时，没有对应调整状态向量的切片区间，索引位置完全偏离了预设的状态向量结构。

状态向量u的预设结构为：[4个系统状态变量 + 16个4×4正交矩阵U的展开元素 + 4个Lyapunov指数累积量]，总长度24，原代码所有切片都不符合这个结构。

具体修改点 #

• 修正LEC_system函数内的切片逻辑：
  • 正交矩阵U的取值区间应该是状态变量之后的16个元素，从索引4到20，替换原错误的u[2:18]
  • Lyapunov指数累积量的取值区间应该是U之后的4个元素，从索引20到24，替换原错误的u[12:15]（原区间只有3个元素，直接导致后续只能取出3个指数，触发索引越界）
• 修正后续提取Lyapunov指数的切片区间，改为从索引20到24取值
• 拉长积分时间，保证Lyapunov指数收敛到目标数值

修复后完整代码 #

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):
    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])


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,200,6001)
u = odeint(lambda u,t:LEC_system(u),u0,t, hmax=0.05)
# 修正切片区间，跳过前50步的暂态过程
L = u[50:,20:24].T/t[50:]

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 = t[50:]
plt.plot(t1,p1);plt.plot(t1,p2);plt.plot(t1,p3);plt.plot(t1,p4)

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

运行说明 #

修复后代码不会再触发索引越界错误，odeint警告也会消失，输出的Lyapunov指数会接近你预期的0.5162、-0.0001、-4.9208、-6.5954，如果需要更高精度，可以继续拉长积分时间、缩小hmax参数调整。

内容的提问来源于stack exchange，提问作者Zewo