Lyapunov Spectrum for known ODEs - Python 3

Question

I want to numerically compute the Lyapunov Spectrum of the Lorenz System by using the standard method which is described in this Paper, p.81.

One basically need integrate the Lorenz system and the tangential vectors (i used the Runge-Kutta method for this). The evolution equation of the tangential vectors are given by the Jacobi matrix of the Lorenz system. After each iterations one needs to apply the Gram-Schmidt scheme on the vectors and store its lengths. The three Lyapunov exponents are then given by the averages of the stored lengths.

I implemented the above explained scheme in python (used version 3.7.4), but I don't get the correct results.

I thing the bug lies in the Rk4-Method for der vectors, but i cannot find any error...The RK4-method for the trajectories x,y,z works correctly (indicated by the plot) and the implemented Gram-Schmidt scheme is also correctly implemented.

I hope that someone could look through my short code and maybe find my error

---

Edit: Updated Code

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) 

The above code is updated according to Lutz suggestions. The results look much better but they are still not 100% accurate.

• Correct Solution (2.16, 0.0, -32.4)

• My solution (1.3540301507934012, -0.0021967491623752448, -16.351653561383387)

The correct solutions are from Wolf's Paper, p.289. On page 290-291 he describes his method, which looks exactly the same as in the paper that i mentioned in the beginning of this post (Paper, p.81).

So there must be another error in my code...

Answer 1

You need to solve the system of point and Jacobian as the (forward) coupled system that it is. In the original source exactly that is done, everything is updated in one RK4 call for the combined system.

So for instance in the second stage, you would mix the operations to have a combined second stage

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

You could also delegate the computation of M inside the g function, as this is the only place where it is needed, and you increase locality in the scope of variables.

Note that I changed the update of d from k1 to k11, which should be the main source of the error in the numerical result.

---

Additional remarks on the last code version (2/28/2021):

As said in the comments, the code looks like it will do what the mathematics of the algorithm prescribes. There are two misreadings that prevent the code from returning a result close to the reference:

• The parameter in the paper is sigma=16.
• The paper uses not the natural logarithm, but the binary one, that is, the magnitude evolution is given as 2^(L_it). So you have to divide the computed exponents by log(2).

Using the method derived in https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent I get the exponents

[2.1531855610566595, -0.00847304754613621, -32.441308372177566]

which is sufficiently close to the reference (2.16, 0.0, -32.4).