Python Nolds: how to get proper value for Lorenz system

currently i use following code to generate lorenz series

def generate(x, stop=10000, s=10, b=8/3, r=28):
    def lor(v):
        return np.array([s * (v[1] - v[0]), v[0] * (r - v[2]) - v[1], v[0] * v[1] - b * v[2]])
    ret = []
    step = 0.1
    xtemp = x.copy()
    for i in range(stop):
        k1 = lor(xtemp)
        k2 = lor(xtemp + step / 2 * k1)
        k3 = lor(xtemp + step /2 * k2)
        k4 = lor(xtemp + step * k3)
        xtemp += step/6 * (k1 + 2 * k2 + 2 * k3 + k4)
        ret.append(xtemp[0])
    return np.array(ret)

but nolds.lyap_r yields invalid value (i assume that valid is 0.91)

import nolds l = generate([1, 0, 0]) nolds.lyap_r(l, tau=0.1, emb_dim=5) 1.0030932070169234

any idea where did i made a mistake?

Upvotes: 3

Views: 930

Answers (1)

The mistake is that you are assuming that the x-coordinate of the Lorenz dynamics corresponds to the first Lyapunov exponent. Observe that you are taking:

ret.append(xtemp[0])

However, the first Lyapunov exponent quantifies the rate of divergence in the more unstable direction of the unstable manifold.

As I can see, you are only estimating the first Lyapunov exponent of the x-coordinate. Moreover, in this approach, for each coordinate {x,y,z}, the Lyapunov exponent will be positive, because this "trivial" decomposition does not capture the stable manifold. Then, you never will find the 3rd Lyapunov exponent (negative) in this way.

The solution is to use the Gram-Schmidt process to get the proper directions of expansion and contraction of your dynamics and therefore calculate all the Lyapunov exponents. The maximum is precisely the Lyapunov exponent that you are looking (approx 0.9). Nevertheless, some papers are more interested in the qualitative result (+,0,-) than the magnitude, so, maybe you can find some other papers showing slightly different values for the maximum Lyapunov exponent.

It is noteworthy that, if we consider the sum of each variable, to build a new signal, the Largest Lyapunov Exponent associated with this new signal, is reached the value that you are expecting. I plot for signal from 8000 points to 10000 points, and obtain the plot attached to this post.

Upvotes: 0

Related Questions