已知ODE的Lyapunov谱-Python 3

Question

我想通过使用此Lorenz System中描述的标准方法来数值计算Paper, p.81的Lyapunov谱。

基本上需要集成Lorenz系统和切向矢量（为此我使用了Runge-Kutta方法）。切向向量的演化方程由Lorenz系统的Jacobi矩阵给出。每次迭代后，需要对向量应用Gram-Schmidt方案并存储其长度。然后，通过存储长度的平均值给出三个Lyapunov指数。

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

我的错误在于der向量的Rk4-方法，但我找不到任何错误...轨迹x，y，z的RK4-方法正确运行（由图表示），并且已实现Gram-施密特方案也得到了正确实施。

我希望有人可以浏览我的短代码并找到我的错误

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

# Lorenz-System: x' = sigma*(y-x) 
#                y' = x*(rho - z) - y
#                z' = x*y - beta*z

# 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 g(d, r):
    dx = d[0]
    dy = d[1]
    dz = d[2]

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

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

# Jacobi Matrix of Lorenz System
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

# Initial conditions
d = array([[1,0,0], [0,1,0], [0,0,1]], float)
r = array([10.0, 1.0, 0.0], float)
sigma, rho, beta = 10, 45.92, 4.0 

T  = 10**5                         # time steps and increment
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):
    k1 = dt * f(r)                #RK4 for trajectories x,y,z
    k2 = dt * f(r + 0.5 * k1)   
    k3 = dt * f(r + 0.5 * k2)
    k4 = dt * f(r + k3)
    r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

    M = jacobian(r)
    k11 = dt * g(d, r)            #RK4 for tangent vectors
    k22 = dt * g(d + 0.5 * k1, r)
    k33 = dt * g(d + 0.5 * k2, r)
    k44 = dt * g(d + k3, r)
    d  += (k11 + 2 * k22 + 2 * k33 + k44) / 6

    xpoints.append(r[0]), ypoints.append(r[1]), zpoints.append(r[2])

    orth_1 = d[0]                 # Gram-Schmidt-Scheme
    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)            
    k2 = dt * f(r + 0.5 * k1)
    k3 = dt * f(r + 0.5 * k2)
    k4 = dt * f(r + k3)
    r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

    k11 = dt * g(d, r)        
    k22 = dt * g(d + 0.5 * k1, r)
    k33 = dt * g(d + 0.5 * k2, r)
    k44 = dt * g(d + k3, r)
    d  += (k11 + 2 * k22 + 2 * k33 + k44) / 6

    orth_1 = d[0]             
    d[0] = orth_1 / norm(orth_1)
    l1 += log(norm(orth_1))              # add lengths up

    orth_2 = d[1] - dot(d[1], d[0]) * d[0]
    d[1] = orth_2 / norm(orth_2)
    l2 += log(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)
    l3 += log(norm(orth_3))    

plt.plot(xpoints, ypoints)
plt.show()

# Correct Solution (2.16, 0.0, -32.4)
lya1 = l1 / (dt * T)
lya2 = l2 / (dt * T) 
lya3 = l3 / (dt * T)

lya1, lya2, lya3

#My result (-3.286149130588001e-15, 12.182263788663551, -9.503287701985853)

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的更新从k1更改为k11，这应该是数值结果中错误的主要来源。