声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5787|回复: 28

[稳定性与分岔] Duffing系统的相图、频谱图、分岔图,poincare等

[复制链接]
发表于 2014-11-28 11:00 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 Rotor2014 于 2014-11-28 11:00 编辑

自己依据书上的例子做的Duffing系统的相图、频谱图、分岔图,poincare等。对于相关的研究有一定的参考作用。
Duffing系统动力学方程.zip (341.92 KB, 下载次数: 170)

点评

赞成: 5.0
赞成: 5
  发表于 2014-11-30 19:49

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2014-11-30 11:39 | 显示全部楼层
等级不够 不能下载。。。。我也是醉了。。。
 楼主| 发表于 2014-11-30 11:50 | 显示全部楼层
本帖最后由 牛小贱 于 2014-11-30 19:50 编辑
学渣硬要撑学霸 发表于 2014-11-30 11:39
等级不够 不能下载。。。。我也是醉了。。。

Duffing系统动力学方程计算程序:
  1. function dx=Duffing(t,x,flag,f); global f; %Duffing方程
  2. a=0.25;
  3. b=-1;
  4. c=1;
  5. w=1;
  6. dx=[x(2);-a*x(2)-b*x(1)-c*x(1)^3+f*cos(x(3));w];
  7. %fftplot幅值谱
  8. function [f,y]=fftplot(x,step)
  9. [r,c] = size(x);
  10. if r == 1
  11.         x = x.';   % make it a column
  12. end;
  13. [n,cc] = size(x);
  14. m = 2^nextpow2(n);
  15. y = fft(x,m);
  16. y = y(1:n,:);
  17. yy=abs(y)*step;
  18. fid1=fopen('fft1.txt','w');
  19. fid2=fopen('aft1.txt','w');
  20. for kk=1:m/2+1;
  21.    f(kk)=(kk-1)/m/step;
  22.    fprintf(fid1,'%8.4f\n',f(kk));
  23.    fprintf(fid2,'%8.4f\n',y(kk));
  24. end
  25. fclose(fid1);
  26. fclose(fid2);
  27. y=yy(1:m/2+1);
  28. %plot(f,y);
  29. %xlabel('Frequency (Hz)');
  30. %ylabel('Fourier Amplitude');
  31. %title('Fouier Transform of the Data')
  32. %求解、绘图:时域波形、相平面、幅值谱、Poincare图、分岔图
  33. clc
  34. global f;
  35. f=0.1
  36. tspan=[0:0.01*2*pi:2000];
  37. options=odeset('RelTol',1e-5,'AbsTol',1e-6);
  38. [t,x]=ode45('Duffing',tspan,[1;1.5;0],options);
  39. x_data=x(:,1);
  40. x_dot=x(:,2);
  41. %时域波形
  42. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  43. hold on;
  44. plot(tspan,x_data,'LineWidth',1.5);
  45. xlim([100 200])
  46. %xlim([200 300])
  47. xlabel('\itt','Fontsize',20,'Fontname','Times new roman');%x轴标注
  48. ylabel('\itx','Fontsize',20,'Fontname','Times new roman','Rotation',0);%y轴标注
  49. grid on
  50. box on
  51. %相平面图
  52. figure(2)
  53. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  54. hold on;
  55. plot(x_data(1000:3000),x_dot(1000:3000),'LineWidth',1.5);
  56. %plot(x_data(10:30000),x_dot(10:30000),'LineWidth',1.5);
  57. xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  58. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  59. grid on
  60. box on
  61. %幅值谱
  62. [freq,amp]=fftplot(x_data,0.01*2*pi);
  63. figure(3)
  64. axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  65. hold on;
  66. xlim([0 1])
  67. ylim([0 200])
  68. plot(freq(1:5000),amp(1:5000),'LineWidth',1.5);
  69. xlabel('/Hz','Fontsize',20,'Fontname','Times new roman');%x轴标注
  70. ylabel('mu','Fontsize',20,'Fontname','Times new roman');%y轴标注
  71. grid on
  72. box on
  73. %poincare图
  74. x(:,2)=mod(x(:,2),2*pi)-pi;
  75. phi0=pi*1;
  76. for k=1:round(max(x(:,3))/2/pi);
  77.     d=x(:,3)-(k-1)*2*pi-phi0;
  78.     [P,K]=sort(abs(d));
  79.     x1l=x(K(1),1);
  80.     x1r=x(K(2),1);
  81.     x2l=x(K(1),2);
  82.     x2r=x(K(2),2);
  83.     x3l=x(K(1),3);
  84.     x3r=x(K(2),3);
  85.     if abs(P(1))+abs(P(2))<3e-16;
  86.         X1(k)=x1l;
  87.         X2(k)=x2l;
  88.     else
  89.         Q=polyfit([x3l,x3r],[x1l,x1r],1);
  90.         X1(k)=polyval(Q,(k-1)*2*pi-phi0);
  91.          Q=polyfit([x3l,x3r],[x2l,x2r],1);
  92.         X2(k)=polyval(Q,(k-1)*2*pi-phi0);
  93.     end
  94. end
  95.     figure(4)
  96.    axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
  97.    plot(X1,X2,'.');
  98.   xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
  99. ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
  100. %分岔图
  101. global f
  102. figure(5)
  103. hold on;
  104. box on;
  105. xlim([0,1]);
  106. ylim([-2,2]);
  107. N=80;
  108. Tn=zeros(1,N);
  109. options=odeset('RelTol',1e-9);
  110. for f=0.05:.05:1;
  111. x0=[0,1.5,0];
  112. for n=1:60;
  113. [t,x]=ode45('Duffing',[0,2*pi],x0,options);
  114. x0=x(end,:);
  115. end
  116. for n=1:N;
  117. [t,x]=ode45('Duffing',[0,pi*2],x0,options);
  118. x0=x(end,:);
  119. xd=x(:,1);
  120. Tn(n)=max(xd);
  121. end
  122. f
  123. plot(f,Tn,'b.','markersize',1);
  124. pause(0.0001);
  125. end
  126. xlabel('\itf','Fontsize',16,'Fontname','Times new roman');%x轴标注
  127. ylabel('\itx','Fontsize',16,'Fontname','Times new roman','Rotation',0);%y轴标注
复制代码




点评

赞成: 5.0
赞成: 5
  发表于 2014-11-30 19:49
发表于 2014-11-30 17:37 | 显示全部楼层
Rotor2014 发表于 2014-11-30 11:50
Duffing系统动力学方程计算程序
function dx=Duffing(t,x,flag,f); global f; %Duffing方程
a=0.25;

不是用Matlab编写的。。。用C语言编的 现在不会弄了。。。愁啊。。。怎么办咧 怎么办咧
发表于 2014-11-30 19:57 | 显示全部楼层
神啊  救救咱啊
发表于 2014-11-30 19:57 | 显示全部楼层
为什么吗不能下载啊  我急需哦
 楼主| 发表于 2014-11-30 20:41 | 显示全部楼层
x589 发表于 2014-11-30 19:57
为什么吗不能下载啊  我急需哦

你级别低吧  我已经粘帖到楼下了
发表于 2014-11-30 20:50 | 显示全部楼层
感谢分享!
发表于 2014-12-1 10:44 | 显示全部楼层
感谢分享!
发表于 2014-12-1 15:04 | 显示全部楼层
多谢楼主  感激涕零哦
发表于 2014-12-1 15:10 | 显示全部楼层
咱中国,还是好人多啊!!!!
发表于 2014-12-1 20:44 | 显示全部楼层
谁用C编的  指导下我呗。。。。
发表于 2014-12-11 22:15 | 显示全部楼层
这程序好像不太好吧
 楼主| 发表于 2014-12-11 22:48 | 显示全部楼层
huadianlizhan 发表于 2014-12-11 22:15
这程序好像不太好吧

《机械振动系统的现代动态设计与分析》 这本书上的
发表于 2014-12-22 11:11 | 显示全部楼层
顶一个   不错的程序
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-5 13:50 , Processed in 0.068844 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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