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