[转帖]绘制庞加莱截面图
Poincare截面在相空间中适当（要有利于观察系统的运动特征和变化，如截面不能与轨线相切，更不能包含轨线）选取一截面，在此截面上某一对共轭变量如x1和x.1取固定值，称此截面为Poincare截面，相空间的连续轨迹与Poincare截面的交点成为截点。通过观察Poincare截面上截点的情况可以判断是否发生混沌：当Poincare截面上有且只有一个不动点或少数离散点时，运动是周期的；当Poincare截面上是一封闭曲线时，运动是准周期的；当Poincare截面上是一些成片的具有分形结构的密集点时，运动便是混沌。
matlab计算程序如下:
文件一，其文件名为Poincare.m
function dx=Poincare(t,x);
% 单摆方程[不显含时间t的自治系统]
% 方程如下:
% dθ/dt=ω,
% dω/dt=-2*β*-ω^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=;
文件二，其文件名为Poincare_section.m
% Poincare_section[绘制庞加莱截面图]
=ode45(@Poincare,,);
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;
=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(,,1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit(,,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);
[ 本帖最后由 suffer 于 2006-10-9 19:29 编辑 ] %%% 另外用下面一个文件也可以实现的
% 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[绘制庞加莱截面图]
=ode45(Poin,,,[],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;
=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(,,1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit(,,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);
[ 本帖最后由 suffer 于 2006-10-9 19:29 编辑 ] 关键是实例麻烦在象我这种的人看不懂啊：（
呵呵，但还是顶，最好有点解释：） 我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了 楼主，我刚接触混沌，请问这个程序能否利用原始的时间序列做出ＰＯＩＮＣＡＲＥ散点图啊？？
另外请教一下，ＭＡＴＬＡＢ中有非线性分析库吗？？
谢谢
好
我感觉这个版块需要的高手要很多最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了 同意7楼的，最好来点注释。
还有闪频法画的相图跟ＰＯＩＮＣＡＲＥ散点图有什么区别？
这个好像有人问了，但是没有看到答案？
谢谢 非常感谢楼主的程序，但还是有些看不懂。怎么把实验数据导入，得到图形，望多多指教，谢谢了！ 这个才是正确的！
那个萝卜网站上的不正确阿！不要上当！ 版主你好！谢谢您辛苦了！我目前在解转子碰摩的微分方程，也做它的poincare映射图，程序还有看不懂的地方，希望能和你一起学习！ 你好！我下载下来了，可以运行，不过我觉得还不是太混沌！不知道你用的参数是进入混沌的参数么，还有
for k=1:round(max(x(:,3))/2/pi);
d=x(:,3)-(k-1)*2*pi-phi0;
=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(,,1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit(,,1);
X2(k)=polyval(Q,(k-1)*2*pi-phi0);
end
end
这段我不明白！真心希望请教你！ 不知道你用的是否是混沌参数么！ 我运行了不过不明白这段是什么意思，想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了 我运行了不过不明白这段是什么意思，
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;
=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(,,1);
X1(k)=polyval(Q,(k-1)*2*pi-phi0);
Q=polyfit(,,1);
X2(k)=polyval(Q,(k-1)*2*pi-phi0);
end
想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了 shuijing :那一段程序是在求曲线和phi=2*pi/3这个截面的交点。
页:
[1]
2