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:
lyap = 0.7764 0.0103 -14.4829 Dky = 2.0536
lyap = 0.8613 0.0294 -14.5915 Dky = 2.0590
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.
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.)