马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
这是书上的正确的图。
这是我编程出来得到的图。
两张图差太多了,不知道是什么原因,会仅仅只是参数初值的影响吗?
下面附一下我的程序。
------------------------------------------------------------
function dx=jeffliewen(t,x)
global w
m=32.1;
c=2100;
e=0.002;
k=2.5e7;
g=9.8;
s0=m*g/k;
omega=sqrt(k/m);
beta=0;
theta=0;
k1=1-beta/pi+sin(4*beta)/(4*pi);%法向刚度
k2=1-beta/pi+2*sin(2*beta)/(3*pi)-sin(4*beta)/(12*pi);%切向刚度
k3=0;%耦合刚度u
k11=k1*(cos(theta))^2+k2*(sin(theta))^2-k3*sin(2*theta);%x 方向刚度
k22=k1*(sin(theta))^2+k2*(cos(theta))^2+k3*sin(2*theta);%y方向刚度
k33=(k1-k2)*sin(theta)*cos(theta)+k3*cos(2*theta);%x,y方向耦合刚度(以上刚度来自中性轴理论)
dx=[x(2);
-x(2)*c/(m*omega)-x(1)*k11/k-x(3)*k33/k+1+cos(theta+beta)*e*w^2/(s0*omega^2);
x(4);
-x(4)*c/(m*omega)-x(3)*k22/k-x(1)*k33/k+sin(theta+beta)*e*w^2/(s0*omega^2)];
----------------------------------------------------------------------
clear;
clc;
hold on;
global w
f=10:3:1200
for i=1:length(f)
disp(f(i));
w=f(i);
T=2*pi;
x0=[0.1;0;0.1;0];
tspan=0:T/100:100*T;
[t,x]=ode45('jeffliewen',tspan,x0);
x0=x(end,:);
plot(f(i),x(4000:100:end,1),'markersize',5);%画分岔图
xlabel ('转速w');
ylabel ('位移x1');
end |