Henon系统 （2）

Henon.m

function newx=Henon(x); % Henon方程[差分方程] % newx=Henon([x;y;p;q]) % % 方程如下: % x(k+1)=-p*x(k)*x(k)+y(k)+1 % y(k+1)=q*x(k) % % Example(函数图象): % % X=[1;0;1.4;.3];Y=[]; % for i=1:1000 % X=feval(@Henon,X); % Y(i,:)=X(1:2,1); % end % plot(Y(:,1),Y(:,2),'.'); % title('Henon映射图'),xlabel('x'),ylabel('y') % % Example(分岔图): % % Z=[]; % for p=linspace(0,1.4,300); % x=[1;0;p;.3]; % for k=1:300; % x=Henon(x); % if k>60 % Z=[Z,p+x(1)*i]; % end % end % end % plot(Z,'.','markersize',1) % title('Henon映射分岔图'),xlabel('p'),ylabel('x') % % Example(最大Lyapunov指数谱图): % % d0=1e-8; % Z=[]; % for p=linspace(0,1.4,300) % le=0; % lsum=0; % x=[1;0;p;.3]; % x1=[1;d0;p;.3]; % for k=1:800 % x=Henon(x); % x1=Henon(x1); % d1=sqrt((x(1)-x1(1))^2+(x(2)-x1(2))^2); % x1=x+(d0/d1)*(x1-x); % if k>100 % lsum=lsum+log(d1/d0); % end % end % le=lsum/(k-100); % Z=[Z,p+le*i]; % end % plot(Z,'-') % title('Henon最大Lyapunov指数谱图'),xlabel('p'),ylabel('lyapunov') % grid on % % Author：yujunjie % Author's email: ustb03-07@yahoo.com.cn newx(1,1)=-x(3)*x(1)*x(1)+x(2)+1; newx(2,1)=x(4)*x(1); newx(3,1)=x(3); newx(4,1)=x(4);