fig4_1.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

b = 100;                 
h = 0.001;            % 积分时间步长

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

rmaxmin=0;  KK1=0; KK2=0;

%brc=zeros(length(200:850),length(1:0.01:1.6));
sigma_ar=200:1:850;
pp_ar=1:0.001:1.6;
kk1=0;
for cc=1:length(sigma_ar);
    sigma=sigma_ar(cc);
    rc=sigma*(sigma+b+3)/(sigma-b-1);
    kk1=kk1+1;
    kk2=0;
    for pp=1:length(pp_ar)
        kk2=kk2+1;
        y = [265,-384,590]./1.65.*pp_ar(pp);     
        r=rc*pp_ar(pp);
        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))));
    %    brc(kk1,kk2) = max(rr(1:floor(k2/2)-1));
%         if brc(kk1,kk2)< rmaxmin
        % plot(r);
%             rmaxmin=brc(kk1,kk2);
%             KK1=kk1;
%             KK2=KK2;
%         end
     Xp=max(X);
     Xm=sqrt(sum((X.^2))/length(X));
     pv=Xp/Xm;
     brc(kk2,kk1) = pv;
    end
    sigma
end
%plot(pp_ar,sigma_ar,brc);
pcolor(sigma_ar,pp_ar,brc)
colorbar;
colormap gray
colormap hot
%colormap(hsv(128))

nn=brc<1.9;
nn=sum(sum(nn))

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;