在 Python 中获取特征李雅普诺夫指数的意外值 - 我的贝纳通算法实现可能出什么问题？

问：

我写了一个代码来计算 Python 中 Lorenz 系统的特征 Lyapunov 指数。我肯定做错了什么，因为CLE还很远，但我无法弄清楚是什么。 我正在使用 Benettin 等人的算法。 我首先随机选择初始状态，然后使用 Gram Schmidt 过程（也计算体积）正交 3 个随机切向量。 初始化后，我开始迭代。我使用欧拉方法演化参考轨迹，并使用切向量演化，其中 L 是 x（t） 处的雅可比矩阵，w（t） 是切向量。我正交每一个并计算对数（卷）的总和。每个CLE都被打印出来。 如果有人知道可能是什么问题，我将不胜感激。 我正在为 Lorenz 系统实现此算法以验证它，以便将其用于不同的系统。w(t+1)=L@w(t)northoiwrite

代码如下：

import numpy as np

def lorenz(x, y, z, sigma, rho, beta):
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

def lorenz_L(x,y,z):
    return np.array([
    [-sigma, sigma, 0],
    [rho - z, -1, -x],
    [y, x, -beta]])   

def gram_schmidt(V):
    """
    Orthogonalize a set of vectors stored as the ROWS of matrix V.
    It also computes the volume generated the first i vectors for i=1,...,dim
    """
    # Get the number of vectors.
    n = V.shape[0]
    W = np.copy(V)
    R = np.zeros(n)
    VOL = np.zeros(n)
    for j in range(n):
        # To orthogonalize the vector in row j with respect to the
        # previous vectors, subtract from it its projection onto
        # each of the previous vectors.
        for k in range(j):
            W[j] -= np.dot(W[k], W[j]) * W[k]
        R[j] = np.linalg.norm(W[j])
        W[j] = W[j] / R[j]
        VOL[0]=R[0]
    for k in range(1,n):
        VOL[k]=VOL[k-1]*R[k]
    return W, VOL

# Parameters
sigma = 16.0
rho = 45.92
beta = 4.0

t_start = 0.0
t_end = 100.0
steps = 10000
northo = 1
iwrite = 1000 
num_points = steps * northo
dt = t_end/num_points
icont=0
dim=3

X = np.random.random(dim)

# Estimate Lyapunov exponents

num_lyapunov_exponents = dim
lyapunov_exponents = np.zeros(num_lyapunov_exponents)

num_vec=dim

sum_volumes = np.zeros(dim)
sum_lyapunov_exponents = np.zeros(dim)
cle = np.zeros(dim)

tangent_vec = np.random.random((dim,dim))
tangent_vec, _ = gram_schmidt(tangent_vec) 

for j in range(steps):
    for is_ in range(northo):
        icont=icont+1
        # velocities
        dX_dt = lorenz(X[0],X[1],X[2],sigma,rho,beta)

        # evolution of the reference trajectory
        for i in range(dim):
            X[i] = X[i] + dX_dt[i] * dt
 
        # stability matrix for Lorenz system
        L = lorenz_L(X[0],X[1],X[2])
        
        # evolve tangent vectors
        tangent_vec = (L@tangent_vec.T).T
        # for ivec in range(num_vec):
        #     tangent_vec[ivec] = L@tangent_vec[ivec]
    
    #orthonormalize
    tangent_vec, volumes = gram_schmidt(tangent_vec) 

    for i in range(dim):
        sum_volumes[i]=sum_volumes[i]+np.log(volumes[i])
    
    # calculate and print results
    if icont%iwrite == 0:
        t = icont*dt
        print(f'time : {t}, dt : {dt}, frame : {icont}')
        sum_lyapunov_exponents = sum_volumes / t
        cle[0]=sum_lyapunov_exponents[0]
        for l in range(1,dim):
            cle[l]=sum_lyapunov_exponents[l]-sum_lyapunov_exponents[l-1]
        for ii in range(dim):
            print(f'{ii+1}   {cle[ii]}')
        print()

我尝试了上面的代码，根据我设置的时间步长值，我得到了 Lyapunov 指数。我希望得到 Kaplan York 维数为 2.07 的理论期望，第一 lyapunov 指数为正 （1.5;2）、第二个零和第三个负（-20;-30）。