最新消息: USBMI致力于为网友们分享Windows、安卓、IOS等主流手机系统相关的资讯以及评测、同时提供相关教程、应用、软件下载等服务。

如何绘制庞加莱映射

IT圈 admin 28浏览 0评论

2024年4月26日发(作者:怀晖)

绘制庞加莱截面图

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)2;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('dtheta/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('dtheta/dt','fontsize',14);

*******************************************************************************

2024年4月26日发(作者:怀晖)

绘制庞加莱截面图

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)2;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('dtheta/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('dtheta/dt','fontsize',14);

*******************************************************************************

发布评论

评论列表 (0)

  1. 暂无评论