% 計算混沌時間序列 Lyapunov 指數 - Henon 序列% 計算 Henon 映射 x(n+1) = 1 - a * x(n)^2 + y(n); y(n+1) = b * x(n)clcclear allclose alldisp('********** Henon ***********')%-----------------------------------------------------------------% 產生 Henon 混沌時間序列a = 1.4;b = 0.3;x0 = 0;y0 = 0;k1 = 20000;                % 前面迭代次數k2 = 5000;                % 后面迭代次數X = zeros(k1+k2,1);for i = 1:k1+k2    x1 = 1 - a * x0^2 + y0;    y1 = b * x0;    x0 = x1;    y0 = y1;    X(i) = x1;endX = X(k1+1:end);%-----------------------------------------------------------------% Henon 入口參數fs = 1;                 % 采樣頻率t = 1;                  % 時延d = 2;                  % 嵌入維tmax = 30;              % 最大離散步進時間p = 100;                % 序列平均周期%-----------------------------------------------------------------% 調用計算函數Y = Lyapunov_rosenstein_2(X,fs,t,d,tmax,p);I = [0:length(Y)-1]';linear = [2:9]';  % 線性區域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)