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

大小的小

1.1参数定义及动力学方程降阶
function [dx,ff1,ff2]=myfun(t,x)
) m# B  C8 P. A/ V1 b4 ]# zt
beita=26;
mn=0.004;
z1=46;                       
: i) C) ~. |4 C. }! b8 a( S9 ^z2=43;                       
z3=122;                     
0 X! t- u/ V, h* N0 _: x1 TT_in=200;
( }+ t% F5 A! gT_out=80;
roug1=7.8E3;            
roug2=7.8E3;            
roug3=7.8E3;            
alphan=20;                                                   
, Q0 Q1 {/ f3 W- Halphat=atand(tand(alphan)/cosd(beita));      
d1=z1*mn/cosd(beita)/1000;                  
1 q# U( G+ M: o8 ]/ H% }db1=d1*cosd(alphat)/1000;                    
* B7 u* g1 H7 cd2=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;                                       
  e& A8 @$ d* E7 m/ M; `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;  
$ j0 T  `; i4 o. t, f9 K( N, |/ A6 q; gI3=((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;                  
  A- s' _5 R5 a7 `+ V* S# q4 Jr1=d1*cosd(alphat)/2;                                          
1 Z" r. D& z+ X3 Br2=d2*cosd(alphat)/2;                                          
: [; f' t! f9 ]r3=d3*cosd(alphat)/2;                                          
0 t2 T8 K; l4 m' }& q$ qfai_sp1x=90;
fai_sp1y=0;
* D; t- @+ e8 O3 hfai_p1rx=-130;
& f& P, m% K' ~  Zfai_p1ry=-220;
kesaiz=0.05;
kesain=0.07;
+ Q' G. G" I4 B. _3 Ikp1x=1e8;
& H; g) ^3 ]- H" }kp1y=kp1x;
! a0 ]& B& R& a5 p" T7 @. R: c7 fcp1x=2*kesaiz*((kp1x*m2)^0.5);
9 X* @' d9 J. Ycp1y=2*kesaiz*((kp1y*m2)^0.5);
  \3 A% ]- t" D- v8 u, tksx=1e8;
% m7 m+ O7 A5 ~8 Iksy=1e8;                                   
/ |+ e: J$ l( T0 T' b7 Scsx=2*kesaiz*((ksx*m1)^0.5);
csy=2*kesaiz*((ksy*m1)^0.5);
krx=1e8;
1 d' z% D5 d" C, _kry=krx;
8 ~9 o# i9 a! ]2 W& g) e. Wcrx=2*kesaiz*((krx*m3)^0.5);   
( u' N% D, H, g# M; ?cry=crx;
Tmesh=2*pi/z1;
kp1r =1e6;
/ L6 J% \4 `. o8 Mcsp1=2*0.07*((ksp1/(1/m1+1/m2))^0.5);
8 j3 t; D( Q* ~$ g' h4 fcp1r=2*0.07*((ksp1/(1/m2+1/m3))^0.5);                 
" n3 [# b: t) M%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为参数定义，可忽略不看，谢谢！
esp1=1e-6;
ep1r=1e-6;
) U8 c# p0 K  @- ~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;
0 a* f$ I2 m" X8 `* @# q9 Ydelta_p1r=((x(7)-x(13))*cosd(fai_p1rx)+(x(9)-x(15))*cosd(fai_p1ry)+r3*x(17)-r2*x(11))*sind(beita)-ep1r;
7 f! I/ O$ P& r! w, S" idelta_p11r=((x(8)-x(14))*cosd(fai_p1rx)+(x(10)-x(16))*cosd(fai_p1ry)+r3*x(18)-r2*x(12))*sind(beita)-ep1r;
%%%%%%%%%%动力学方程
9 l& \* i7 P3 W2 j4 xdx=zeros(18,1);
: |8 z4 G/ T: {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;
1 e7 R0 @6 l. _2 C) M. P% H( cdx(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;
5 |% q# ?; ]6 ~. F- Fdx(7)=x(8);
$ b; q0 L7 I3 @' k  ~, C# R" T! _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方向
: `8 Z' G0 k* @+ \/ adx(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方向
3 u5 [/ r: G  n( {1 {, I6 rdx(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);
; ?0 |; Y) J3 i9 V9 Y4 I" rdx(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;

+ Y0 o& c- ^+ Q' p- D7 H2 u3 ~/ E
! K2 v& d- i5 ^! L7 t) G
5 b$ n0 C( j& a" k6 ~. [1.2 ode程序
3 q8 |; y$ r  }# |4 F2 Vclc;
clear all
% b% P3 j8 t5 ~: a  `1 Z' k, Dx0=zeros(18 ,1)
[t,x] = ode45('myfun',[0:0.0001:10],x0);
figure
" s6 u2 w9 v$ A  o4 ^& w: gplot(t,x(:,1))
+ x5 J, o$ x9 K% e* A
( J$ I! g' X* h4 Y& c  {# L
3 ~9 l6 Y/ Q6 V3.绘图结果如下，为什么画出来是一条直线，而且图中结果没有计算到规定的时间
+ T; b5 D* Z* }3 u: Z

/ b" F" R' F( f# z# Q' T
# F& f6 J/ Z+ z3 ]
+ a! k' u# e4 F' L) @- V9 C: I

Uifhjvv

帮你顶一下

IBB-EUT

你这个自定义函数不完整吧，找不到对 ksp1 这个变量的有效赋值，csp1等等需要调用ksp1的赋值也没法做

shuddkk

function [dx,ff1,ff2]=myfun(t,x)
t
; A; W5 w1 w8 V. l, vbeita=26;
+ E* V6 q& b: M9 M) kmn=0.004;
) p: A+ S' M" xz1=46;
. s2 ~% g4 ^0 o& j2 i, kz2=43;
& r+ Y! V1 U# }8 `0 c9 Q& ]. `z3=122;
T_in=200;
T_out=80;
roug1=7.8E3;
roug2=7.8E3;
0 L7 K8 I0 e; I" Z- }3 @( \$ qroug3=7.8E3;
- z2 s6 u! T1 Y( Galphan=20;
alphat=atand(tand(alphan)/cosd(beita));
0 Q* m/ d* k4 Z" Id1=z1*mn/cosd(beita)/1000;
db1=d1*cosd(alphat)/1000;
& g, t  r( m) q% S# \d2=z2*mn/cosd(beita)/1000;
db2=d2*cosd(alphat)/1000;
d3=z3*mn/cosd(beita)/1000;
; x$ y* M$ z- ]1 j6 ?db3=d3*cosd(alphat)/1000;
* o) v; C9 ^& B8 s; kbp1=116/1000;
bp2=116/1000;
1 ]( @/ d% C+ bbp3=116/1000;
/ h% X+ o3 `5 g7 ^$ Zbp=116/1000;
I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;
6 Y% S3 k, b1 e$ tI2=((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;
% G1 V% E# Q( C7 O/ G# Cm2=roug2*pi*(d2/2)^2*bp2;
# d. h7 C/ v1 V  km3=roug3*pi*((d3)-(d1+d2))^2*bp3;
& T9 U! K' {2 W  ~9 [+ b' er1=d1*cosd(alphat)/2;
: n3 K: D$ b* F! j4 b9 Mr2=d2*cosd(alphat)/2;
r3=d3*cosd(alphat)/2;
. |( Y! L5 P4 @7 k4 x, A$ afai_sp1x=90;
fai_sp1y=0;
$ x6 ~7 [+ o) w. \7 Zfai_p1rx=-130;
fai_p1ry=-220;
kesaiz=0.05;
3 h& `% O8 a+ ]kesain=0.07;
1 F! j- `: `1 @2 Bksp1 =1e6;
, f! S6 y8 r/ R! O* b4 U$ S/ zkp1r =1e6;
kp1x=1e8;
kp1y=kp1x;
cp1x=2*kesaiz*((kp1x*m2)^0.5);
cp1y=2*kesaiz*((kp1y*m2)^0.5);
  X7 L& {) T+ R% a& Hksx=1e8;
ksy=1e8;
csx=2*kesaiz*((ksx*m1)^0.5);
  y" \( ?3 ^- |8 n$ a3 @csy=2*kesaiz*((ksy*m1)^0.5);
  L/ E' x3 f8 Rkrx=1e8;
3 C0 K" Y+ H2 {( Q. \kry=krx;
crx=2*kesaiz*((krx*m3)^0.5);
  q8 X) b, W% `2 v+ Kcry=crx;
& t  r- L9 H% Z+ c. ~# bTmesh=2*pi/z1;
kp1r =1e6;
% H( W$ y0 d2 X. T+ H! scsp1=2*0.07*((ksp1/(1/m1+1/m2))^0.5);
cp1r=2*0.07*((ksp1/(1/m2+1/m3))^0.5);
: G  Y7 x' F& Z%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为参数定义，可忽略不看，谢谢！
+ a' h* R6 B+ Y7 ~/ Sesp1=1e-6;
ep1r=1e-6;
3 r. [* V1 p& Fdelta_sp1=((x(1)-x(7))*cosd(fai_sp1x)+(x(3)-x(9))*cosd(fai_sp1y)+r1*x(5)+r2*x(11))*cosd(beita)+esp1;
  r7 o5 c9 V4 Y# Y9 _1 T2 r0 Ndelta_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;
2 [. m) e6 |5 Z) Y. S%%%%%%%%%%动力学方程
dx=zeros(18,1);
7 Q- w- z' S. Ldx(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;
3 L- K4 o6 U, C1 ?6 x0 Z& _* X* Jdx(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);
4 y6 I! g$ E- \- E5 adx(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;
7 [; F; R, V3 U$ k3 _& m- @/ Ydx(17)=x(18);
dx(18)=(80-(cp1r*delta_p11r+kp1r*delta_p1r)*cosd(beita)*r3-T_out)/I3;

lkk

这么专业的问题还是求助淘宝吧