burif.m

这些东西相信对作转子动静碰磨得同学有帮助

clear all
tic

m=4       %      20;  %4.0;                 %kg
k=0.25*10^6       %    6*10^6; 
kc=6.0*10^7       %4*10^7                %Nm-1
c=1200;                %Nsm-1
f=0.2    %0.06;
e=0.00006    %0.05;          %m
jianxi=0.00015          %1*10^(-3);   %0.169;    %设置间隙值
w0s=k/m;              %w0平方
w0=sqrt(w0s);          %求w0
kcc=c/(2*m*w0);        %阻尼比
g=9.8; 

X0=[0.001;0.001;0;0]; %[0,0.1,2.1,0.2];
rate=[0.001:.01:7];

offset=1;
i=1;
for i=1:length(rate)
  
G=rate(i);
tspan=[0 2040*pi/rate(i)];
T=[0:pi/(1024*rate(i)):2040*pi/rate(i)];
  w=rate(i)*w0;
[T,Y]=ode45(@TESTXdot,T,X0',[],m,k,kc,c,f,e,rate(i),jianxi,w0s,w0,g);
t=T';

        cutindexes=find(T>2000*pi/rate(i));
        T=T(cutindexes);
        Ytmp1=Y(cutindexes,1);
        Ytmp2=Y(cutindexes,2);
        Ytmp1(:,2)=Ytmp2(:,1);
        Y=Ytmp1;
        t=t(cutindexes);

tmp=length(Y);
m(offset:(tmp+offset-1),1)=g(i);
m(offset:(tmp+offset-1),2)=Y(:,2);
offset=tmp+offset;
end
figure
        plot(m(:,1),m(:,2),'k.','markersize',4)
        title('Bifurcation Diagramm')
        xlabel('s')
        ylabel('Y')
toc;