油膜力的程序是不对的，总有虚数部分泄露、
请参考
function oilforce=oilx(x,y,dx,dy,wi)
%R轴承半径，L轴承长度，wi激励频率，miu油膜黏度 dert为油膜平均厚度
R=0.025; L=0.012; miu=0.018;dfai=1; dert=0.00011;
e=sqrt(x*x+y*y);
%delta为Sommerfeld修正系数
delta=miu*wi*R*L*(R/dert)*(R/dert)*(L/2.0/R)*(L/2.0/R);
%
ppp1=(y+2.0*dx)/(x-2.0*dy);
sign1=sign(ppp1);
ppp2=y+2.0*dx;
sign2=sign(ppp2);
alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((y*cos(alpha)-x*sin(alpha))/sqrt(abs(1.0-abs(x*x)-abs(y*y))));
fg=2.0*(pi/2.0+alphaa)/sqrt(abs(1.0-abs(x*x)-abs(y*y)));
fv=(2.0+(y*cos(alpha)-x*sin(alpha))*fg)/(1.0-abs(x*x)-abs(y*y));
fs=(x*cos(alpha)+y*sin(alpha))/(1.0-abs((x*cos(alpha)+y*sin(alpha))*(x*cos(alpha)+y*sin(alpha))));
f1=sqrt(abs(abs((x-2.0*dy)*(x-2.0*dy))+abs((y+2.0*dx)*(y+2.0*dx))))/(1.0-abs(x*x)-abs(y*y));
fx=-1.0*f1*(3.0*x*fv-sin(alpha)*fg-2.0*cos(alpha)*fs);
fy=-1.0*f1*(3.0*y*fv+cos(alpha)*fg-2.0*sin(alpha)*fs);
得出的结果与原文的差距有点离谱，是不是应该换种积分方法？
我也正做转子动力学呢 才列出来动力学方程 还没无量钢化呢啊
{:{40}:}{:{40}:}
学习了，最近也在研究这方面
这类转子油膜碰摩的文献里，大多数根据别人论文里的参数却做不出同样的结果来。
这种现象我遇见的太多了。真不知道怎么去解释好。同样的参数，按照别人同样的处理计算方法，同样的取点，就是做不出他们的那种结果来。可能是我自己有些问题，但有一些人的论文里就是有很明显错误的。
愿与广大同学共勉！
可能是我自己有些问题，但有一些人的论文里就是有很明显错误的。
hsfy919 发表于 2010-3-16 23:09 static/image/common/back.gif
帮你调试了一下，程序通了，但不知道是不是你要的
%%%%%%%%%%%%%%%%%%%函数文件
function dx=fangcheng_f ...
21楼中定义的global w f f=400:5:600和前面函数中f=0.1会不会有冲突？
另外不知道你的分岔图怎么样了？现在对了吗？
“我自己认为是先做参数随频率的分岔图，然后在分岔出现的地方，取不同频率，再做轨迹图和映射图分析，不知道对否”
这个想法是对的，但是有时候不好选择参数的范围的时候，可以先尝试着计算一些点的轨迹图和映射图