0

我想通过使用本文第 81 页中描述的标准方法来数值计算洛伦兹系统的 Lyapunov 谱。

一个基本上需要整合洛伦兹系统和切向向量（我为此使用了龙格-库塔方法）。切向向量的演化方程由洛伦兹系统的雅可比矩阵给出。在每次迭代之后，需要对向量应用 Gram-Schmidt 方案并存储其长度。然后由存储长度的平均值给出三个 Lyapunov 指数。

我在 python 中实现了上面解释的方案（使用版本 3.7.4），但我没有得到正确的结果。

我认为错误在于 der 向量的 Rk4 方法，但我找不到任何错误......轨迹 x、y、z 的 RK4 方法工作正常（由图指示）和实施的 Gram-Schmidt 方案也正确执行。

我希望有人可以查看我的短代码并可能找到我的错误

编辑：更新的代码

from numpy import array, arange, zeros, dot, log
import matplotlib.pyplot as plt
from numpy.linalg import norm

# Evolution equation of tracjectories and tangential vectors
def f(r):
    x = r[0]
    y = r[1]
    z = r[2]

    fx = sigma * (y - x)
    fy = x * (rho - z) - y
    fz = x * y - beta * z

    return array([fx,fy,fz], float)

def jacobian(r):
    M = zeros([3,3])
    M[0,:] = [- sigma, sigma, 0]
    M[1,:] = [rho - r[2], -1, - r[0] ]
    M[2,:] = [r[1], r[0], -beta]

    return M

def g(d, r):
    dx = d[0]
    dy = d[1]
    dz = d[2]

    M = jacobian(r)

    dfx = dot(M, dx)
    dfy = dot(M, dy)
    dfz = dot(M, dz)

    return array([dfx, dfy, dfz], float)

# Initial conditions
d = array([[1,0,0], [0,1,0], [0,0,1]], float)
r = array([19.0, 20.0, 50.0], float)
sigma, rho, beta = 10, 45.92, 4.0 

T  = 10**5                         # time steps 
dt = 0.01                          # time increment
Teq = 10**4                        # Transient time

l1, l2, l3 = 0, 0, 0               # Lengths

xpoints, ypoints, zpoints  = [], [], []

# Transient
for t in range(Teq):
    # RK4 - Method 
    k1  = dt * f(r)                 
    k11 = dt * g(d, r)

    k2  = dt * f(r + 0.5 * k1)
    k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)

    k3  = dt * f(r + 0.5 * k2)
    k33 = dt * g(d + 0.5 * k22, r + 0.5 * k2)

    k4  = dt * f(r + k3)
    k44 = dt * g(d + k33, r + k3)

    r  += (k1  + 2 * k2  + 2 * k3  + k4)  / 6
    d  += (k11 + 2 * k22 + 2 * k33 + k44) / 6

    # Gram-Schmidt-Scheme
    orth_1 = d[0]                    
    d[0] = orth_1 / norm(orth_1)

    orth_2 = d[1] - dot(d[1], d[0]) * d[0]
    d[1] = orth_2 / norm(orth_2)

    orth_3 = d[2] - (dot(d[2], d[1]) * d[1]) - (dot(d[2], d[0]) * d[0]) 
    d[2] = orth_3 / norm(orth_3)

for t in range(T):
    k1  = dt * f(r)                 
    k11 = dt * g(d, r)

    k2  = dt * f(r + 0.5 * k1)
    k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)

    k3  = dt * f(r + 0.5 * k2)
    k33 = dt * g(d + 0.5 * k22, r + 0.5 * k2)

    k4  = dt * f(r + k3)
    k44 = dt * g(d + k33, r + k3)

    r  += (k1  + 2 * k2  + 2 * k3  + k4)  / 6
    d  += (k11 + 2 * k22 + 2 * k33 + k44) / 6

    orth_1 = d[0]                    # Gram-Schmidt-Scheme
    l1 += log(norm(orth_1))
    d[0] = orth_1 / norm(orth_1)

    orth_2 = d[1] - dot(d[1], d[0]) * d[0]
    l2 += log(norm(orth_2))
    d[1] = orth_2 / norm(orth_2)

    orth_3 = d[2] - (dot(d[2], d[1]) * d[1]) - (dot(d[2], d[0]) * d[0]) 
    l3 += log(norm(orth_3))
    d[2] = orth_3 / norm(orth_3)

# Correct Solution (2.16, 0.0, -32.4)

lya1 = l1 / (dt * T)
lya2 = l2 / (dt * T)  - lya1
lya3 = l3 / (dt * T) -  lya1 - lya2 

lya1, lya2, lya3

# my solution T = 10^5 : (1.3540301507934012, -0.0021967491623752448, -16.351653561383387)

上面的代码是根据 Lutz 的建议更新的。结果看起来好多了，但它们仍然不是 100% 准确的。

  • 正确解 (2.16, 0.0, -32.4)

  • 我的解决方案（1.3540301507934012，-0.0021967491623752448，-16.351653561383387）

正确的解决方案来自Wolf's Paper, p.289。在第 290-291 页，他描述了他的方法，该方法看起来与我在本文开头提到的论文中完全相同（论文，第 81 页）。

所以我的代码中一定还有另一个错误......

4

1 回答 1

1

您需要将点系统和雅可比系统求解为（前向）耦合系统。在原始源代码中正是这样做的，所有内容都在一次RK4调用组合系统中更新。

因此，例如在第二阶段，您将混合操作以具有组合的第二阶段

   k2 = dt * f(r + 0.5 * k1)   
   M = jacobian(r + 0.5 * k1)
   k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)

您还可以委托函数M内部的计算g，因为这是唯一需要它的地方，并且您增加了变量范围内的局部性。

请注意，我将dfrom的更新更改k1为k11，这应该是数值结果中错误的主要来源。

关于最后一个代码版本（2/28/2021）的补充说明：

正如评论中所说，代码看起来像算法的数学规定的那样。有两个误读会阻止代码返回接近引用的结果：

  • 论文中的参数是sigma=16。
  • 本文使用的不是自然对数，而是二进制对数，即幅度演化为2^(L_it)。所以你必须将计算出的指数除以log(2)。

使用https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent中派生的方法，我得到了指数

[2.1531855610566595, -0.00847304754613621, -32.441308372177566]

这是足够接近参考(2.16, 0.0, -32.4)。

于 2020-02-21T12:43:24.637 回答