请教齿轮动力学问题

tianlei2005

我在学齿轮动力学的仿真，一个3自由度的扭转振动方程。挺简单的一个方程，可是每次运算都花费好几个小时，
采用的是4阶定步长Runge-Kutta积分方法，所以请问是不是ode45算法中的参数设置的问题还是其它原因，特请教各位高人，
下面是两个文件的代码
1.TorsionVibration.m文件
function f=TorsionVibration(t,y);
%齿轮系统参数
Jm=0.1;Jp=1.772;Jg=1.873;Jf=0.15;
Rp=0.1255;Rg=0.256;R1=0.0335;R2=0.0575;
cp=5.55;cm=504.75;cg=7.18;
kp=3.0869e6;km=1.12e8;kg=5.5998e6;
Me=Jp*Jg/(Rp^2*Jg+Rg^2*Jp);
M=diag([Jm;Jp;Jg;Jf]);
Tm=5000;Tf=5000*47*0.9/23;
a11=-cp*(1/Jm+1/Jp);a12=-kp*(1/Jm+1/Jp);a13=km*R1*Rp/Jp;a14=cm*R1*Rp/Jp;
a21=kp*Rp/(Jp*R1);a22=cp*Rp/(Jp*R1);a23=-km/Me;a24=-cm/Me;a25=kp*Rp/(Jg*R2);a26=cp*Rp/(Jg*R2);
a35=-cg*(1/Jg+1/Jf);a36=-kg*(1/Jg+1/Jf);a33=km*R2*Rg/Jg;a34=cm*R2*Rg/Jg;
U=[0 1 0 0 0 0;
    a11 a12 a13 a14 0 0;
    0 0 0 1 0 0;
    a21 a22 a23 a24 a25 a26;
    0 0 0 0 0 1;
    0 0 a33 a34 a35 a36];
f=U*y+[0;Tm*R1/Jm;0;0;0;Tf*R2/Jf];

2.TorsionVibrationResult.m文件
[t,y]=ode45('TorsionVibration',[0:0.001:10],[1.2e-3 0 2.56e-4 0 1.45e-3 0]');
u1=y(:,1);
u2=y(:,3);
u3=y(:,5);
u4=y(:,6);
figure(1);
plot(t,u1);
figure(2);
plot(t,u2);
figure(3);
plot(t,u3);
figure(4);
plot(t,u4);

zhlong

能否把相关的文字说明也发上来？
动力学方程和图形什么的。