怎么样才能绘制出最准确、最精确的分岔图？？

liliangbiao

回复 30 # 410610288 的帖子

要是方便的话，吧你看的资料给我发一份！
lihaimin1955@yahoo.com.cn

410610288

回复 31 # liliangbiao 的帖子

已发送，期待liliangbiao先生的及时回复与宝贵建议。

paopaotaiqiang

学习下

a139

很详细，顶起来

wangu

向高手学习

lihaitao123

回复 1 # yufeiyfyf 的帖子

为什么要   
          ;？
     rosser  里  x0=[0.5,0.4,0.3]'  y0=[0.5;0.3;0.2];
    而x0与y0相同

lihaitao123

回复 1 # yufeiyfyf 的帖子

楼主，我用这个程序，套在我的系统中，出现同样情况，
??? Subscript indices must either be real positive integers or logicals.
不知楼主解决了没有？？

boom

学长。能不能帮我查看一下程序的错误，是关于碰摩转子的分岔图的，万分感谢！
function zdot=pengmo(tt,z,flag,w)
zdot=zeros(8,1);


m1=4.0;%转子在轴承处集中质量kg
m2=32.1;%转子在圆盘处等效的集中质量kg
b1=0.00005;%偏心量m
R=0.025;%轴承半径m
L=0.012;%轴承长度m
c=0.00011;%平均油膜厚度m
u=0.018;%润滑油黏度N*s/m
%w=900;%转子转动角速度
f=0.1;
c1=1050;%转子在轴承出阻尼系数N*S/m
c2=2100;%转子圆盘轴承系数N*s/m
k=2.5*10^7;%静子刚度N/m
g=9.8;%N/kg
delta=0.005;%m
%delta=0.000016;
k1= 3.5*10^6;
G1=g/(c*w^2);
b=b1/c;

tau=w*tt;
P=m2*g/2;%转子圆盘重量的一半
s=(u*w*R*L*((R/c)^2)*(L/(2*R))^2)/P; % s为Sommerfeld修正系数，c平均油膜厚度
e=sqrt(z(5)^2+z(7)^2);

if e>=delta
Px=-c*((e-delta)/e)*k1*(z(5)-f*z(7));
Py=-c*((e-delta)/e)*k1*(z(5)*f+z(7));

else
Px=0;
Py=0;
end
   

alpha=atan((z(3)+2*z(2))/(z(1)-2*z(4)))-pi/2*sin((z(3)+2*z(2))/(z(1)-2*z(4)))-pi/2*sin(z(3)+2*z(2));
S=(z(1)*cos(alpha)+z(3)*sin(alpha))/(1-abs((z(1)*cos(alpha)+z(3)*sin(alpha))^2));
G=(2/abs(sqrt(1-abs(z(1)^2)-abs(z(3)^2))))*(pi/2+atan((z(3)*cos(alpha)-z(1)*sin(alpha))/abs(sqrt(1-abs(z(1)^2)-abs(z(3)^2)))));
V=(2+(z(3)*cos(alpha)-z(1)*sin(alpha))*G)/abs((1-abs(z(1)^2)-abs(z(3)^2)));

fx=-abs((abs((z(1)-2*z(4))^2)+abs((z(3)+2*z(2))^2))^(1/2))/(1-abs(z(1)^2)-abs(z(3)^2))*(3*z(1)*V-G*sin(alpha)-2*S*cos(alpha));
fy=-abs((abs((z(1)-2*z(4))^2)+abs((z(3)+2*z(2))^2))^(1/2))/(1-abs(z(1)^2)-abs(z(3)^2))*(3*z(3)*V-G*cos(alpha)-2*S*sin(alpha));%油膜力

zdot(1)=z(2);
zdot(2)=-c1*z(2)/(w*m1)-k*(z(1)-z(5))/(w^2*m1)+s*P*fx/(w^2*m1*c)+b*cos(tau);
zdot(3)=z(4);
zdot(4)=-c1*z(4)/(w*m1)-k*(z(3)-z(7))/(w^2*m1)+s*P*fy/(w^2*m1*c)+b*sin(tau)-G1;
zdot(5)=z(6);
zdot(6)=-c2*z(6)/(w*m2)-2*k*(z(5)-z(1))/(w^2*m2)+Px/(c*m2*w^2);

clc
clear
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=0:25:1000;           %转子的转速是以rad/s为单位的转速


%位移、速度、转角、转角速度的初值
z1=1e-005;   z2=1e-005;   z3=-1e-5;     z4=3e-005;      z5=2e-005;   z6=-3e-5;    z7=-1e-005;   z8=1e-005;     

for  n=1:length(w);
      T=2*pi;
      ts=[0:T/200:T*100];                          % 时间区间/w(n)
      z0=[z1;z2;z3;z4;z5;z6;z7;z8;];       % xo包含的三十个向量           
      [tt,z]=ode45('pengmo',ts,z0,[],w(n));
   
      figure(1)
      plot(w(n),z(100:100:end,1),'*');
      xlabel('\fontsize{18}\omega');
      ylabel('\fontsize{18}x');grid
      hold on
      
  end
toc

      





zdot(7)=z(8);
zdot(8)=-c2*z(8)/(w*m2)-2*k*(z(7)-z(3))/(w^2*m2)+Py/(c*m2*w^2)-G1;

分岔图

这个是分岔图？？？可以参考下 prof。chip的论文，定义的分插图。这个由一定的争议。

挑战振动

tspan=[0 200];zhang=[zhang;fradata1],时间步长太小运行会出现r里面有负的，那个；你打成：了。改了这两处就可以运行了！估计你用不上了，我最近学习分岔的，刚看到