谢谢无水1324 的回复
分岔图还是没做好，轨迹图和映射图也没做出来，到底怎么回事啊:'( ，可能还是前面方程求解有问题，还请各位多多指点 恩，是不是那个积分进行不下去，你要仔细检查一下程序，或者调整一下参数。回复 17楼 无水1324 的帖子
是的，在做轨迹图和时程图是就发现警告说积分最小误差不能满足，程序应该没问题，估计是里面的参数有问题，参数不知道怎么调好，从去年折腾到今年 w=400:5:600;for i=1:length(w)
disp(w(i));
T=2*pi/w(i); %既然你是非自治系统，里面是sin（t）、cos（t），所以T=2*pi
tspan=;
x0=;%初值
=ode45(@fangcheng,tspan,x0,[ ],w(i));
%y0=x(end,:);%这里不明白 %应该是x0=y（end,:），是指每计算一个周期就把该周期的最后一次迭代数据作为初值带入下一个周期，你再试试
plot(w(i),y(4000:100:end,1),'markersize',2);%这里应该是频率随位移变化的图形 %我计算了一下，用龙格库塔法解得y值只有几十个，这一步是画不出来的，估计你函数有问题，你再检查一下
hold on;
xlabel('频率');
ylabel('位移');
end
回复 19楼 hsfy919 的帖子
谢谢fsfy919的指点，不过问题还是多多，图形出不来，以下是出现的问题，还请多多赐教400
Warning: Failure at t=1.760105e+001.Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-014) at time t.
> In ode45 at 371
In myxiugai0316 at 14
??? Error using ==> plot
Vectors must be the same lengths.
w=400:5:600;
for i=1:length(w)
disp(w(i));
T=2*pi;
tspan=;
x0=;
=ode45(@fangcheng,tspan,x0,[ ],w(i));
x0=y(end,:);
plot(w(i),y(4000:100:end,1),'markersize',2);
hold on;
xlabel('频率');
ylabel('位移');
end
回复 20楼 yangxiaokang 的帖子
帮你调试了一下，程序通了，但不知道是不是你要的%%%%%%%%%%%%%%%%%%%函数文件
function dx=fangcheng_fun(t,x)
%参数
%global w
w=760;
kc=3.5e6;
m1=4;%轴承处集中质量
m2=32.1;%圆盘处集中质量
R=0.025;%轴承半径
L=0.012;%轴承长度
c=0.00011;%油膜厚度
mu=0.018;%油膜黏度
f=0.1;%摩擦系数
c1=1050;%轴承阻尼
c2=2100;%圆盘阻尼
k=2.5e7;%弹性轴刚度
r=0.00005;%圆盘质量偏心距
delta=0.000016;%转子与定子间隙
g=9.8;
d=r/c;
G1=g/(c*w^2);
G2=g/(c*w^2);
e=sqrt(x(5)^2+x(7)^2); %圆盘中心的径向位移
%判断是否发生碰摩
H=1/2*(sign(abs(e-delta))+sign(e-delta)); %判断是否发生碰摩，sign正数返回1，负数返回-1
%碰摩力
Px=-c*(e-delta)*kc/e*(x(5)-f*x(7))*H;
Py=-c*(e-delta)*kc/e*(f*x(5)+x(7))*H;
%非线性油膜力
s=mu*w*R*L*(R/c)^2*(L/(2*R))^2;%短轴承油膜力的sommerfeld数
a=atan((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign((x(3)+2*x(2))/(x(1)-2*x(4)))-pi/2*sign(x(3)+2*x(2));
G=2/sqrt(1-x(1)^2-x(3)^2)*(pi/2+atan((x(3)*cos(a)-x(1)*sin(a))/sqrt(1-x(1)^2-x(3)^2)));
S=(x(1)*cos(a)+x(3)*sin(a))/(1-(x(1)*cos(a)+x(3)*sin(a))^2);
V=(2+(x(3)*cos(a)-x(1)*sin(a))*G)/(1-x(1)^2-x(3)^2);
fx=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(1)*V-sin(a)*G-2*cos(a)*S);
fy=sqrt((x(1)-2*x(4))^2+(x(3)+2*x(2))^2)/(1-x(1)^2-x(3)^2)*(3*x(3)*V+cos(a)*G-2*sin(a)*S);
%方程组
dx=[x(2);
-c1/(w*m1)*x(2)-k/(w^2*m1)*(x(1)-x(5))+s/(c*w^2*m1)*fx;
x(4);
-c1/(w*m1)*x(4)-k/(w^2*m1)*(x(3)-x(7))+s/(c*w^2*m1)*fy-G1;
x(6);
-c2/(w*m2)*x(6)-2*k/(w^2*m2)*(x(5)-x(1))+Px/(c*w^2*m2)+d*cos(t);
x(8);
-c2/(w*m2)*x(8)-2*k/(w^2*m2)*(x(7)-x(3))+Py/(c*w^2*m2)+d*sin(t)-G2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%下面是主程序
clear
clc
hold on;
global w f
f=400:5:600;
x0=;
for i=1:length(f)
disp(f(i));
w=f(i);
T=2*pi;
tspan=0:T/100:100*T;
options = odeset('RelTol', 1e-1,'AbsTol',1e-1);
=ode45('fangcheng_fun',tspan,x0,options,w);
x0=x(end,:);
plot(f(i),x(8000:100:end,1),'markersize',2);
xlabel('频率');
ylabel('位移');
end
%%%%%%%%%%顺便问一下，你是搞转子那部分的
回复 20楼 yangxiaokang 的帖子
哦，对不起，我又试了一下，主程序中options项还可以设置更精确些'RelTol', 1e-2,'AbsTol',1e-6
当然结果越精确越好了，呵呵
另外plot(f(i),x(8000:100:end,1),'markersize',2);
改为plot(f(i),x(4000:100:end,1),'markersize',2); %%%%%%尽量点取多点
回复 22楼 hsfy919 的帖子
能否说明一下，程序那里做了修改哈哈！回复 21楼 hsfy919 的帖子
非常感谢hsfy919的热心帮助，小弟在此拜谢了。我是做转子碰摩的，我们这边就我一个在做，没人商量探讨，一个人搞了很长时间也没结果，几乎灰心了。回复 24楼 yangxiaokang 的帖子
另外还要感谢无水大哥一如既往的关注和帮助，指点了很多，谢谢怎么我运算的还是有问题，x有虚部
我照着hsfy919修改的程序，运算了几遍，警告说x虚部被忽略，得到的图形也有问题。分叉图计算有些难，可不可以先做时程图和轴心轨迹图，还请各位帮忙:handshake回复 26楼 yangxiaokang 的帖子
那些警告不用理会，意思是在计算中得到了共轭复根，但在画图中用实数轴不能表示，故忽略，在实数域不影响绘图，况且复数域并没有研究价值。%%%%%%%%%%%%%%%%%这是时间响应曲线
clear
clc
tspan=; %%%%%%%%%%%%%时间范围可以自己改变
x0=;
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
=ode45('fangcheng_fun',tspan,x0,options);
plot(t(10000:end),x(10000:end,1)) %%%%%%%%%%%%%%我没见你的图，也不清楚你的x1表示什么，所以先取了x1，你可以更改
%%%%%%%%%%%%%%%%%%%%%%%%%函数文件见21楼，不变
我做的图如下
非常感谢hsfy919帮助
按你说的可以计算时程图，跟你计算出来图形一样，不过时间范围变大如tspan=;图形就有问题了%%%%%%%%%%%%%%%%%这是时间响应曲线
clear
clc
tspan=; %%%%%%%%%%%%%时间范围可以自己改变
x0=;
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
=ode45('fangcheng_fun',tspan,x0,options);
plot(t(10000:end),x(10000:end,1)) %%%%%%%%%%%%%%我没见你的图，也不清楚你的x1表示什么，所以先取了x1，你可以更改。我这里的x1表示左轴颈中心横坐标的位移，x3是纵坐标位移；x5是圆盘中心横坐标的位移，x7是纵坐标位移。可以看看我前面的模型。
%%%%%%%%%%%%以下是我的计算程序，方程程序没变
clear
clc
tspan=;
x0=;
options = odeset('RelTol', 1e-2,'AbsTol',1e-6);
=ode45('fangcheng_fun',tspan,x0,options);
figure(1)%左轴颈x1的时程图
plot(t(4000:end),x(4000:end,1));title('时程图'); xlabel('t');ylabel('x1');
figure(2)%左轴颈中心轨迹图
plot(x(4000:end,1),x(4000:end,3)); title('轨迹图');xlabel('x1');ylabel('y1');
figure(3)%圆盘中心轨迹图
plot(x(4000:end,5),x(4000:end,7)); title('轨迹图');xlabel('x2');ylabel('y2');
回复 28楼 yangxiaokang 的帖子
呵呵，我有检查了一下程序，觉得计算程序没有问题，因为同样的程序对上田振子的庞加莱截面都是正确的，请仔细检查模型，以及方程的输入。当然还有一种情况，就是方程本身有某些特性，比如不稳定性，病态等，那样的话就不能用matlab自带的解法了，要自己编程 不知道有没有其他收藏的方法
我只好回帖了