3

目前我使用以下代码生成洛伦兹级数

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)

但nolds.lyap_r产生无效值（我假设有效值为 0.91）

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

知道我在哪里犯了错误吗？

4

1 回答 1

0

错误在于您假设 Lorenz 动力学的 x 坐标对应于第一个 Lyapunov 指数。观察你正在服用：

ret.append(xtemp[0])

然而，第一个 Lyapunov 指数量化了不稳定流形的更不稳定方向上的发散率。

如我所见，您只是在估计 x 坐标的第一个 Lyapunov 指数。此外，在这种方法中，对于每个坐标 {x,y,z}，Lyapunov 指数将为正，因为这种“平凡”的分解不会捕获稳定流形。然后，您永远不会以这种方式找到第三个李雅普诺夫指数（负）。

解决方案是使用 Gram-Schmidt 过程来获得动力学膨胀和收缩的正确方向，从而计算所有 Lyapunov 指数。最大值正是您正在寻找的 Lyapunov 指数（大约 0.9）。尽管如此，有些论文对定性结果 (+,0,-) 比对量级更感兴趣，因此，也许您可​​以找到一些其他论文显示最大 Lyapunov 指数的值略有不同。

值得注意的是，如果我们考虑每个变量的总和，以构建一个新信号，则与这个新信号相关的最大李雅普诺夫指数达到您所期望的值。我绘制了从 8000 点到 10000 点的信号，并获得了附在这篇文章中的图。

于 2019-10-16T04:42:43.433 回答