Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Larges Lyapunov Exponent from time series
function D = dist(X)
s = length(X(1,:));
n = length(X(:,1));
D = zeros(s,s);
for i=1:s
for j=1:i
if i == j
D(i,j) = 0;
else
D(i,j) = sqrt(sum(power(X(:,i)-X(:,j),2)));
endif
endfor
endfor
D = D + D' + eye(s,s)*1e6;
endfunction
%
% findLyapunov()
% funkcja znajduje najwiekszy wykladnik Lapunowa dla trajektorii zanurzonej w 3 wymiarach
% x, y, z - wspolrzedne punktow odleglych o jednakowy odstep czasu
% autor: Mateusz Galazyn
% oparte na algorytmie z:
% WOLF, SWIFT, SWINNEY, AND VASTANO "Determining lyapunov exponents from a time series."
% Physica D: Nonlinear Phenomena 16, 3 (July 1985), 285–317.
%
function L = findLyapunov(x,y,z)
s = length(x);
D = dist([x' ; y' ; z' ]);
L = 0;
counter = 0; % licznik punktow dla ktorych zostala obliczona ekspansja
for i=1:floor(s-1)
[v, iv] = min(D(i,:)); % najblizszy sasiad
if iv < s
distancePrime = dist( [x(i+1) x(iv+1) ; y(i+1) y(iv+1) ; z(i+1) z(iv+1) ] )(1,2); % L prime
distance = dist( [x(i) x(iv) ; y(i) y(iv) ; z(i) z(iv) ] )(1,2); % L
if distance > 0 && distancePrime > 0
++counter;
L += log2( distancePrime / distance );
endif
endif
endfor
if counter == 0
L = 0; % jezeli orbita okresowa, zaden punkt nie zostal uzyty do obliczen
else
L /= counter;
endif
endfunction
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment