设为首页收藏本站联系我们

声振论坛

 找回密码
 我要加入

QQ登录

只需一步，快速开始

声振论坛 › 展示 ›振动理论› 非线性振动 › 查看内容

庞加莱截面图绘制matlab代码

2015-10-23 18:31| 发布者: aspen| 查看: 3372| 评论: 23|原作者: suffer|来自: 声振论坛

摘要: Poincare截面 　　在相空间中适当（要有利于观察系统的运动特征和变化，如截面不能与轨线相切，更不能包含轨线）选取一截面，在此截面上某一对共轭变量如x1和x.1取固定值，称此截面为Poincare截面，相空间的连续轨 ...
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];

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

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

  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 编辑 ]

%%% 另外用下面一个文件也可以实现的
  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;

  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 编辑 ]
收藏 分享 邀请
上一篇：关于最大Lyapunov指数计算的几个疑问下一篇：线性与非线性的理解
发表评论

最新评论

引用 suffer 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;

  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 2006-4-8 14:27
关键是实例麻烦在象我这种的人看不懂啊：（
呵呵，但还是顶，最好有点解释：）
引用 cdwxg 2006-4-9 00:03
我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了
引用 juan8909 2006-4-9 14:09
楼主，我刚接触混沌，请问这个程序能否利用原始的时间序列做出ＰＯＩＮＣＡＲＥ散点图啊？？
另外请教一下，ＭＡＴＬＡＢ中有非线性分析库吗？？
谢谢
引用 coldspring 2006-4-10 10:14
我感觉这个版块需要的高手要很多
最起码要有很多的注释才行
或者等到高手多了后再开？或者合并到哪个地方
现在开需要你们付出太多时间和精力来打里
因为如果只给例子没注释，对我们水平差的太痛苦了
引用 无水1324 2006-4-16 16:58
同意7楼的，最好来点注释。
还有闪频法画的相图跟ＰＯＩＮＣＡＲＥ散点图有什么区别？
这个好像有人问了，但是没有看到答案？
谢谢
引用 yfh2006 2006-4-26 16:52
非常感谢楼主的程序，但还是有些看不懂。怎么把实验数据导入，得到图形，望多多指教，谢谢了！
引用 ThomasLee 2006-5-6 21:16
这个才是正确的！
那个萝卜网站上的不正确阿！不要上当！
引用 shuijing 2006-5-12 08:33
版主你好！谢谢您辛苦了！我目前在解转子碰摩的微分方程，也做它的poincare映射图，程序还有看不懂的地方，希望能和你一起学习！
引用 shuijing 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
这段我不明白！真心希望请教你！
引用 shuijing 2006-5-15 16:48
不知道你用的是否是混沌参数么！
引用 shuijing 2006-5-17 08:48
我运行了不过不明白这段是什么意思，想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了
引用 shuijing 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

想请教一下，还有我觉得你选的参数是混沌的参数，可是出来的图并不是太明显，不知道是我运行过程中出错了
引用 toes 2006-5-18 14:24
shuijing :那一段程序是在求曲线和phi=2*pi/3这个截面的交点。
引用 百分百 2006-6-3 12:34
楼主,我想问一下:为什么把你的程序运用到我的系统上却运行不出来,我不知道为什么,特向你请教.我的系统为什么:

dx(1)/dt= -x(1)+tanh(x(1))+u*[ v*tanh(x(1))+2*tanh(x(2))+tanh(x(3)) ]
dx(2)/dt= -x(2)+tanh(x(2))+u*[ -2*tanh(x(1))+v*tanh(x(2)) ]
dx(3)/dt= -x(3)+tanh(x(3))+u*[ 4*tanh(x(1))-4*tanh(x(2))+tanh(x(3)) ]
其中u,v是参数. 当u=1, v=0.56
引用 cleverxue 2013-3-27 20:22
请问下实测数据怎么画庞加莱截面啊，相空间已经重构了，可是如何取截面啊
引用 天涯青驹 2013-4-7 16:17
引用 飞翼 2013-7-13 09:27
冒昧的问一句，这是定相位法？
引用 earth123654 2013-7-21 19:04
很不错，学习一下

查看全部评论(23)

相关分类

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-30 00:05 , Processed in 0.055748 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部