Rotor2014 发表于 2014-11-28 11:00

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

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

自己依据书上的例子做的Duffing系统的相图、频谱图、分岔图，poincare等。对于相关的研究有一定的参考作用。


学渣硬要撑学霸 发表于 2014-11-30 11:39

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

Rotor2014 发表于 2014-11-30 11:50

本帖最后由 牛小贱 于 2014-11-30 19:50 编辑

学渣硬要撑学霸 发表于 2014-11-30 11:39
等级不够 不能下载。。。。我也是醉了。。。
Duffing系统动力学方程计算程序：function dx=Duffing(t,x,flag,f); global f; %Duffing方程
a=0.25;
b=-1;
c=1;
w=1;
dx=;
%fftplot幅值谱
function =fftplot(x,step)
= size(x);
if r == 1
      x = x.';   % make it a column
end;
= 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=;
options=odeset('RelTol',1e-5,'AbsTol',1e-6);
=ode45('Duffing',tspan,,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()
%xlim()
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
%幅值谱
=fftplot(x_data,0.01*2*pi);
figure(3)
axes('FontSize',16,'FontName','Times New Roman','LineWidth',1.5);
hold on;
xlim()
ylim()
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;
    =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(,,1);
      X1(k)=polyval(Q,(k-1)*2*pi-phi0);
         Q=polyfit(,,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();
ylim([-2,2]);
N=80;
Tn=zeros(1,N);
options=odeset('RelTol',1e-9);
for f=0.05:.05:1;
x0=;
for n=1:60;
=ode45('Duffing',,x0,options);
x0=x(end,:);
end
for n=1:N;
=ode45('Duffing',,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轴标注



学渣硬要撑学霸 发表于 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语言编的 现在不会弄了。。。愁啊。。。怎么办咧 怎么办咧

x589 发表于 2014-11-30 19:57

神啊救救咱啊

x589 发表于 2014-11-30 19:57

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

Rotor2014 发表于 2014-11-30 20:41

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

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

frankfei111 发表于 2014-11-30 20:50

感谢分享！

无水1324 发表于 2014-12-1 10:44

感谢分享！

x589 发表于 2014-12-1 15:04

多谢楼主感激涕零哦

x589 发表于 2014-12-1 15:10

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

学渣硬要撑学霸 发表于 2014-12-1 20:44

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

huadianlizhan 发表于 2014-12-11 22:15

这程序好像不太好吧

Rotor2014 发表于 2014-12-11 22:48

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

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

wenglei907 发表于 2014-12-22 11:11

顶一个   不错的程序
页: [1] 2
查看完整版本: Duffing系统的相图、频谱图、分岔图，poincare等