ode45求滚动轴承支撑下碰摩转子动力学方程，程序报错。

LebRonPeng

%%%%动力学方程程序
function dx=gd_pm1(t,x)
dx=zeros(12,1);
global w

%%求滚动轴承非线性作用力
global Nb Kb
Nb=8;                %滚珠个数
Kb=13.34e+09;        %轴承滚珠和滚道的接触刚度系数  单位 N/m^3/2    8.9176
rb=0.0401;           %轴承内圈半径  单位m
Rb=0.0639;           %轴承外圈半径  单位m
w0=2*pi*(w/60);      %转速和角速度换算公式
wb=w0*rb*(rb+Rb);    %轴承保持架的角速度
c=0.000008;          %轴承游隙

%左侧轴承1的非线性作用力；轴承1几何中心坐标（x2,y2）写为（x（5）,x（7））
[Fx_2 Fy_2] = Force_bearing(t,x(5),x(7),c,wb);
%右侧轴承2的非线性作用力；轴承2几何中心坐标（x3,y3）写为（x（9）,x（11））
[Fx_3 Fy_3] = Force_bearing(t,x(9),x(11),c,wb);  


%%系统相关系数
m1=32.1;            %圆盘处集中质量
m2=4;               %轴承1处、轴承2处集中质量
g=9.8;
c1=2100;            %圆盘处阻尼
c2=1050;            %轴承处阻尼
k0=0.85e7;          %轴系刚度系数
e=0.00005;          %质量偏心距离

%%无量纲化
% wn=(k0/m1)^0.5;
kesai1=c1/(m1*w);   %无量纲阻尼
kesai2=c2/(m2*w);
numda1=k0/(m1*w^2);
numda2=k0/(m2*w^2);
gg=g/(c*w^2);
U=e/c;
% omega=w/wn;
beta=0;

%%求碰摩力
ku=0.1;                 %摩擦系数
kc=3.5e7;               %碰摩刚度系数
r=sqrt(x(1)^2+x(3)^2);  %圆盘几何中心位移
deta=0.00002;           %碰摩间隙
H=0.5*(sign(abs(r-deta))+sign(r-deta));%判断是否发生碰摩，sigh正数为1，负数为-1
px=H*(-c*(r-deta)*kc*(x(1)-ku*x(3))/r);
py=H*(-c*(r-deta)*kc*(ku*x(1)+x(3))/r);

%%方程
dx(1)=x(2);
dx(2)=-kesai1*x(2)-numda1*(x(1)-x(5))-numda1*(x(1)-x(9))+U*cos(t-beta)+px/(w^2*m1*c);
dx(3)=x(4);
dx(4)=-kesai1*x(4)-numda1*(x(3)-x(7))-numda1*(x(3)-x(11))+U*sin(t-beta)+py/(w^2*m1*c)-gg;

dx(5)=x(6);
dx(6)=-kesai2*x(6)-numda2*(x(5)-x(1))+Fx_2/(w^2*m2*c);
dx(7)=x(8);
dx(8)=-kesai2*x(8)-numda2*(x(7)-x(3))+Fy_2/(w^2*m2*c)-gg;

dx(9)=x(10);
dx(10)=-kesai1*x(10)-numda1*(x(9)-x(1))+Fx_3/(w^2*m2*c);
dx(11)=x(12);
dx(12)=-kesai1*x(12)-numda1*(x(11)-x(3))+Fy_3/(w^2*m2*c)-gg;


%%%%求解程序
tic
clear all
clc

global w

w=400;


T=2*pi;
delt_T=T/100;
x0=[0.00001;0;0.00001;0;0.00001;0;0.00001;0;0.00001;0;0.00001;0];
tspan=[0:delt_T:100*T];      
options = odeset('RelTol', 1e-6,'AbsTol',1e-6);
[t,x]=ode45('gd_pm1',tspan,x0,options);
x0=x(end,:);
tspan=[100*T:delt_T:200*T];
[t x]=ode45('gd_pm1',tspan,x0,options);  

figure
%%time domain plot
subplot(4,1,1)
u=[];
u=x(1:end,1)-mean(x(1:end,1));
plot(t(1:end),u(1:end,1));%,set(gca,'xlim',[0,500])
xlabel('t/s');
ylabel('x');
title('(a)')

subplot(4,1,2)  
%%axis orbit
plot(x(1:end,1),x(1:end,3))
xlabel('x');
ylabel('y');
title('(b)')

%% FFT
fs=100/(2*pi);
N = length(x(1:end,1));
df = 0:fs/N:(N-1)*fs/N;
a= abs(fft(x(1:end,1)-mean(x(1:end,1))))*2/N;
DF=2*pi*df;

subplot(4,1,3)
plot(DF,a),set(gca,'xlim',[0,5])
xlabel('f/fr');
ylabel('\rmAmplitude');
title('(c)')

%%poincare map
m=[];
n=[];
m=x(1:100:end,1);
n=x(1:100:end,2);
subplot(4,1,4)
plot(m,n,'k.')
xlabel('x');
ylabel('x’');
title('(d)')


%%%%轴承力程序
function [Fx_bearing Fy_bearing]= Force_bearing(t,x,y,c,wb)  %该函数用于求滚动轴承的非线性支撑力

global  Nb Kb

theta_ball=zeros(1,Nb);
for i=1:Nb
    theta_ball(i) = 2*pi*(i-1)/Nb;    %各转动体初始角
end
theta_change=wb*t;    %转动体角度变动量
theta_ball=theta_ball+theta_change;  %各滚珠的各时变转角
sigma_ball= x*cos(theta_ball)+y*sin(theta_ball)-c;  %各滚珠的法向接触变形量
Heaviside_1 =  heaviside(sigma_ball);   % heaviside 函数用于筛选 sigma_ball>0的变形量
useful_sigma = sigma_ball.*Heaviside_1;
sigma_1=useful_sigma.^3;
sigma = sqrt(sigma_1);
fix_bearing = (sigma.*cos(theta_ball));  %轴承支撑力的x方向分量
fiy_bearing = (sigma.*sin(theta_ball));  %轴承支撑力的y方向分量

Fx_bearing=Kb * sum(fix_bearing);
Fy_bearing=Kb * sum(fiy_bearing);

end


将三段程序分别保存为3个m文件，运用求解程序，出现报错。求助坛内大神，小弟应该做如何修改。
先行谢过了！！！