lyapunov.m

非线性动力学中lyapunov、duffing等方程的求解

function lyapunov(filename)

% load in data file
data = load(filename);
% calculate number of data points
N = length(data);
% create another array for return map
nextdata = zeros(N,1);
for (i=1:N-1)
  nextdata(i)=data(i+1);
end
nextdata(N)=nextdata(N-1);
% now plot it (except last point)
plot(data(1:N-1),nextdata(1:N-1),'.');
% give it a title
title('return map');
display('press any key to continue');
% wait for user input ...
pause;

% calculate the fast fourier transform
fftcoeff = fft(data,N);
t = 2:N;
% plot the absolute value fourier coefficients
plot(t,abs(fftcoeff(2:N)));
% give it a title
title('Fourier Coefficients');
display('press any key to continue');
% wait for user input
pause;

N2 = floor(N/2);
N4 = floor(N/4);
% find mid point of orbit sequence
k=N2;
% create space for exponents
exponent = zeros(N4,1);
% look at 1/4 of the points
for (j=1:N4)
    % set distance initially
    d = abs(data(k+1)-data(k));
    index = k+1;
    for (i=2:N-1)
        % see if there is a closer point
        if (i ~= k) && (abs(data(i)-data(k)))<d
            % if so, store index and distance
            d = abs(data(i)-data(k));
            index = i;
        end
    end
    % write log of quotient as difference of logs to get better accuracy!
    if (data(k) ~= data(index)) && (data(k+1) ~= data(index+1))
        exponent(j) = log( abs(data(k+1)-data(index+1)))-log(abs(data(k)-data(index)));
    end
    % repeat with the next point
    k = k+1;
end

% now plot the lyapunov exponents
t = 1:N4;
lyapunov = exponent(1:N4);
exp_avg = 0.0;
% find the average value for lyapunov exponent
for (i=1:N4)
    exp_avg = exp_avg + exponent(i);
end
% plot the exponents, the average, and the baseline
exp_avg = exp_avg/N4;
plot(t,lyapunov,t,0,t,exp_avg);
display 'average value for lyapunov exponent is:'
exp_avg