求助：关于一个10自由度齿轮系统振动求解

tianxia01

function  dx=testfun(t,x)
%齿轮系统的基本参数
r1=0.4785;r2=0.253/2;r3=r2;r4=r2;r5=0.555/2;r6=r5;r7=r5;r8=0.145/2;r9=0.706/2;r10=0.149/2;
I1=189.7569;I2=1.2765;I3=I2;I4=I2;
I5=13.6031;I6=I5;I7=I5;I8=0.1320;I9=20.0138+2.5093;I10=0.09108;
m1=1003;m2=196;m3=m2;m4=m2;m5=263;m6=263;m7=263;m8=34.5;m9=367;m10=45.5;
z1=87;z2=23;z3=z2;z4=z2;z5=92;z6=z5;z7=z5;z8=24;z9=99;z10=21;
kn1=0.1363e9;kn2=0.01433e9;cn1=2*0.04*((kn1/((1/1.2765)+(1/13.6031)))^0.5);
cn2=2*0.04*((kn2/((1/0.132)+(1/I9)))^0.5);
%输入和输出的转矩
%Ts=9.88e6+1.236e6*(sin(2.361*t)^3);Tc=4954.47+495.45*sin(161.42*t);
%参数常数化尝试
Ts=9.88e6;Tc=4954.47;


%km矩阵表示的是齿轮啮合的刚度矩阵，cm矩阵表示的是齿轮啮合的阻尼矩阵；
%刚度
%km1=5.543e9/2+(6/(1*pi)*sin(4*1*pi/7)*cos(2*pi*1/0.03058)*t+(4/(1*pi)*cos(4/7*1*pi)-6/(1*pi))*sin((1*2*pi/0.03058)*t)...
%       +6/(2*pi)*sin(4*2*pi/7)*cos(2*pi*2/0.03058)*t+(4/(2*pi)*cos(4/7*2*pi)-6/(2*pi))*sin((2*2*pi/0.03058)*t)...
  %      +6/(3*pi)*sin(4*3*pi/7)*cos(2*pi*3/0.03058)*t+(4/(3*pi)*cos(4/7*3*pi)-6/(3*pi))*sin((3*2*pi/0.03058)*t)...
   %     +6/(4*pi)*sin(4*4*pi/7)*cos(2*pi*4/0.03058)*t+(4/(4*pi)*cos(4/7*4*pi)-6/(4*pi))*sin((4*2*pi/0.03058)*t)...
    %    +6/(5*pi)*sin(4*5*pi/7)*cos(2*pi*5/0.03058)*t+(4/(5*pi)*cos(4/7*5*pi)-6/(5*pi))*sin((5*2*pi/0.03058)*t))*(10^9);
    km1=5.543e9/2;
   
% km2=3.7e9/2+(0.6/(1*pi)*sin(2*pi*1/0.007646)*t-0.6/(1*pi)*cos((1*2*pi/0.007646)*t)...
      %  +0.6/(2*pi)*sin(2*pi*2/0.007646)*t-0.6/(2*pi)*cos((2*2*pi/0.007646)*t)...
       % +0.6/(3*pi)*sin(2*pi*3/0.007646)*t-0.6/(3*pi)*cos((3*2*pi/0.007646)*t)...
       % +0.6/(4*pi)*sin(2*pi*4/0.007646)*t-0.6/(4*pi)*cos((4*2*pi/0.007646)*t)...
       % +0.6/(5*pi)*sin(2*pi*5/0.007646)*t-0.6/(5*pi)*cos((5*2*pi/0.007646)*t))*(10^9);
    km2=3.7e9/2
   
km3=km2;
%阻尼
cm1=2*0.1*r1*r2*(km1*I1*I2/(r1^2*I1+r2^2*I2))^0.5;
cm2=2*0.1*r7*r8*(km2*I7*I8/(r1^2*I7+r2^2*I8))^0.5;
cm3=2*0.1*r9*r10*(km3*I9*I10/(r1^2*I9+r2^2*I10))^0.5;

%误差
%em1=0.02405*sin(2*pi*t/0.03058);
%em2=0.02055*sin(2*pi*t/0.001854);
%em3=0.01816*sin(2*pi*t/0.001854);
%em11=diff(em1);
%em22=diff(em2);
%em33=diff(em3);

%尝试常数化
em1=0.02405;
em2=em1;em3=em1;em11=em1;em22=em1;em33=em1;


%em矩阵表示的误差激励的矩阵；
%kn1表示的是第一根齿轮轴的刚度，kn2表示的第二根轴的刚度；
%cn1表示的是第一根齿轮轴的阻尼，cn2表示的第二根轴的阻尼；
%响应函数
%dx=zeros(20,1);
dx(1)=x(11);dx(2)=x(12);
dx(3)=x(13);dx(4)=x(14);
dx(5)=x(15);dx(6)=x(16);
dx(7)=x(17);dx(8)=x(18);
dx(9)=x(19);dx(10)=x(20);
dx(11)=(Ts-3*cm1*(r1^2)*x(11)+cm1*r1*r2*x(12)+cm1*r1*r3*x(13)+cm1*r1*r4*x(14)-3*cm1*r1*em11...
    -3*km1*(r1^2)*x(1)+km1*r1*r2*x(2)+km1*r1*r2*x(3)+km1*r1*r4*x(4)+3*km1*r1*em1)/I1;
dx(12)=(-cm1*r1*r2*x(11)-(cn1-cm1*(r2^2))*x(12)+cn1*x(15)-cm1*em11+km1*r1*r2*x(1)-(kn1-km1*r2*r2)...
   *x(2)+kn1*x(5)-km1*em1)/I2;
dx(13)=(-cm1*r1*r3*x(11)-(cn1-cm1*(r3^2))*x(13)+cn1*x(16)-cm1*em11+km1*r1*r3*x(1)-(kn1-km1*r3*r3)...
   *x(3)+kn1*x(6)-km1*r3*em1)/I3;
dx(14)=(-cm1*r1*r4*x(11)-(cn1-cm1*(r4^2))*x(14)+cn1*x(17)-cm1*em11+km1*r1*r4*x(1)-(kn1-km1*r4*r4)...
   *x(4)+kn1*x(7)-km1*r4*em1)/I4;
dx(15)=(cn1*x(12)-(cn1+cm2*r5*r5)*x(15)+cm2*r5*r8*x(18)+cm2*r5*em22+kn1*x(2)-(kn1+km2*r5*r5)*x(5)+...
   km2*r5*r8*x(8)+km2*r5*em2)/I5;
dx(16)=(cn1*x(13)-(cn1+cm2*r6*r6)*x(16)+cm2*r5*r8*x(18)+cm2*r6*em22+kn1*x(3)-(kn1+km2*r6*r6)*x(6)+...
   km2*r6*r8*x(8)+km2*r6*em2)/I6;
dx(17)=(cn1*x(14)-(cn1+cm2*r7*r7)*x(17)+cm2*r5*r8*x(18)+cm2*r7*em22+kn1*x(4)-(kn1+km2*r7*r7)*x(7)+...
   km2*r7*r8*x(8)+km2*r7*em2)/I7;
dx(18)=(-cm2*r5*r8*x(15)-cm2*r6*r8*x(16)-cm2*r7*r8*x(17)-(cn2-3*cm2*r8^2)*x(18)+cn2*x(19)-3*r8*cm2*em2-...
   km2*r5*r8*x(5)-km2*r6*r8*x(6)-km2*r7*r8*x(7)-(kn2-3*km2*r8*r8)*x(8)+kn2*x(9)-3*km2*r8*em2)/I8;
dx(19)=(cn2*x(18)-(cm3*r9*r9+cn2)*x(19)+cm3*r9*r10*x(20)+cm3*r9*em33+kn2*x(8)-(km3*r9*r9-kn2)*x(9)+...
   km3*r9*r10*x(10)+km3*r9*em3)/I9;
dx(20)=(-cm3*r9*r10*x(19)+cm3*r10*r10*x(20)-cm3*r10*em33-km3*r9*r10*x(9)+km3*r9*r9*x(10)-km3*r9*em3-Tc)/I10;
return

求解：
[t,x]=ode45('testfun',[0:pi/64:2*pi],[0 0.0007 0.0007 0.0007 0.0005 0.0005 0.0005 0.001 0.0026 0.013;0 0 0 0 0 0 0 0 0 0]);

出现的问题：

??? Error using ==> funfun\private\odearguments at 113
TESTFUN must return a column vector.

Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

[ 本帖最后由 tianxia01 于 2008-5-12 09:45 编辑 ]

tianxia01

参数是时变的，但现在即使把方程中的参数变成常数好像还是解不了，希望各位前辈能帮忙看一下问题出在哪

tianxia01

方程的基本格式就如下，一共有10个方程

咕噜噜

方程可以求解，子程序改一下，但是可能是你的参数或者初始条件有问题，算出来很多是无穷大
%dx=zeros(20,1);中%去掉
return去掉

无水1324

咕噜是对的
但是你这个方程最后求解的时候发散了