|
楼主 |
发表于 2014-11-30 11:50
|
显示全部楼层
本帖最后由 牛小贱 于 2014-11-30 19:50 编辑
Duffing系统动力学方程计算程序:- 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轴标注
复制代码
|
|