用ODE解龙格库塔方程遇到的问题
我要用ODE解个16自由度的常微分方程，然后画图，画出来应该是周期的图，可不知道那里错了，根本出不来周期，那为能帮下忙啊方程的M文件如下
function dK=chilun2222(T,K)
c=10.078019;
q=3.14/180;
c=c*q;
a =-0.6419
oa=104.53;
oae=49.03;
oae=oae*q;
bao=130.97;
bao=bao*q;
aoi=-35.7794;
aoi=aoi*q;
aoe=69.6775;
aoe=aoe*q;
oea=61.28;
oea=oea*q;
qea=67.97;
qea=qea*q;
pe=14.057;
ae=111.767;
qae=56.24;
qae=qae*q;
oi=78.93;
oc=oa*sin(oae)/sin(bao-a);
ac=oc*sin(a)/sin(oae);
b=sqrt((oc/oi)^2-1)-atan(sqrt((oc/oi)^2-1))-aoi;
of=oc*sin(bao+a);
ec=oc*sin(aoe-a)/sin(oea);
d=atan(ac*sin(aoe)*sin(qae)/(ec*sin(qea)+ac*sin(qae)*cos(aoe)));
qc=ac*sin(qae)/sin(d);
qh=qc*cos(d+qae);
Il=5.0291288e+05;
Ip=1.2070271e+04;
Ig=5.0291288e+04;
Im=1.2070271e+05;
Kb=7.5;
Kc=7.5;
Ks=7.5;
M=15;
Mp=3.2194583e+00;
Mg=6.5646591e+00;
Qb=3.2;
Qc=3.2;
Km=1739358.884722;
Rp=84*cos(20*q);
Rg=120*cos(20*q);
Kmb=Km/Rp^2;
e=0.1;
Qmb=2*e*sqrt((Km*Rp^2*Rg^2*Ip*Ig)/(Rp^2*Ip+Rg^2*Ig));
Fc=Kmb*(Qmb*0.2-Rp*142+Rg*139);
U=0.3;
Ff=U*Fc;
w1=141;
w2=141;
w3=96;
w4=96;
v1=0.5;
v2=0.5;
v3=0.5;
v4=0.5;
v5=0.5;
v6=0.5;
v7=0.5;
v8=0.5;
v9=0.5;
v10=0.5;
v11=0.5;
v12=0.5;
Tin=1193.5;
Tout=1193.5*0.97;
dK=[K(2);
(Tin-Qc*(K(2)-K(4))-Kc*(K(1)-K(3)))/Im;
K(4);
(Kc*(K(3)-K(1))+Rp*Kmb*(Rp*K(3)-Rg*K(5)-K(11)+K(17))+Qc*(K(4)-K(2))+Rp*Qmb*(Rp*K(4)-Rg*K(6)-K(12)+K(18))+Ff*of)/(-Ip);
K(6);
(Ff*qh-Kc*(K(5)-K(7))-Rg*Kmb*(Rg*K(5)-Rp*K(3)-K(17)+K(11))-Qc*(K(6)-K(8))-Rg*Qmb*(Rg*K(6)-Rp*K(4)-K(18)+K(12)))/Ig;
K(8);
(Tout+Qc*(K(8)-K(6))+Kc*(K(7)-K(5)))/(-Il);
K(10);
(Qb*K(10)+Kb*K(9)+Ks*(K(9)-K(11)))/(-M);
K(12);
(Ks*(K(11)-K(9))+Ks*(K(11)-K(13))+Kmb*(K(11)-K(17)-Rp*K(3)+Rg*K(5))+Qmb*(K(12)-K(18)-Rp*K(4)+Rg*K(6)))/(-Mp);
K(14);
(Qb*K(14)+Kb*K(13)+Ks*(K(13)-K(11)))/(-M);
K(16);
(Qb*K(16)+Kb*K(15)+Ks*(K(15)-K(17)))/(-M);
K(18);
(Ks*(K(17)-K(15))+Ks*(K(17)-K(19))+Kmb*(K(17)-K(11)-Rg*K(5)+Rp*K(3))+Qmb*(K(18)-K(12)-Rg*K(6)+Rp*K(4)))/(-Mg);
K(20);
(Qb*K(20)+Kb*K(19)+Ks*(K(19)-K(17)))/(-M);
K(22);
(Qb*K(22)+Kb*K(21)+Ks*(K(21)-K(23)))/(-M);
K(24);
(Ff-Ks*(K(23)-K(21))-Ks*(K(23)-K(25)))/Mp;
K(26);
(Qb*K(26)+Qb*K(25)+Ks*(K(25)-K(23)))/(-M);
K(28);
(Qb*K(28)+Kb*K(27)+Ks*(K(27)-K(29)))/(-M);
K(30);
(Ks*(K(29)-K(27))+Ks*(K(29)-K(31))+Ff)/(-Mg);
K(32);
(Qb*K(32)+Kb*K(31)+Ks*(K(31)-K(29)))/(-M)]
命令文件如下
K0=[1141;0.0099;141;-5.4278e+006;96;1.7502e+006;96;-0.0023;5;
-1.0667;5;2.0465e+008;5;-1.0667;5;-1.0667;5;-1.0036e+008;5;-1.0667;5;-1.0667;5;
1.1373e+007;5;-1.0667;5;-1.0667;5;-5.5774e+006;5;-1.0667]
[T,Y]=ode15s('chilun22223',[0, 2000],K0);
plot(Y(:,1),Y(:,2))
