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

Question

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.)

Answer 1

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.