matlab ode45求解齿轮动力学,使用ODE45求解齿轮系统动力学方程后结果发散

function [dx,ff1,ff2]=myfun(t,x)

t

beita=26;

mn=0.004;

z1=46;

z2=43;

z3=122;

T_in=200;

T_out=80;

roug1=7.8E3;

roug2=7.8E3;

roug3=7.8E3;

alphan=20;

alphat=atand(tand(alphan)/cosd(beita));

d1=z1*mn/cosd(beita)/1000;

db1=d1*cosd(alphat)/1000;

d2=z2*mn/cosd(beita)/1000;

db2=d2*cosd(alphat)/1000;

d3=z3*mn/cosd(beita)/1000;

db3=d3*cosd(alphat)/1000;

bp1=116/1000;

bp2=116/1000;

bp3=116/1000;

bp=116/1000;

I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;

I2=((roug2*pi*(d2/2)^2*bp2)*(d2/2)^2)/2;

I3=((roug3*pi*(d3/2)^2*bp3)*(d1/2)^2)/2;

m1=roug1*pi*(d1/2)^2*bp1;

m2=roug2*pi*(d2/2)^2*bp2;

m3=roug3*pi*((d3)-(d1+d2))^2*bp3;

r1=d1*cosd(alphat)/2;

r2=d2*cosd(alphat)/2;

r3=d3*cosd(alphat)/2;

fai_sp1x=90;

fai_sp1y=0;

fai_p1rx=-130;

fai_p1ry=-220;

kesaiz=0.05;

kesain=0.07;

kp1x=1e8;

kp1y=kp1x;

cp1x=2*kesaiz*((kp1x*m2)^0.5);

cp1y=2*kesaiz*((kp1y*m2)^0.5);

ksx=1e8;

ksy=1e8;

csx=2*kesaiz*((ksx*m1)^0.5);

csy=2*kesaiz*((ksy*m1)^0.5);

krx=1e8;

kry=krx;

crx=2*kesaiz*((krx*m3)^0.5);

cry=crx;

Tmesh=2*pi/z1;

kp1r =1e6;

csp1=2*0.07*((ksp1/(1/m1+1/m2))^0.5);

cp1r=2*0.07*((ksp1/(1/m2+1/m3))^0.5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为参数定义，可忽略不看，谢谢！

esp1=1e-6;

ep1r=1e-6;

delta_sp1=((x(1)-x(7))*cosd(fai_sp1x)+(x(3)-x(9))*cosd(fai_sp1y)+r1*x(5)+r2*x(11))*cosd(beita)+esp1;

delta_sp11=((x(2)-x(8))*cosd(fai_sp1x)+(x(4)-x(10))*cosd(fai_sp1y)+r1*x(6)+r2*x(12))*cosd(beita)+esp1;

delta_p1r=((x(7)-x(13))*cosd(fai_p1rx)+(x(9)-x(15))*cosd(fai_p1ry)+r3*x(17)-r2*x(11))*sind(beita)-ep1r;

delta_p11r=((x(8)-x(14))*cosd(fai_p1rx)+(x(10)-x(16))*cosd(fai_p1ry)+r3*x(18)-r2*x(12))*sind(beita)-ep1r;

%%%%%%%%%%动力学方程

dx=zeros(18,1);

dx(1)=x(2);

dx(2)=(1500*cosd(fai_sp1x)-(csx*x(2))-ksx*x(1)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x))/m1;

dx(3)=x(4);

dx(4)=(1500*cosd(fai_sp1y)-(csx*x(4))-ksy*x(2)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y))/m1;

dx(5)=x(6);

dx(6)=(400-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*r1+T_in)/I1;

dx(7)=x(8);

dx(8)=((-300*cosd(fai_sp1x)+200*cosd(fai_p1rx)-(cp1x*x(8))-kp1x*x(7)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x)-(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m2;                %第一个行星齿轮x方向

dx(9)=x(10);

dx(10)=((-300*cosd(fai_sp1y)+200*cosd(fai_p1ry)-(cp1y*x(10))-kp1y*x(9)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m2;                 %第一个行星齿轮y方向

dx(11)=x(12);

dx(12)=(120-((csp1*delta_sp11+ksp1*delta_sp1)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*r2)/I2;

dx(13)=x(14);

dx(14)=((-200*cosd(fai_p1rx)-(crx*x(14))-krx*x(13)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m3;

dx(15)=x(16);

dx(16)=((-200*cosd(fai_p1ry)-(crx*x(16))-krx*x(15)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m3;

dx(17)=x(18);

dx(18)=(80-(cp1r*delta_p11r+kp1r*delta_p1r)*cosd(beita)*r3-T_out)/I3;