% 計算混沌時間序列 Lyapunov 指數 - Lorenz 吸引子%clcclearclose alldisp('********** Lorenz ***********')%-----------------------------------------------------------------% 產生 Lorenz 混沌時間序列% dx/dt = sigma*(y-x)% dy/dt = r*x - y - x*z% dz/dt = -b*z + x*ysigma = 16;             % Lorenz 方程參數 a = 16 | 10b = 4;                  %                 b = 4 | 8/3r = 45.92;              %                 c = 45.92 | 28           y = [-1,0,1];           % 起始點 (1 x 3 的行向量)h = 0.01;               % 積分時間步長k1 = 30000;             % 前面的迭代點數k2 = 5000;              % 后面的迭代點數Z = LorenzData(y,h,k1+k2,sigma,r,b);X = Z(k1+1:end,1);      % 時間序列%-----------------------------------------------------------------% Lorenz 入口參數fs = 1/h;               % 采樣頻率t = 10;                 % 時延d = 3;                  % 嵌入維tmax = 300;             % 最大離散步進時間p = 50;                 % 序列平均周期%-----------------------------------------------------------------% 調用計算函數Y = Lyapunov_rosenstein_2(X,fs,t,d,tmax,p);I = [0:length(Y)-1]';linear = [60:180]';  % 線性區域F1 = polyfit(I(linear),Y(linear),1);Y1 = polyval(F1,I(linear),1);%-----------------------------------------------------------------% 結果顯示figuresubplot(211); plot(I,Y,'k-'); grid; xlabel('i'); ylabel('y(i)'); hold on;plot(I(linear),Y1,'r-'); hold off;subplot(212); plot(I(1:end-1),diff(Y),'k-'); grid; xlabel('n'); ylabel('slope');Lyapunov1_2 = F1(1)Lyapunov1_e = F1(1)*log(2)