已知 ODE 的李亚普诺夫谱 - Python 3

Question

我想用数值计算 Lorenz System by using the standard method which is described in this Paper, p.81 的李雅普诺夫谱。

基本上需要对洛伦兹系统和切向向量进行积分（为此我使用了 Runge-Kutta 方法）。切向矢量的演化方程由洛伦兹系统的雅可比矩阵给出。每次迭代后，需要对向量应用 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 页，他描述了他的方法，该方法看起来与我在本文开头提到的论文完全相同 post（论文，第 81 页）。

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

Answer 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 函数，因为这是唯一需要它的地方，并且您增加了变量范围的局部性。

注意我把d的update由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).

足够接近

已知 ODE 的李亚普诺夫谱 - Python 3

Lyapunov Spectrum for known ODEs - Python 3

python

ode

chaos