LyapSpec.zip

Calculates full spectrum of Lyapunov exponents or k first exponents

This method was presented on 11th Workshop on Optimal Control, Dynamic Games and
Nonlinear Dynamics in Amsterdam, 2010 by Anton O. Беляков
function [LE,trJ,x] = LyapSpec(fun,T,x0,r,k)
% Calculates Lyapunov exponents (LE) with Gram-Schmidt ortonormalization at each step of second order solver comprising Heun's method for linearized system and Leapfrog (Verlet) method for state equation.
% fun is the function of a dynamical system (n) and its Jacobian matrix (n x n)
% T is the time interval: T(1) is the start time, T(2) is the end time,
% x0 is the vector of the initial state variables (n),
% r is the sampling rate
% k is the number of LE to be calculated, (1)
% Returns:
% LE is the vector of Lyapunov exponents (k),
% trJ is sum of all LE (1),
% x is the final state variables (n).

%Example of function fun:

function [dz,dzz] = Roessler(t,z)
% Rossler system
global a b c d
% function of the right-hand side
dz = [-(z(2) + z(3))
z(1) + a * z(2) + z(4)
b + z(1) * z(3)
c * z(4) - d * z(3)];
% Jacobian of the right-hand side function
dzz = [0, -1, -1, 0
1, a, 0, 1
z(3), 0, z(1), 0
0, 0, -d, c];
% t is not used and introduced for compatibility with other solvers
end

% main file to calculate Lyapunov spectrum of Rössler system with hyperchaos for different sampling rates
clc;
clear;
global a b c d
a = 0.25;
b = 3.0;
c = 0.05;
d = 0.5;

r = 100; % initial sampling rate
T = 100000; % timespan of calculation
x0 = [-20;0;0;15]; % initial conditions
% tic
% [LE,trJ,x0] = LyapSpec(@Roessler,[0 5000],x0,r,0); % approach the attractor
% toc
fprintf('T = %16.0f\n',T);
NumLE = length(x0); % number of Lyapunov Exponents
fprintf('dt');
for k = 1:NumLE, fprintf('\t LE_%d',k); end
fprintf('\t traceJ\n');

tic
for k = 1:4
[LE,trJ] = LyapSpec(@Roessler,[0 T],x0,r,NumLE);
fprintf('%16.8f',1/r);
fprintf('\t %16.8f',[LE', trJ]);
fprintf('\n');
r = r*2; % doubling the sampling rate
end
toc

Result:

T = 100000
dt LE_1 LE_2 LE_3 LE_4 traceJ
0.01000000 0.11144250 0.02013357 -0.00010669 -23.29383296 -24.78187348
0.00500000 0.11141259 0.02143191 -0.00006508 -24.67123402 -24.86829197
...

Notice, that the third Lyapunov exponent (LE_3) theoretically should be zero. Thus, its deviation from 0 reflects precision of the method in hand.