受迫duffing振子的吸引盆程序求助

北瓜头

我在做一个受简谐激励力作用下duffing方程的研究，需要绘制其吸引域。目前查书和文献了解了一下吸引域的绘制方法，但是对编程还是没有头绪，能否提供一个duffing方程吸引域绘制的程序或者类似的程序让我参考一下？
我所做的duffing方程为
dx=[x(2);
        F*cos(Omega*t)-(0.02*x(2)+x(1)+0.54/(1-0.54)*x(1)*(1-1/(sqrt(0.54^2+x(1).^2))))]
F=0.01; Omega=0.54.
不胜感激各位高手的帮助！！！

北瓜头

求大神伸出援手

北瓜头

求大神解救

Catsayer

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轴标注

复制代码

北瓜头

Catsayer 发表于 2017-3-9 11:11

谢谢你，不过这些程序我都编出来了，能否提供一下吸引域的编程参考一下？

Accompany

这方面论文不少  代码没见过

北瓜头

Accompany 发表于 2017-3-10 09:07
这方面论文不少  代码没见过

是这样的，我查过不少论文，不过由于编程水平有限，没啥头绪。所以想找一两个代码来参考参考。

xinhu9

北瓜头 发表于 2017-3-10 14:15
是这样的，我查过不少论文，不过由于编程水平有限，没啥头绪。所以想找一两个代码来参考参考。

请问您找到了吗？

Kevin_

北瓜头 发表于 2017-3-10 14:15
是这样的，我查过不少论文，不过由于编程水平有限，没啥头绪。所以想找一两个代码来参考参考。

请问能否交流一下 吸引域 的画法？