# Numerical computation of Lyapunov exponents: how to find convergence or non-convergence efficiently?

I am wondering what are the standards for convergence of Lyapunov exponents (and Kaplan-Yorke dimension)? For example, I have a MATLAB code to calculate Lyapunov exponents for the classic Lorenz attractor using QR decomposition:

clear all; close all;

% parameters
par = [10; 28; 8/3];

% initial conditions
x0 = [1; 1; 1];

% time
dt = 0.001;
tspan = dt:dt:50;

%% time integration of the system
options = odeset('RelTol', 1e-12, 'AbsTol',1e-12*ones(1,3));
[t,x]=ode45(@(t,x) lorenz(t,x,par), tspan, x0,options);

plot3(x(:,1),x(:,2),x(:,3),'b','LineWidth',1.5);
xlabel('x')
ylabel('y')
zlabel('z')
set(gca,'FontSize',15)

%% find lyapunov exponent with QR decompostion
% initial deformation matrix, I
M = eye(3);
% empty vector to store Lyapunov exponents
lyap = [];
% total number of iterations
N = length(t);
% store convergence information
lyap_conv = zeros(3,N);

for k = 1:N
x0 = x(k,:);  % new x0 after each step
Mn = (eye(3)+ jac_lorenz(x0,par)*dt) * M;

[Q,R] = qr(Mn);
lyap = [lyap log(abs(diag(R)))];

if k == 1
lyap_conv(:,k) = lyap(:,k);
else
lyap_conv(:,k) = 1/(k*dt)*sum(lyap')';
end

M = Q;
end

L = 1/(N*dt)*sum(lyap',1)

%%% convergence plot
figure(2)
plot(t,lyap_conv(1,:),'r-','LineWidth',2)
hold on
plot(t,lyap_conv(2,:),'k-','LineWidth',2)
hold on
plot(t,lyap_conv(3,:),'b-','LineWidth',2)
xlabel('time (s)')
ylabel('Lyapunov exponents')
set(gca,'FontSize',15)

% Kaplan-Yorke dimension
Dky = 2+lyap(1)/abs(lyap(3))

function xdot = lorenz(t,x,par)
xdot = [par(1)*(x(2)-x(1));
x(1)*(par(2)-x(3))-x(2);
x(1)*x(2)-par(3)*x(3)];

function J = jac_lorenz(x0,par)
J = [-par(1) par(1) 0;
par(2)-x0(3) -1 -x0(1);
x0(2) x0(1) -par(3)];



When I change the final time in tspan, say from 50 to 150, I get the following results:

for 50:

lyap =

0.7764    0.0103  -14.4829

Dky =

2.0536


for 100:

lyap =

0.8613    0.0294  -14.5915

Dky =

2.0590


for 150:

lyap =

0.8830    0.0379  -14.6224

Dky =

2.0604


These results are similar, but I do not think each individual Lyapunov exponent is converging to a single value. In particular, the 2nd Lyapunov exponent, which I believe should be 0, is actually increasing as the time becomes longer... So if this was a completely new system, how do I know this second exponent is actually zero, rather than the system being hyper-chaotic (having more than 1 Lyapunov exponents > 0)?

I am therefore curious to learn if I were to report the 3 Lyapunov exponents for this particular system, i.e. equations, parameters (sigma, r, b, for the Lorenz system), initial conditions FIXED, what values should I report?

Say if I want to do a parameter space study. For example, I want to learn how the system changes (via Lyapunov exponents) as I vary only sigma. Now I need to run a large number of simulations (one for each sigma value I want to test), and therefore I would like to run each of these simulations with the shortest time, biggest time step size, and an initial condition that is as close to the attractor as possible. Basically, how can I find the convergence of Lyapunov exponents efficiently?

Lastly, I have not tried this myself, but I have been thinking about the cases where Lyapunov exponents are not defined (e.g. for a fixed point or a limit cycle). For example, for a 3D system with a stable/unstable spiral fixed point, I assume the 3 Lyapunov exponents will keep decreasing/increasing as the time becomes longer and longer, i.e. similar to what I showed above with the Lorenz attractor, so how can I tell here that these Lyapunov exponents are actually NOT converging, and therefore are undefined. Again, I would like to know this with the least amount of computing (by decreasing final simulation time, etc.) without compromising reliability.

As an update, I made the corresponding convergence plots for the 3 different final times:   These plots do show convergence! Hence I think I can apply some standard convergence test on these and report the values as soon as the convergence test is satisfied. (PS. In a previous version, the plots were not converging. This was because at each step the solution should be weighted by its step number, not the total step number. I have corrected this.)

The object to analyze is the step Jacobian $$I+J_k\,dt.$$ $$J_k$$ should correctly be $$\nabla_x\Phi_F(t_k,x_k,dt)$$, the Jacobian of the step update $$x_k+\Phi_F(t_k,x_k,dt)\,dt$$, but in a pinch can be approximated with the Jacobian $$\nabla_x F(t_k,x_k)$$ of the ODE system function.

As $$J_{k+1}-J_k=O(dt)$$ if the ODE is smooth enough, the step Jacobian is slowly changing, with difference $$O(dt^2)$$ from one step to the next. Over $$O(1/dt)$$ steps it could be considered constant. Then the operations amount to the QR algorithm $$A_k=M_k^T(I+J_k\,dt)M_k\\ Q_k·R_k=A_k\\ M_{k+1}=M_kQ_k$$ so that, approximately $$A_{k+1}=Q_k^TM_k^T(I+J_{k+1}\,dt)M_kQ_k\approx Q_k^TA_kQ_k=R_kQ_k$$ The structure of the matrices $$A_k$$ converge to an upper triangular form with a block diagonal. Diagonal blocks result from complex eigenvalues and eigenvalues with multiplicity. Inside non-trivial blocks there will be no convergence. The diagonal of $$R_k$$ will, hopefully and on average, correspond to the absolute values of the eigenvalues.

With all the handwaving involved, it should become clear that an error of $$O(dt)$$ from some ideal result is to be expected.

Also, to get sensible results the initial matrix $$M_1$$ should be the result of several steps of the QR algorithm for the matrix $$I+J_1\,dt$$. Otherwise the initial convergence phase can greatly delay the convergence of the average.

The initial point for the ODE solution should be close to the attractor, which the point $$[1,1,1]$$ is not. This also influences the result, initializing the averaging process with the result of a different dynamic. Compute a throw-away intro segment to get close to the attractor, start close to one of the saddle points or at the midpoint between them.

• I am still wondering: 1. what values should I report based on my computation? 2._ "Inside non-trivial blocks there will be no convergence."_ I am not really understanding this, but I think you are talking about how to find non-convergence here, so if you could elaborate more on this, it will be very much appreciated! 3.Regarding your last point, I want to confirm this only increases the speed of convergence, and doesn't influence much the final values right (assuming large enough final T)? For a new system (high-dimensional), is there a general way to find the rough location of an attractor? Jun 20 at 0:18
• 1.) It seems that only one decimal after the dot is reliable. With the proposed corrections, this could be increased to 3 or 4 decimal places. Note that visual convergence usually amount to about 1% error. 2.) The diagonal blocks correspond to invariant subspaces. There is no unique basis, and the QR process will change the basis vectors from step to step. 3.) The last points serve to cut away the initial transients. The exponents do not change. The diagonal entries of R will oscillate along the path, this evens out in the average, but that is a slow convergence. Jun 20 at 3:47
• For the location of an attractor, you could check if the evolution of the solution becomes finally constrained to a small volumen of the phase space. Sometimes a larger such volume can be derived from the equations. I'm afraid this becomes very empirical and heuristical. Jun 20 at 3:54
• Thanks for the help! I am still unclear about these: 4. it sounds like the uncertainties come from linearization of the fiducial line (i.e. the use of the Jacobian of the system; correct me if I misunderstood you), but not Euler method? If I only increase the order for the integrator and leave the Jacobian alone, will I improve my accuracy? More so, I am not seeing in your proposed method that you are including higher order terms for the linearization, so I am afraid I am still missing something. Jun 20 at 5:07
• 5. regarding non-convergence. If I have an n-dimensional system that has a 1D curve as an attractor (a limit cycle), I should have Lyapunov spectrum of one 0 and (n-1) number of undefined/non-converging values. When can I declare that these are indeed non-converging (and the attractor is a 1D limit cycle), while ruling out the possibility that I just didn't integrate long enough... sorry, I am not seeing a practical way of doing this based on your discussion of the structure of the R matrix. Jun 20 at 5:13