poincare_section.m

来自「分形计算的几个程序原码」· M 代码 · 共 27 行

% Poincare_section[绘制庞加莱截面图] 
[t,x]=ode45(@Poincare,[0,2800],[0,1.5,0]); 
x(:,2)=mod(x(:,2),2*pi)-pi; 

phi0=pi*2/3; % 选择phi=2*pi/3这个截面 
for k=1:round(max(x(:,3))/2/pi); 
   d=x(:,3)-(k-1)*2*pi-phi0; 
   [P,K]=sort(abs(d)); 
   x1l=x(K(1),1); 
   x1r=x(K(2),1); 
   x2l=x(K(1),2); 
   x2r=x(K(2),2); 
   x3l=x(K(1),3); 
   x3r=x(K(2),3); 
   if abs(P(1))+abs(P(2))<3e-16; 
       X1(k)=x1l; 
       X2(k)=x2l; 
   else 
       Q=polyfit([x3l,x3r],[x1l,x1r],1); 
       X1(k)=polyval(Q,(k-1)*2*pi-phi0); 
       Q=polyfit([x3l,x3r],[x2l,x2r],1); 
       X2(k)=polyval(Q,(k-1)*2*pi-phi0); 
   end 
end 
plot(X1,X2,'.'); 
xlabel('\theta','fontsize',14); 
ylabel('d\theta/dt','fontsize',14);