computing Lyapunov exponents numerically

Question

I am trying to compute numerically the Lyapunov exponents of an ODE. I follow the method described in Parker, Chua "Practical Numerical Algorithms for Chaotic systems"

There is also relevant answer here Numerical computation of Lyapunov exponent

Below is my code. It gives wrong values for the Lorenz system at classical values. There must be something wrong with it, but I do not see what. Please help me to find the mistake.

#!/usr/bin/python3
import numpy as np
import math
sigma = 10.
rho = 28.0
beta = 8./3.
def lor(x):
    y = []
    y.append(sigma(x[1]-x[0]))
    y.append(rhox[0]-x[1]-x[0]x[2])
    y.append(x[0]x[1]-beta*x[2])
    return np.array(y)
def lorj(x):
    z = []
    y1 =[]
    y1.append(-sigma)
    y1.append(sigma)
    y1.append(0.)
    y2 = []
    y2.append(rho - x[2])
    y2.append(-1.)
    y2.append(-x[0])
    y3 = []
    y3.append(x[1])
    y3.append(x[0])
    y3.append(-beta)
    z.append(y1)
    z.append(y2)
    z.append(y3)
    return np.array(z)
def osel(y,h,f,j,u):
    k1 = f(y)
    k2 = f(y+0.5hk1)
    k3 = f(y+0.5hk2)
    k4 = f(y+hk3)
    h6 = h/6.0
    for i in range(len(y)):
        k1[i]= h6
        k2[i]= 2.h6
        k3[i]= 2.h6
        k4[i]= h6
    J = j(y)
    v = np.dot(J,u)
    u = u+hv
    return y+k1+k2+k3+k4,u,v
def classic(y,h,f,j,u,T):
    for t in range(T):
        y,u,v = osel(y,h,f,j,u)
    return y,u
def log(x):
    y = []
    for i in x:
        if i>0:
            y.append(math.log(i))
    return np.array(y)
def dot(x,y):
    r = 0.
    for i in range(len(x)):
        r+=x[i]*y[i]
    return r
def orth(x):
    u = []
    v = []
    for i in range(len(x)):
        y = [k for k in x[i]]
        for j in range(i):
            c = dot(x[i],u[j])
            for k in range(len(y)):
                y[k] -= cu[j][k]
        v.append([k for k in y])
        a = 1./math.sqrt(dot(y,y))
        for k in range(len(y)):
            y[k] = a
        u.append([k for k in y])
    return v,u
h= 0.0001
y = [3.,1.,1.]
u = [[1,0,0],[0,1,0],[0,0,1]]
L = [0. for i in range(len(y))]
K0 = 1000
for K in range(K0):
    y,u, = classic(y,h,lor,lorj,u,int(1./h))
    v,u = orth(u)
    s= ""
    for i in range(len(y)):
        L[i] += math.log(math.sqrt(dot(v[i],v[i])))
        s+= " "+str(L[i]/(K+1))
    print(K,s)

EDIT:---------------------

I agree with the criticism in the comments. I improved the question as follows.

The algorithm in Parker, Chua pp. 73-81, in particular formula (3.39). This is so-called Benettin algorithm. This is the formula I implemented.

Let me summarize.

Together with the system \begin{equation} dx_i = f_i(x)dt \end{equation}

we solve the system \begin{equation} du = J(x(t))u dt \end{equation} where $J(x)$ is the Jacobian $\partial_i f_j$, $ u(t) $ is a matrix, initialized by an orthonormal system of vectors.

Then the computation consists in iteration of the following steps

1. evolve the system $(x,u)$ for time $T$. Here $u = [U_1,...,U_n]$, $ U_i$ are n-dimensional vectors, initialized at $(U_i)_j =\delta_{ij}$
2. then we orthonormalize the system

\begin{gather} v_1 = u_1 \\ \nonumber w_1 = v_1/|v_1| \\ \nonumber v_2 = u_2 - <u_2,w_1>w_1 \\ \nonumber w_2 = v_2/|v_2| \\ \nonumber .... \\ \nonumber v_i = u_i - <u_i,w_1>w_1-....-<u_i,w_{i-1}>w_{i-1} \\ \nonumber w_i = v_i/|v_i|\\ ... \end{gather}

We use $[w_1,...,w_n]$ as initial conditions for $u$ at the next step.

The Lyapunov spectrum is computed as \begin{equation} \lambda_i = \frac{1}{KT} \sum_{k=1}^K \ln|v^{(k)}_i| \end{equation} where $v^{(k)}_i$ is the vector $v_i$ computed at $k$-th iteration.

My original implementation of the algorithm has a mistake. I was orthonormalizing the rows of matrix $u$, while I should orthonormalize columns, as this is what is used in the evolution equations for the nearby trajectories.

I rewrote the code, as suggested in the comments. I attach the code in c below. I also wrote a python code that uses numpy arrays. But this code works horribly slow. Please suggest how to improve it so that it has a better performance.

I also tested statements made in the comments that RK4 can be trusted as low as h=0.03, see attached plots for the difference between trajectories for the time interval [0:5]. This seems to be a correct statement. The next statement is that we can trust the simulatiton until time ~5. This also seems to be true, see the second plot. After T~10, the trajectories initially separated by 1e-5 start to have distance >1.

The code as it stands now seems to be in qualitative agreement with the literature. But still, I am observing ~10% discrepancy, which is very unsatisfactory.

As a benchmark, I use the following paper

Marek Balcerzak · Danylo Pikunov · Artur Dabrowski

"The fastest, simplified method of Lyapunov exponents spectrum estimation for continuous-time dynamical systems"

The relevant plot from this paper is attached below.

I obtain similar looking plots, but there is some "noise" on them. What is the origin of this "noise"?

Question: How to improve the code so that it reproduces the calculations in the paper?

#include <stdio.h>
#include <math.h>
float sigma = 10.0;
float rho =30;//2.5;//28.0;
float beta = 8.0/3.0;
/*
float sigma = 16.0;
float rho = 45.92;
float beta = 4.;
*/
void lor(float f,float x){
    f[0] = sigma(x[1]-x[0]);
    f[1] = rhox[0] - x[1] - x[0]x[2];
    f[2] = x[0]x[1] - beta*x[2];
}
void lorj(float* J,float x){
    J[13+0] = rho - x[2];
    J[1*3+2] = -x[0];
    J[2*3+0] = x[1];
    J[2*3+1] = x[0];
}
void mul_m(floatz,floatx, float* y,int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            float u=0;
            for(int k=0;k<n;k++){
                u+=x[i*n+k]y[kn+j];
            }
            z[i*n+j]=u;
        }
    }
}
void mul_v(floatz,float h,floatx,int n){
    for(int i=0;i<n;i++){
        z[i] = h*x[i];
    }
}
void add_v(floatz,floatx,float*y,int n){
    for(int i=0;i<n;i++){
        z[i] = x[i]+y[i];
    }
}
void mul_ms(floatz,float h,floatx,int n){
    for(int i=0;i<nn;i++){
        z[i]=hx[i];
    }
}
void add_m(floatz,floatx,floaty,int n){
    for(int i=0;i<nn;i++){
        z[i]=x[i]+y[i];
    }
}
void copy(floaty,floatx,int n){
    for(int i=0;i<n;i++){
        y[i]=x[i];
    }
}
void print_v(float*x,int n){
    for(int i=0;i<n;i++){
        printf("%f ",x[i]);
    }
    printf("\n");
}
void rk4(int n,float* x,float* y,float* tmp,float* k1,float* k2,float* k3,float* k4,float* a1,float* a2,float* a3,float w,floatv,float h,void f(float,float),void J(float,float),float u){
    f(k1,x);
    mul_v(tmp,0.5h,k1,n);
    add_v(a1,x,tmp,n);
    f(k2,a1);
    mul_v(tmp,0.5h,k2,n);
    add_v(a2,x,tmp,n);
    f(k3,a2);
    mul_v(tmp,h,k3,n);
    add_v(a3,x,tmp,n);
    f(k4,a3);
    float h6 = h/6;
    for(int i=0;i<n;i++){
        y[i] = x[i] + h6k1[i]+2h6k2[i]+2h6k3[i]+h6*k4[i];
    }
    J(w,x);
    mul_m(v,w,u,n);
    mul_ms(v,h,v,n);
    add_m(u,u,v,n);
}
void orth(floatu,floatv,floatx,int n){
    for(int i=0;i<n;i++){
        for(int k=0;k<n;k++){
            u[kn+i]=x[k*n+i];
        }
        for(int j=0;j<i;j++){
            float c =0;
            for(int k=0;k<n;k++){
                c+=x[k*n+i]v[kn+j];
            }
            for(int k=0;k<n;k++){
                u[k*n+i]-=cv[kn+j];
            }
        }
        float c=0;
        for(int k=0;k<n;k++){
            c+=u[k*n+i]u[kn+i];
        }
        c = sqrt(c);
        if(c==0){
            printf("c=0\n");
        }else{
            for(int k=0;k<n;k++){
                v[k*n+i]=u[k*n+i]/c;
            }
        }
    }
}
void evol(int n,float* x,float* y,float* tmp,float* k1,float* k2,float* k3,float* k4,float* a1,float* a2,float* a3,float w,floatv,float h,void f(float,float),void J(float,float),float *u,int T){
    for(int t=0;t<T;t++){
        rk4(n,x,y,tmp,k1,k2,k3,k4,a1,a2,a3,w,v,h,f,J,u);
        copy(x,y,n);
    }
}
void benettin(int n,float* x,float* y,float* tmp,float* k1,float* k2,float* k3,float* k4,float* a1,float* a2,float* a3,float w,floatv,float h,void f(float,float),void J(float,float),float u,int T,int Q,floatr,floats,floatl,floatl1){
    for(int q=0;q<Q;q++){
        evol(n,x,y,tmp,k1,k2,k3,k4,a1,a2,a3,w,v,h,f,J,u,T);
        orth(r,s,u,n);
        for(int i=0;i<n;i++){
            float z=0;
            for(int j=0;j<n;j++){
                z+=r[jn+i]r[jn+i];
            }
            if(z==0){
                printf("z==0\n");
            }else{
                l[i] += 0.5*log(z);
            }
        }
    copy(u,s,n*n);
    mul_v(l1,1./(float(1+q)*h*float(T)),l,n);
    //printf(&quot;q=%i  &quot;,q);
    //print_v(l1,n);
}    

}
void call(floatp,float h, int T,int Q){
    int n=3;
    float x = new floatn;
    float* y = new floatn;
    float* tmp = new floatn;
    float* k1 = new floatn;
    float* k2 = new floatn;
    float* k3 = new floatn;
    float* k4 = new floatn;
    float* a1 = new floatn;
    float* a2 = new floatn;
    float* a3 = new floatn;
    float* u = new floatn*n;
    float* v = new floatn*n;
    float* w = new floatn*n;
    float* r = new floatn*n;
    float* s = new floatn*n;
    float* l = new floatn;
    float* l1 = new floatn;
w[0*3+0] = -sigma;
w[0*3+1] = sigma;
w[0*3+2] = 0;
w[1*3+1] = -1.0;
w[2*3+2] = -beta;

x[0] = 1;
x[1] = 1;
x[2] = 1;

for(int i=0;i&lt;n;i++){
    for(int j=0;j&lt;n;j++){
        u[i*n+j] =0;
    }
}
for(int i=0;i&lt;n;i++){
    u[i*n+i] = 1;
}

benettin(n,x,y,tmp,k1,k2,k3,k4,a1,a2,a3,w,v,h,lor,lorj,u,T,Q,r,s,l,l1);

copy(p,l1,n);

delete x;
delete y;
delete tmp;
delete k1;
delete k2;
delete k3;
delete k4;
delete a1;
delete a2;
delete a3;
delete u;
delete v;
delete w;
delete r;
delete s;
delete l;
delete l1;

}
int main(){
    float dr = 0.001;
    float* p = new float3;
    FILE*F;
    F = fopen("r-dependence_dr=0.001_h=0.03_Q=10.txt","w");
for(int r=0;r&lt;100./dr;r++){
    rho = dr*float(r);
    call(p,0.03,33,10);
    fprintf(F,&quot;%f %f %f %f\n&quot;,rho,p[0],p[1],p[2]);

    printf(&quot;-----&gt; %f %f %f %f\n&quot;,rho,p[0],p[1],p[2]);
}
fclose(F);
delete p;

}

import numpy as np
import math


sigma = 10.
rho = 28.0
beta = 8./3.


def lor(x):
    return np.array([sigma*(x[1]-x[0]),rho*x[0]-x[1]-x[0]*x[2],x[0]*x[1]-beta*x[2]])

def lorj(x):
    return np.matrix([[-sigma,sigma,0.],[rho - x[2],-1.,-x[0]],[x[1],x[0],-beta]])

def rk4(x,u,h,f,jac):
    k1 = f(x)
    k2 = f(x+0.5*h*k1)
    k3 = f(x+0.5*h*k2)
    k4 = f(x+h*k3)
    y = x+(h/6.)*(k1+2.*k2+2.*k3+k4)
    J = jac(x)
    v = u + h*J@u
    return y,v



def evol(x,u,h,f,jac,T):
    for t in range(T):
        x,u = rk4(x,u,h,f,jac)
    return x,u

def orth(x):
    n = len(x)
    u = np.matrix([[0. for i in range(n)] for j in range(n)])
    v = np.matrix([[0. for i in range(n)] for j in range(n)])
    for i in range(n):
        for k in range(n):
            u[k,i] = x[k,i]
        for j in range(i):
            c = 0. 
            for k in range(n):
                c += x[k,i]*v[k,j]
            for k in range(n):
                u[k,i] -= c*v[k,j]
        c =0. 
        for k in range(n):
            c += u[k,i]*u[k,i]
        c = math.sqrt(c)
        for k in range(n):
            v[k,i] = u[k,i]/c
    return u,v         


def benettin(x,u,h,f,jac,T,Q):
    n = len(x)
    L = np.array([0. for k in range(len(x))])
    for q in range(Q):
        x,u = evol(x,u,h,f,jac,T)        
        u,v = orth(u)
        for i in range(n):
            c = 0.
            for k in range(n):
                c += u[k,i]*u[k,i]
            c = math.sqrt(c)
            c = math.log(c)
            L[i] += c

        u = v
        print(q,(1./((1.+q)*h*T))*L)



x=np.array([1.,1.,1.])
u=np.matrix([[1,0,0],[0,1,0],[0,0,1]])


benettin(x,u,0.0001,lor,lorj,10000,1000)

---

"Dependence on RK step size"

"Separation between nearby trajectories as function of time"

"Lyapunov exponents from the paper"

"Lyapunov exponents for total time 100"

"Lyapunov's for total time 300"

Answer 1

You do not reach any kind of convergence in running this code, as your total time span is T=K0*h=1e+3*1e-4=0.1, which is much too short. You will need at least T=5 to reliably reach the attractor. You will need T>50 to start to even out the oscillations of the local exponents. The local exponents are a quasi-periodic oscillation, it is only their time-averages that converge with speed 1/t, which is very slow.

Runge-Kutta for the Lorenz system works for step sizes down from h=0.05, with step size h=0.0001 you are in the region of optimal accuracy for 64bit floats. For Lyapunov exponents you do not need optimal accuracy, only a guarantee that the solution remains qualitatively valid. Such a step size could be h=0.005 to h=0.01.