[转帖]绘制庞加莱截面图

suffer

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

matlab计算程序如下:

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

1. function dx=Poincare(t,x);
   
2. % 单摆方程[不显含时间t的自治系统]
   
3. % 方程如下:
   
4. % dθ/dt=ω,
   
5. % dω/dt=-2*β*[dθ/dt]-ω^2*sin(θ)+F*cos(vt)
   
6. % dψ/dt=v
   
7. betaa=0.25;
   
8. F=1.093;
   
9. v=2/3;
   
10. P2=-2*betaa*x(2)-x(2).^2.*sin(x(1))+F*cos(v*t);
    
11. dx=[x(2);P2;v];
    
12. 
    
13. 文件二，其文件名为Poincare_section.m
    
14. 
    
15. % Poincare_section[绘制庞加莱截面图]
    
16. [t,x]=ode45(@Poincare,[0,2800],[0,1.5,0]);
    
17. x(:,2)=mod(x(:,2),2*pi)-pi;
    
18. 
    
19. phi0=pi*2/3; % 选择phi=2*pi/3这个截面
    
20. for k=1:round(max(x(:,3))/2/pi);
    
21. d=x(:,3)-(k-1)*2*pi-phi0;
    
22. [P,K]=sort(abs(d));
    
23. x1l=x(K(1),1);
    
24. x1r=x(K(2),1);
    
25. x2l=x(K(1),2);
    
26. x2r=x(K(2),2);
    
27. x3l=x(K(1),3);
    
28. x3r=x(K(2),3);
    
29. if abs(P(1))+abs(P(2))<3e-16;
    
30. X1(k)=x1l;
    
31. X2(k)=x2l;
    
32. else
    
33. Q=polyfit([x3l,x3r],[x1l,x1r],1);
    
34. X1(k)=polyval(Q,(k-1)*2*pi-phi0);
    
35. Q=polyfit([x3l,x3r],[x2l,x2r],1);
    
36. X2(k)=polyval(Q,(k-1)*2*pi-phi0);
    
37. end
    
38. end
    
39. plot(X1,X2,'.');
    
40. xlabel('\theta','fontsize',14);
    
41. ylabel('d\theta/dt','fontsize',14);

复制代码

[ 本帖最后由 suffer 于 2006-10-9 19:29 编辑 ]

suffer

%%% 另外用下面一个文件也可以实现的

1. % Poincare_section[绘制庞加莱截面图]
   
2. betaa=0.25;
   
3. F=1.093;
   
4. v=2/3;
   
5. Poin=inline(['[x(2);',...
   
6. '-2*betaa*x(2)-x(2).^2.*sin(x(1))+F*cos(v*t);',...
   
7. 'v]'],...
   
8. 't','x','flag','betaa','F','v');
   
9. % Poincare_section[绘制庞加莱截面图]
   
10. [t,x]=ode45(Poin,[0,2800],[0,1.5,0],[],betaa,F,v);
    
11. x(:,2)=mod(x(:,2),2*pi)-pi;
    
12. 
    
13. phi0=pi*2/3; % 选择phi=2*pi/3这个截面
    
14. for k=1:round(max(x(:,3))/2/pi);
    
15. d=x(:,3)-(k-1)*2*pi-phi0;
    
16. [P,K]=sort(abs(d));
    
17. x1l=x(K(1),1);
    
18. x1r=x(K(2),1);
    
19. x2l=x(K(1),2);
    
20. x2r=x(K(2),2);
    
21. x3l=x(K(1),3);
    
22. x3r=x(K(2),3);
    
23. if abs(P(1))+abs(P(2))<3e-16;
    
24. X1(k)=x1l;
    
25. X2(k)=x2l;
    
26. else
    
27. Q=polyfit([x3l,x3r],[x1l,x1r],1);
    
28. X1(k)=polyval(Q,(k-1)*2*pi-phi0);
    
29. Q=polyfit([x3l,x3r],[x2l,x2r],1);
    
30. X2(k)=polyval(Q,(k-1)*2*pi-phi0);
    
31. end
    
32. end
    
33. plot(X1,X2,'.');
    
34. xlabel('\theta','fontsize',14);
    
35. ylabel('d\theta/dt','fontsize',14);

复制代码

[ 本帖最后由 suffer 于 2006-10-9 19:29 编辑 ]

cdwxg

关键是实例麻烦在象我这种的人看不懂啊：（
呵呵，但还是顶，最好有点解释：）

cdwxg

我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了

juan8909

楼主，我刚接触混沌，请问这个程序能否利用原始的时间序列做出ＰＯＩＮＣＡＲＥ散点图啊？？
另外请教一下，ＭＡＴＬＡＢ中有非线性分析库吗？？
谢谢

coldspring

我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了

无水1324

同意7楼的，最好来点注释。
还有闪频法画的相图跟ＰＯＩＮＣＡＲＥ散点图有什么区别？
这个好像有人问了，但是没有看到答案？
谢谢

yfh2006

非常感谢楼主的程序，但还是有些看不懂。怎么把实验数据导入，得到图形，望多多指教，谢谢了！

ThomasLee

这个才是正确的！
那个萝卜网站上的不正确阿！不要上当！

shuijing

版主你好！谢谢您辛苦了！我目前在解转子碰摩的微分方程，也做它的poincare映射图，程序还有看不懂的地方，希望能和你一起学习！

shuijing

你好！我下载下来了，可以运行，不过我觉得还不是太混沌！不知道你用的参数是进入混沌的参数么，还有
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
这段我不明白！真心希望请教你！

shuijing

不知道你用的是否是混沌参数么！

shuijing

我运行了不过不明白这段是什么意思，想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了

shuijing

我运行了不过不明白这段是什么意思，
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

想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了

toes

shuijing :那一段程序是在求曲线和phi=2*pi/3这个截面的交点。