转子碰摩Matlab'程序及分岔图

Rotor2014

自己学习用到的碰摩转子的matlab程序及分岔图程序，可以通过修改实现自己的用途。有没有人做裂纹转子的啊，自己编的程序画出来的分岔图不对，会的交流一下，谢谢

function zdot=rubbing(t,z)

global w;

global u;

global beta;

zdot=zeros(4,1);

%parameter

kesi=0.8372;

beta=12.25;

delta=1e-4;

u=4e-5;

U=u/delta;

g=9.8;

f0=173.2;

f=0.2;

G=g/(delta*(f0)^2);

E=sqrt(z(1)^2+z(3)^2);

%dimensionless equations

zdot=[z(2);

   -2*kesi*z(2)-z(1)-(1-1/E)*(beta^2)*(z(1)-f*z(3))*1/2*(sign(abs(E-1))+sign(E-1))+U*w^2*cos(w*t);

   z(4);

   -2*kesi*z(4)-z(3)-(1-1/E)*(beta^2)*(f*z(1)+z(3))*1/2*(sign(abs(E-1))+sign(E-1))+U*w^2*sin(w*t)-G];

clear all;

clc;

global w;

global u;

%global E;

%global beta;

u=0.00007;

f0=250;

range=[2.5:0.005:4];

k=1;

YY=[];%初始化YY

z0=[0 0.001 1.27e-3 0.001];%设置初值

for w=range

  period=2*pi/w;

  step=period/512;

    w

   %j=j+1;

    % discard the first *** periodic data;

    %除去前面2000个周期的数据，并将最后的结果作为下一次积分的初值

    tspan=[0:step:2000*period];

    %options=odeset('RelTol',10^-3,'AbsTol',10^-5);

   [t,Y]=ode45('rubbing',tspan,z0);

   z0=Y(end,:)

   %j=1;

   %for i=2000:2200

       tspan=[2000*period:step:2100*period];

       %options=odeset('RelTol',10^-3,'AbsTol',10^-5);

       [t,Y]=ode45('rubbing',tspan,z0);

       YY(k,:)=Y(1:512:end,1);   % getthe omega data from every period end

       %j=j+1;               

       z0=Y(end,:);

       k=k+1;

   %end

end

bifdata=YY(:,end-51:end);

plot(range,bifdata,'k.','LineWidth',1);

xlabel('频率比w/w0');ylabel('x');title('随频率比变化的分岔图');

%xlabel('刚度比kc/k0');ylabel('x');title('随刚度比变化的分岔图');

无水1324

是什么问题呢？

Rotor2014

无水1324 发表于 2014-12-1 12:52
是什么问题呢？

没问题，是网上找到的程序，分享给大家看看的

wenglei907

顶一个   不错的程序

liyurun

mark一下，慢慢学习

1713573225

学习了好好研究下

Overdosed

先学习一下。
论坛里关于分岔的新帖子不多了，是都不搞这方面的研究了么。

Triste

程序里的z是什么  是数据文件吗  能否给一个数字文件的例子  不然运行不了

sizhiyuan2006

谢谢楼主分享