# computing Lyapunov exponents numerically

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(rho*x[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.5*h*k1)
k3 = f(y+0.5*h*k2)
k4 = f(y+h*k3)
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+h*v
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] -= c*u[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 $$$$dx_i = f_i(x)dt$$$$

we solve the system $$$$du = J(x(t))u dt$$$$ 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 - w_1 \\ \nonumber w_2 = v_2/|v_2| \\ \nonumber .... \\ \nonumber v_i = u_i - w_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 $$$$\lambda_i = \frac{1}{KT} \sum_{k=1}^K \ln|v^{(k)}_i|$$$$ 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] = rho*x[0] - x[1] - x[0]*x[2];
f[2] = x[0]*x[1] - beta*x[2];
}

void lorj(float* J,float *x){
J[1*3+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(float*z,float*x, 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[k*n+j];
}
z[i*n+j]=u;
}
}
}

void mul_v(float*z,float h,float*x,int n){
for(int i=0;i<n;i++){
z[i] = h*x[i];
}
}

for(int i=0;i<n;i++){
z[i] = x[i]+y[i];
}
}

void mul_ms(float*z,float h,float*x,int n){
for(int i=0;i<n*n;i++){
z[i]=h*x[i];
}
}

for(int i=0;i<n*n;i++){
z[i]=x[i]+y[i];
}
}

void copy(float*y,float*x,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,float*v,float h,void f(float*,float*),void J(float*,float*),float *u){
f(k1,x);
mul_v(tmp,0.5*h,k1,n);
f(k2,a1);
mul_v(tmp,0.5*h,k2,n);
f(k3,a2);
mul_v(tmp,h,k3,n);
f(k4,a3);
float h6 = h/6;
for(int i=0;i<n;i++){
y[i] = x[i] + h6*k1[i]+2*h6*k2[i]+2*h6*k3[i]+h6*k4[i];
}
J(w,x);
mul_m(v,w,u,n);
mul_ms(v,h,v,n);
}

void orth(float*u,float*v,float*x,int n){
for(int i=0;i<n;i++){
for(int k=0;k<n;k++){
u[k*n+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[k*n+j];
}
for(int k=0;k<n;k++){
u[k*n+i]-=c*v[k*n+j];
}
}
float c=0;
for(int k=0;k<n;k++){
c+=u[k*n+i]*u[k*n+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,float*v,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,float*v,float h,void f(float*,float*),void J(float*,float*),float *u,int T,int Q,float*r,float*s,float*l,float*l1){
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[j*n+i]*r[j*n+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("q=%i  ",q);
//print_v(l1,n);
}
}

void call(float*p,float h, int T,int Q){
int n=3;
float* x = new float[n]();
float* y = new float[n]();
float* tmp = new float[n]();
float* k1 = new float[n]();
float* k2 = new float[n]();
float* k3 = new float[n]();
float* k4 = new float[n]();
float* a1 = new float[n]();
float* a2 = new float[n]();
float* a3 = new float[n]();
float* u = new float[n*n]();
float* v = new float[n*n]();
float* w = new float[n*n]();
float* r = new float[n*n]();
float* s = new float[n*n]();
float* l = new float[n]();
float* l1 = new float[n]();

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<n;i++){
for(int j=0;j<n;j++){
u[i*n+j] =0;
}
}
for(int i=0;i<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 float[3]();
FILE*F;
F = fopen("r-dependence_dr=0.001_h=0.03_Q=10.txt","w");

for(int r=0;r<100./dr;r++){
rho = dr*float(r);
call(p,0.03,33,10);
fprintf(F,"%f %f %f %f\n",rho,p[0],p[1],p[2]);

printf("-----> %f %f %f %f\n",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"
• Lots of code, no comments, no explanation of what the algorithm it implements should be, no discussion of what you expect the answer to be, nor what you actually get, no discussion of simple test cases you have already tried => Little chance of a good answer. Jul 13, 2023 at 4:41
• Welcome to Scicomp! If I were you, i would take a step back and calculate the exponent for a simple 1D toy problem first. Have you checked: scicomp.stackexchange.com/questions/36013 ? Jul 13, 2023 at 6:10
• @MPIchael: While you can define and calculate Lyapunov exponents for lower-dimensional systems, many problems only arise for oscillatory and chaotic systems and for those you need at least two or three dimensions, respectively. Also, in 1D many potential pitfalls vanish. While I agree with the other criticism on this question, starting out in 3D is okay. Jul 13, 2023 at 6:35
• If you want to code in C++, then code in C++. If you use Python with numpy, then really use numpy, such as using np.dot, np.dot, the .copy method of arrays, the scalar multiplication of arrays, etc. This is just so inconsistent, for instance in osel in the slope computation you use the scalar vector multiplication, but next you make a loop to multiply element wise instead of just using dx=(k1+2*k2+2*k3+k4)*h/6. The shorter code would then make the question more readable. Jul 13, 2023 at 7:01
• Have you looked into using an integrator from something like SciPy with adaptive time stepping so you don't have to worry about the dependence on step size and you can just set the error tolerances yourself? For long-time integration, your method with a fixed step size may be drifting your solution and causing large errors in asymptotic behavior. It also removes one more source of user error Jul 17, 2023 at 21:24

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.