做过转子碰摩的前辈请进来帮忙看下程序
%非线性刚度转子系统碰摩的混沌行为%罗跃纲.2002.东北大学学报.vol.23，no.9
function dx=rubbing(t,x);
global gamma;
dx=zeros(4,1);
%parameter
xi=0.12;
alpha=0.5;
beta=3;
f=0.12;
delta=0.16;
b=0.136;
g=9.8;
f0=25;
w0=2*pi*f0;
e=sqrt(x(1)^2+x(3)^2);
H=1/2*(sign(abs(e-1))+sign(e-1));%判断是否发生碰磨
%dimensionless equations
dx=[x(2);
-2*xi*x(2)-x(1)-alpha*(x(1)^2+x(3)^2)*x(1)-beta*(1-1/e)*(x(1)-f*x(3))*H+b*gamma^2*cos(gamma*t);
x(4);
-2*xi*x(4)-x(3)-alpha*(x(1)^2+x(3)^2)*x(3)-beta*(1-1/e)*(f*x(1)+x(3))*H+b*gamma^2*sin(gamma*t)-g/(delta*w0^2)];
主程序
clear all;
clc;
global gamma;
range=;
period=2*pi;
k=0;
YY1=[];
step=period/512;
for gamma=range
x0=;
gamma
k=k+1;
% discard the first *** periodic data;
tspan=;
options=odeset('RelTol',10^-3,'AbsTol',10^-5);
=ode45('rubbing',tspan,x0,options);
x0=x(end,:);
j=1;
for i=2000:2200
tspan=;
options=odeset('RelTol',10^-3,'AbsTol',10^-5);
=ode45('rubbing',tspan,x0,options);
YY1(k,j)=x(end,1); % get the omega data from every period end
j=j+1;
x0=x(end,:);
end
end
bifdata=YY1(:,end-51:end);
plot(range,bifdata,'k.','markersize',1);
xlabel('频率比w/w0');ylabel('x');title('随频率比变化的分岔图');
我得到的图是几根线，跟原文中的分岔图差的太多了，不知道错在了哪里？请帮忙看下
原文分岔图 学习了，正想涉及非线性
页:
[1]