查看: 13299|回复: 23

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

[复制链接]
发表于 2006-4-7 11:27 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
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. 文件二,其文件名为Poincare_section.m

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

  16. phi0=pi*2/3; % 选择phi=2*pi/3这个截面
  17. for k=1:round(max(x(:,3))/2/pi);
  18. d=x(:,3)-(k-1)*2*pi-phi0;
  19. [P,K]=sort(abs(d));
  20. x1l=x(K(1),1);
  21. x1r=x(K(2),1);
  22. x2l=x(K(1),2);
  23. x2r=x(K(2),2);
  24. x3l=x(K(1),3);
  25. x3r=x(K(2),3);
  26. if abs(P(1))+abs(P(2))<3e-16;
  27. X1(k)=x1l;
  28. X2(k)=x2l;
  29. else
  30. Q=polyfit([x3l,x3r],[x1l,x1r],1);
  31. X1(k)=polyval(Q,(k-1)*2*pi-phi0);
  32. Q=polyfit([x3l,x3r],[x2l,x2r],1);
  33. X2(k)=polyval(Q,(k-1)*2*pi-phi0);
  34. end
  35. end
  36. plot(X1,X2,'.');
  37. xlabel('\theta','fontsize',14);
  38. ylabel('d\theta/dt','fontsize',14);
复制代码

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

评分

1

查看全部评分

回复
分享到:

使用道具 举报

 楼主| 发表于 2006-4-7 11:28 | 显示全部楼层
%%% 另外用下面一个文件也可以实现的
  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. phi0=pi*2/3; % 选择phi=2*pi/3这个截面
  13. for k=1:round(max(x(:,3))/2/pi);
  14. d=x(:,3)-(k-1)*2*pi-phi0;
  15. [P,K]=sort(abs(d));
  16. x1l=x(K(1),1);
  17. x1r=x(K(2),1);
  18. x2l=x(K(1),2);
  19. x2r=x(K(2),2);
  20. x3l=x(K(1),3);
  21. x3r=x(K(2),3);
  22. if abs(P(1))+abs(P(2))<3e-16;
  23. X1(k)=x1l;
  24. X2(k)=x2l;
  25. else
  26. Q=polyfit([x3l,x3r],[x1l,x1r],1);
  27. X1(k)=polyval(Q,(k-1)*2*pi-phi0);
  28. Q=polyfit([x3l,x3r],[x2l,x2r],1);
  29. X2(k)=polyval(Q,(k-1)*2*pi-phi0);
  30. end
  31. end
  32. plot(X1,X2,'.');
  33. xlabel('\theta','fontsize',14);
  34. ylabel('d\theta/dt','fontsize',14);
复制代码


[ 本帖最后由 suffer 于 2006-10-9 19:29 编辑 ]
发表于 2006-4-8 14:27 | 显示全部楼层
关键是实例麻烦在象我这种的人看不懂啊:(
呵呵,但还是顶,最好有点解释:)
发表于 2006-4-9 00:03 | 显示全部楼层
我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开?或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释,对我们水平差的太痛苦了
发表于 2006-4-9 14:09 | 显示全部楼层
楼主,我刚接触混沌,请问这个程序能否利用原始的时间序列做出POINCARE散点图啊??
另外请教一下,MATLAB中有非线性分析库吗??
谢谢
发表于 2006-4-10 10:14 | 显示全部楼层

我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开?或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释,对我们水平差的太痛苦了
发表于 2006-4-16 16:58 | 显示全部楼层
同意7楼的,最好来点注释。
还有闪频法画的相图跟POINCARE散点图有什么区别?
这个好像有人问了,但是没有看到答案?
谢谢
发表于 2006-4-26 16:52 | 显示全部楼层
非常感谢楼主的程序,但还是有些看不懂。怎么把实验数据导入,得到图形,望多多指教,谢谢了!
发表于 2006-5-6 21:16 | 显示全部楼层
这个才是正确的!
那个萝卜网站上的不正确阿!不要上当!
发表于 2006-5-12 08:33 | 显示全部楼层
版主你好!谢谢您辛苦了!我目前在解转子碰摩的微分方程,也做它的poincare映射图,程序还有看不懂的地方,希望能和你一起学习!
发表于 2006-5-15 16:44 | 显示全部楼层
你好!我下载下来了,可以运行,不过我觉得还不是太混沌!不知道你用的参数是进入混沌的参数么,还有
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
这段我不明白!真心希望请教你!
发表于 2006-5-15 16:48 | 显示全部楼层
不知道你用的是否是混沌参数么!
发表于 2006-5-17 08:48 | 显示全部楼层
我运行了不过不明白这段是什么意思,想请教一下,还有我觉得你选的参数是混沌的参数,可是出来的图并不是太明显,不知道是我运行过程中出错了
发表于 2006-5-17 08:51 | 显示全部楼层
我运行了不过不明白这段是什么意思,
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

想请教一下,还有我觉得你选的参数是混沌的参数,可是出来的图并不是太明显,不知道是我运行过程中出错了
发表于 2006-5-18 14:24 | 显示全部楼层
shuijing :那一段程序是在求曲线和phi=2*pi/3这个截面的交点。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

快速回复 返回顶部 返回列表