main_lorenz_2.m

用matlab编程产生非线性动力系统的时间序列

clc
clear all
close all

%--------------------------------------------------------------------------
% 方程表达式
% dx/dt = sigma*(y-x)
% dy/dt = r*x - y - x*z
% dz/dt = -b*z + x*y

Z=1;

sigma = 441 *Z;             % Lorenz 方程参数 a = 16 | 10
b = 165 *Z;                  %                 b = 4 | 8/3
rc=sigma*(sigma+b+3)/(sigma-b-1);
r = rc*1.1;              %                 c = 45.92 | 28    
r=982 *Z;

%y = [265,-384,590];        % 起始点 (1 x 3 的行向量)

y = [265,-384-50+41+1,1966-50+75+1];        % 起始点 (1 x 3 的行向量)

h = 0.001;            % 积分时间步长

k1 = 10000;           % 前面的迭代点数
k2 = 1000;           % 后面的迭代点数

rmaxmin=-10;
kkkk=0;

for ll=1:100;
   % y(2)=y(2)+1;
for kk=1:100
   %  y(3)=y(3)+1;

    data = LorenzData(y,h,k1+k2,sigma,r,b);
    data = data(k1+1:end,:);

    %--------------------------------------------------------------------------

    X = data(:,1);
    Y = data(:,2);
    Z = data(:,3);

    figure
    plot3(Z,Y,X);
    xlabel('Z');ylabel('Y');zlabel('X');
    title('Lorenz attractor');
%     
%     figure(2)
%     xx=abs(fft(X(1:1:end)));
%     plot(20*log10(xx/max(xx)));

    rr=xcorr(X,X);
    rr=20*(log10(abs(rr)/max(abs(rr))));
    % plot(r);
    rmm(ll,kk) = max(rr(1:floor(k2/2)-1));
    if rmm(ll,kk)<rmaxmin
        rmaxmin=rmm(ll,kk);
        kkkk=kk;
        jjjj=ll;
    end
end
ll
end
plot(rmm);
figure
pcolor(rmm)
colorbar;
colormap bone
colormap hot
colormap(hsv(256))

load x
figure(3)
xx=abs(fft(x1(1:1:end)));
plot(20*log10(xx/max(xx)),'r--');

hold on
xx=abs(fft(x2(1:1:end)));
plot(20*log10(xx/max(xx)),'r--');

xx=abs(fft(x3(1:1:end)));
plot(20*log10(xx/max(xx)),'k-');

xx;