绘制庞加莱截面图.txt

绘制庞加莱截面图,在相空间中适当（要有利于观察系统的运动特征和变化

Poincare截面 
　　在相空间中适当（要有利于观察系统的运动特征和变化，如截面不能与轨线相切，更不能包含轨线）选取一截面，在此截面上某一对共轭变量如x1和x.1取固定值，称此截面为Poincare截面，相空间的连续轨迹与Poincare截面的交点成为截点。通过观察Poincare截面上截点的情况可以判断是否发生混沌：当Poincare截面上有且只有一个不动点或少数离散点时，运动是周期的；当Poincare截面上是一封闭曲线时，运动是准周期的；当Poincare截面上是一些成片的具有分形结构的密集点时，运动便是混沌。 

matlab计算程序如下: 

文件一，其文件名为Poincare.m 

function dx=Poincare(t,x); 
% 单摆方程[不显含时间t的自治系统] 
% 方程如下: 
% dθ/dt=ω, 
% dω/dt=-2*β*[dθ/dt]-ω^2*sin(θ)+F*cos(vt) 
% dψ/dt=v 
betaa=0.25; 
F=1.093; 
v=2/3; 
P2=-2*betaa*x(2)-x(2).^2.*sin(x(1))+F*cos(v*t); 
dx=[x(2);P2;v]; 

文件二，其文件名为Poincare_section.m 

% 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); 

%%% 另外用下面一个文件也可以实现的
% Poincare_section[绘制庞加莱截面图] 
betaa=0.25; 
F=1.093; 
v=2/3; 
Poin=inline(['[x(2);',... 
       '-2*betaa*x(2)-x(2).^2.*sin(x(1))+F*cos(v*t);',... 
       'v]'],... 
       't','x','flag','betaa','F','v'); 
% Poincare_section[绘制庞加莱截面图] 
[t,x]=ode45(Poin,[0,2800],[0,1.5,0],[],betaa,F,v); 
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);