[未答复] 正确绘制庞加莱截面,取点间隔要为一个激励周期!

[复制链接]
一蓑烟_MCQ7L 发表于 2021-9-29 00:32:49
本帖最后由 一蓑烟_MCQ7L 于 2021-9-29 23:33 编辑
  1. <div class="blockcode"><blockquote>%% % © Copyright 一蓑烟雨任平生
  2. %绘制受驱单摆庞加莱截面 Poincare Portrait
  3. clc;clear all; close all;
  4. %% 参数
  5. Q=4;       %阻尼无量纲参数  Q=1/2β=mω0/γ
  6. omiga=1; %驱动频率无量纲参数  Ω=ν/ω0
  7. f=1.4;  %驱动力无量纲参数 f=F/mg
  8. %% 数值求解
  9. T=2*pi/omiga; %初始条件 ZQ:周期
  10. tspan=0:T/100:1000*T;
  11. y0=[0 0.1];
  12. [t,y]=ode45(@(t,y) pendulum_nolinear(t,y,Q,f,omiga),tspan,y0);%解方程 在一个周期ZQ中取100个点
  13. %% 画位移曲线
  14. figure
  15. plot(t(1000:3000),y(1000:3000,1));title('位移曲线');
  16. %% 画相图   
  17. figure
  18. %subplot(2,1,1)
  19. plot(y(1000:end,1),y(1000:end,2));hold on
  20. xtick=get(gca,'XTick');
  21. ytick=get(gca,'YTick');
  22. %f=num2str(f);s2='_.jpg';s11=strcat('f=',f,s2);
  23. %saveas(gcf,s11);
  24. a=size(y,1);
  25. %% 画庞加菜截面图
  26. %subplot(2,1,2)
  27. phi=0.5*pi/omiga; %相角(截面位置)
  28. T=2*pi/omiga;   %时间周期
  29. dt=t(2)-t(1);   %计算步长间隔
  30. dn=round (T/dt);%一个外激励周期内对应的间隔数
  31. nn=size(t,1);
  32. j0=dn*500+round(phi/dt);
  33. j=j0:dn:nn;
  34. plot(y(j,1),y(j,2),'k.','markersize',20);

  35. function dydt=pendulum_linear(t,y,Q,f,omiga)  
  36. dydt=[y(2);-y(2)/Q-y(1)+f*cos(omiga*t)]; %驱动线性摆方程:y(1)表示角度,y(2)表示角速度
  37. end

  38. function dydt=pendulum_nolinear(t,y,Q,f,omiga)  
  39. dydt=[y(2);-y(2)/Q-sin(y(1))+f*cos(omiga*t)]; %驱动非线性摆方程:y(1)表示角度,y(2)表示角速度
  40. end
复制代码



参数:Q=2;omiga=2/3; F=1.21; tspan=0:T/100:100000*T;y0=[0.2 0];

参数:Q=2;omiga=2/3; F=1.21; tspan=0:T/100:100000*T;y0=[0.2 0];
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表