%非线性弹性支承转子系统
function dy=diffequat1(t,y)
global theta e kz k w W r2 M m g D beta;
e=0.0002;
theta=1.415263;
r2=1000.00;
w=0.25;
W=500;
kz=19304448;
M=50;
m=10;
D=0.008;
g=9.800;
k=1930444.8;
beta=1;
dy=[y(2);
e*kz/M/g*cos(t-theta)-r2*kz*w^2/M/g*sin(theta)-kz/M/W^2*cos(theta)-2*D*(kz/M)^0.5/W*(y(2)-y(4))-kz/M/W^2*(y(1)-y(3));
y(4);
-r2*kz*w^2/W^2/m/g*sin(theta)-kz/m/W^2*cos(theta)+2*D*M*(kz/M)^0.5/W/m*(y(2)-y(4))+kz/m/W^2*(y(1)-y(3))-2*D*M*(kz/M)^0.5/W/m*y(4)-k/m/W^2*y(3)-beta*k/m/W^2*y(3)^3;
y(6);
e*kz/M/g*sin(t-theta)-r2*kz*w^2/M/g*cos(theta)-kz/M/W^2*cos(theta)-2*D*(kz/M)^0.5/W*(y(6)-y(8))-kz/M/W^2*(y(5)-y(7));
y(8);
-r2*kz*w^2/W^2/m/g*cos(theta)-kz/m/W^2*sin(theta)+2*D*M*(kz/M)^0.5/W/m*(y(6)-y(8))+kz/m/W^2*(y(5)-y(7))-2*D*M*(kz/M)^0.5/W/m*y(8)-k/m/W^2*y(7)-beta*k/m/W^2*y(7)^3];
绘制分岔图程序:
clear all;
clc;
W=100:10:1300;
for h=1:length(W)
T=2pi;
[t,x]=ode45(‘diffequat1’,[1:T/100:1000*T],[11 1 1 1 1 1 1]),[],W(h));
plot(W(h),x(80000:100:end,1),’k.’);
hold on;
end
绘制的图出来感觉是不对的。