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

北瓜头 · 发表于 2017-3-8 16:27

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

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

x

我在做一个受简谐激励力作用下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.
不胜感激各位高手的帮助！！！

北瓜头 · 发表于 2017-3-8 16:36

求大神伸出援手

北瓜头 · 发表于 2017-3-8 21:56

求大神解救

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

function dx=Duffing(t,x,flag,f); global f; %Duffing方程
a=0.25;
b=-1;
c=1;
w=1;
dx=[x(2);-a*x(2)-b*x(1)-c*x(1)^3+f*cos(x(3));w];
%fftplot幅值谱
function [f,y]=fftplot(x,step)
[r,c] = size(x);
if r == 1
x = x.'; % make it a column
end;
[n,cc] = size(x);
m = 2^nextpow2(n);
y = fft(x,m);
y = y(1:n,:);
yy=abs(y)*step;
fid1=fopen('fft1.txt','w');
fid2=fopen('aft1.txt','w');
for kk=1:m/2+1;
f(kk)=(kk-1)/m/step;
fprintf(fid1,'%8.4f\n',f(kk));
fprintf(fid2,'%8.4f\n',y(kk));
end
fclose(fid1);
fclose(fid2);
y=yy(1:m/2+1);
%plot(f,y);
%xlabel('Frequency (Hz)');
%ylabel('Fourier Amplitude');
%title('Fouier Transform of the Data')
%求解、绘图：时域波形、相平面、幅值谱、Poincare图、分岔图
clc
global f;
f=0.1
tspan=[0:0.01*2*pi:2000];
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
[t,x]=ode45('Duffing',tspan,[1;1.5;0],options);
x_data=x(:,1);
x_dot=x(:,2);
%时域波形
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
plot(tspan,x_data,'LineWidth',1.5);
xlim([100 200])
%xlim([200 300])
xlabel('\itt','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('\itx','Fontsize',20,'Fontname','Times new roman','Rotation',0);%y轴标注
grid on
box on
%相平面图
figure(2)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
plot(x_data(1000:3000),x_dot(1000:3000),'LineWidth',1.5);
%plot(x_data(10:30000),x_dot(10:30000),'LineWidth',1.5);
xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
grid on
box on
%幅值谱
[freq,amp]=fftplot(x_data,0.01*2*pi);
figure(3)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
xlim([0 1])
ylim([0 200])
plot(freq(1:5000),amp(1:5000),'LineWidth',1.5);
xlabel('/Hz','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('mu','Fontsize',20,'Fontname','Times new roman');%y轴标注
grid on
box on
%poincare图
x(:,2)=mod(x(:,2),2*pi)-pi;
phi0=pi*1;
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
figure(4)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
plot(X1,X2,'.');
xlabel('\itx','Fontsize',20,'Fontname','Times new roman');%x轴标注
ylabel('d\itx/d\itt','Fontsize',20,'Fontname','Times new roman');%y轴标注
%分岔图
global f
figure(5)
hold on;
box on;
xlim([0,1]);
ylim([-2,2]);
N=80;
Tn=zeros(1,N);
options=odeset('RelTol',1e-9);
for f=0.05:.05:1;
x0=[0,1.5,0];
for n=1:60;
[t,x]=ode45('Duffing',[0,2*pi],x0,options);
x0=x(end,:);
end
for n=1:N;
[t,x]=ode45('Duffing',[0,pi*2],x0,options);
x0=x(end,:);
xd=x(:,1);
Tn(n)=max(xd);
end
f
plot(f,Tn,'b.','markersize',1);
pause(0.0001);
end
xlabel('\itf','Fontsize',16,'Fontname','Times new roman');%x轴标注
ylabel('\itx','Fontsize',16,'Fontname','Times new roman','Rotation',0);%y轴标注

复制代码

北瓜头 · 发表于 2017-3-9 20:31

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

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

Accompany · 发表于 2017-3-10 09:07

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

北瓜头 · 发表于 2017-3-10 14:15

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

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

xinhu9 · 发表于 2017-6-8 11:51

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

请问您找到了吗？

Kevin_ · 发表于 2021-3-26 19:43

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

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

账号		自动登录	找回密码
密码			我要加入

[分形与混沌] 受迫duffing振子的吸引盆程序求助

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