Duffing系统的相图、频谱图、分岔图，poincare等

Rotor2014

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

2014-11-28 10:59 上传

点击文件名下载附件
相关程序以及图

学渣硬要撑学霸

等级不够 不能下载。。。。我也是醉了。。。

Rotor2014

学渣硬要撑学霸 发表于 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轴标注

复制代码

学渣硬要撑学霸

Rotor2014 发表于 2014-11-30 11:50
Duffing系统动力学方程计算程序
function dx=Duffing(t,x,flag,f); global f; %Duffing方程
a=0.25;

不是用Matlab编写的。。。用C语言编的 现在不会弄了。。。愁啊。。。怎么办咧 怎么办咧

x589

神啊  救救咱啊

x589

为什么吗不能下载啊  我急需哦

Rotor2014

x589 发表于 2014-11-30 19:57
为什么吗不能下载啊  我急需哦

你级别低吧  我已经粘帖到楼下了

frankfei111

感谢分享！

无水1324

感谢分享！

x589

多谢楼主  感激涕零哦

x589

咱中国，还是好人多啊！！！！

学渣硬要撑学霸

谁用C编的  指导下我呗。。。。

huadianlizhan

这程序好像不太好吧

Rotor2014

huadianlizhan 发表于 2014-12-11 22:15
这程序好像不太好吧

《机械振动系统的现代动态设计与分析》 这本书上的

wenglei907

顶一个   不错的程序