主程序调用子程序提示以下错误

dj15175771160

Jeffcott刚性转子，采用Capone油膜力，分析它的时间历程
子程序如下：
function dy=fuhuan140922(t,x,flag,wj,p)

uo=0.018;         %外油膜粘度,pa.s
ui=0.018;         %内油膜粘度,pa.s
co=0.000114;      %外油膜间隙,m
ci=0.000093;      %内油膜间隙,m
Rro=0.018830;     %浮环外径,m
Rri=0.014093;     %浮环内径,m
Rj=0.014;         %轴颈半径,m
L=0.018;          %轴承宽度,m
mj=0.018257/2;    %轴子质量之半,kg
mr=0.000079/2;    %浮环质量之半,kg
g=9.8;            %重力加速度,N/kg(m/s^2)
simga=0.018;      %润滑油的动力粘度
wr=wj/(1+(uo/ui)*(ci/co)*(Rro/Rri)^3);             %浮环的转速
P=mj*g;
Po=mr*g;
detai=simga*wj*L*Rj*(Rj/ci)^2*(L/(2*Rj))^2;        %转子Sommerfeld修正数
detao=simga*wr*L*Rro*(Rro/co)^2*(L/(2*Rro))^2;     %浮环Sommerfeld修正数  
Mj=mj*ci*wj.*wj./detai;
Mr=mr*co*wr.*wr./detao;
Mk=mr*co*wr.*wr./detai;
Gi=g./(ci*wj.^2);
Go=g./(co*wr.^2);

Xj=x(1);dXj=x(2);Yj=x(3);dYj=x(4);
Xr=x(5);dXr=x(6);Yr=x(7);dYr=x(8);

%以下为求解“外层”油膜力在x,y方向上的无量纲分力——————————————————
ppp1=(Yr+2.00*dXr)/(Xr-2.00*dYr);
sign1=sign(ppp1);  
ppp2=Yr+2.00*dXr;  
sign2=sign(ppp2);

alpha=atan(ppp1)-pi/2.0*(sign1+sign2);
alphaa=atan((Yr*cos(alpha)-Xr*sin(alpha))/sqrt(abs(1.00-abs(Xr*Xr)-abs(Yr*Yr))));
fgo=2.00*(pi/2.00+alphaa)/sqrt(abs(1.00-abs(Xr*Xr)-abs(Yr*Yr)));
fvo=(2.00+(Yr*cos(alpha)-Xr*sin(alpha))*fgo)/(1.00-abs(Xr*Xr)-abs(Yr*Yr));
fso=(Xr*cos(alpha)+Yr*sin(alpha))/(1.00-abs((Xr*cos(alpha)+Yr*sin(alpha))*(Xr*cos(alpha)+Yr*sin(alpha))));

fao=sqrt(abs(abs((Xr-2.00*dYr)*(Xr-2.00*dYr))+abs((Yr+2.00*dXr)*(Yr+2.00*dXr))))/(1.00-abs(Xr*Xr)-abs(Yr*Yr));
fxo=1.00*fao*(3.00*Xr*fvo-sin(alpha)*fgo-2.00*cos(alpha)*fso);
fyo=1.00*fao*(3.00*Yr*fvo+cos(alpha)*fgo-2.00*sin(alpha)*fso);
%以上为求解“外层”油膜力在x,y方向上的无量纲分力——————————————————

%以下为求解“内层”油膜力在x,y方向上的无量纲分力——————————————————
ppp3=(1.34*Yj+2.00*dXj)/(1.34*Xj-2.00*dYj);
sign3=sign(ppp3);  
ppp4=1.34*Yj+2.00*dXj;  
sign4=sign(ppp4);

alpha1=atan(ppp3)-pi/2.00*(sign3+sign4);
alphaa1=atan((Yj*cos(alpha1)-Xj*sin(alpha1))/sqrt(abs(1.00-abs(Xj*Xj)-abs(Yj*Yj))));
fgi=2.00*(pi/2.00+alphaa1)/sqrt(abs(1.00-abs(Xj*Xj)-abs(Yj*Yj)));
fvi=(2.00+(Yj*cos(alpha1)-Xj*sin(alpha1))*fgi)/(1.00-abs(Xj*Xj)-abs(Yj*Yj));
fsi=(Xj*cos(alpha1)+Yj*sin(alpha1))/(1.00-abs((Xj*cos(alpha1)+Yj*sin(alpha1))*(Xj*cos(alpha1)+Yj*sin(alpha1))));

fai=sqrt(abs(abs((Xj-2.00*dYj)*(Xj-2.00*dYj))+abs((Yj+2.00*dXj)*(Yj+2.00*dXj))))/(1.00-abs(Xj*Xj)-abs(Yj*Yj));
fxi=-1.00*fai*(3.00*Xj*fvi-sin(alpha1)*fgi-2.00*cos(alpha1)*fsi);
fyi=-1.00*fai*(3.00*Yj*fvi+cos(alpha1)*fgi-2.00*sin(alpha1)*fsi);
%以上为求解“内层”油膜力在x,y方向上的无量纲分力——————————————————

fxi=fxi./detai;  fyi=fyi./detai;
fxo=fxo./detao;  fyo=fyo./detao;

% 以下为浮环轴承系统动力学模型进行无量纲化处理———————————————————

dy=zeros(8,1);
dy=[x(2);
   fxi./Mj+p*sin(t);
   x(4);
   fyi./Mj+p*cos(t)+Gi;
   x(6);
   fxo./Mr-fxi./Mk;
   x(8);
   fyo./Mr-fyi./Mk+Go];

主程序如下：
function  fuhuanjeffcott_solution
p=0;            %偏心率
n=12000;         %转子的转速
wj=2*pi*n/60;   %转子的角频率
%uo=0.018;       %外油膜粘度,pa.s
%ui=0.018;       %内油膜粘度,pa.s
%co=0.000114;    %外油膜间隙,m
%ci=0.000093;    %内油膜间隙,m
%Rro=0.018830;   %浮环外径,m
%Rri=0.014093;   %浮环内径,m
%Rj=0.014;       %轴颈半径,m
%L=0.018;        %轴承宽度,m
%wr=wj/(1+(uo/ui)*(ci/co)*(Rro/Rri)^3);   %浮环的转速

x0=[-0.2 0 -0.2 0 -0.2 0 -0.2 0];        %赋初值,另一组初值[0 0.1 0 0.1];
options=odeset('RelTol',1.e-7,'AbsTol',[1.e-8 1.e-8 1.e-8 1.e-8 1.e-8 1.e-8 1.e-8 1.e-8] );  %控制精度

T=2*pi;            %一个周期
T1=2*pi/wj;

%
figure           % 1时间历程
%[t,X]=ode45('fuhuan140922',[0:T1/100:14000*T1],x0,options，wj,p);  
[t,X]=ode45('fuhuan140922',[0:0.01:2000],x0,options,wj,p);  
subplot(121)   
plot(t(100000:end),X(100000:end,1))
xlabel('t')
ylabel('Yj')     %  grid on
title('轴颈时间历程x-t')

提示错误如下：
Warning: Failure at t=6.250756e-004.  Unable to meet integration tolerances without
reducing the step size below the smallest value allowed (1.734723e-018) at time t.
> In ode45 at 371
  In fuhuanjeffcott_solution at 24

meiyongyuandeze

可能是你option中的收敛误差设置的太小了，积分残差不满足设置的值！