F25=0.0;

FM25=-RM25.*9.8.*U.*VP1./OMGA;

%弹筒26的运动规律，认为弹筒26在90°周期内作平动：

RM26=68.8398;

RJ26=0.294227;

U=0.15;

RJ26=RM26.*(VP1./OMGA).^2;

F26=0.0;

FM26=-RM26.*9.8.*U.*VP1./OMGA;

%弹筒27的运动规律，认为弹筒27在90°周期内作平动：

RM27=68.8398;

RJ27=0.294227;

U=0.15;

RJ27=RM27.*(VP1./OMGA).^2;

F27=0.0;

FM27=-RM27.*9.8.*U.*VP1./OMGA;

%弹筒28的运动规律，认为弹筒28在90°周期内作平动：

RM28=68.8398;

RJ28=0.294227;

U=0.15;

RJ28=RM28.*(VP1./OMGA).^2;

F28=0.0;

FM28=-RM28.*9.8.*U.*VP1./OMGA;

%弹筒29的运动规律，认为弹筒29在90°周期内作平动：

RM29=68.8398;

RJ29=0.294227;

U=0.15;

RJ29=RM29.*(VP1./OMGA).^2;

F29=0.0;

FM29=-RM29.*9.8.*U.*VP1./OMGA;

%弹筒30的运动规律，认为弹筒30在90°周期内作平动：

RM30=68.8398;

RJ30=0.294227;

U=0.15;

RJ30=RM30.*(VP1./OMGA).^2;

F30=0.0;

FM30=-RM30.*9.8.*U.*VP1./OMGA;

%弹筒31的运动规律，认为弹筒31在90°周期内作平动：

RM31=68.8398;

RJ31=0.294227;

U=0.15;

RJ31=RM31.*(VP1./OMGA).^2;

F31=0.0;

FM31=-RM31.*9.8.*U.*VP1./OMGA;

EJ=RJ1+RJ2+RJ3+RJ4+RJ5+RJ6+RJ7+RJ8+RJ9+RJ10+RJ11+RJ12+RJ13+RJ14+RJ15+RJ16+RJ17+RJ18+RJ19+RJ20+RJ21+RJ22+RJ23+RJ24+RJ25+RJ26+RJ27+RJ28+RJ29+RJ30+RJ31;

EF=F1+F2+F3+F4+F5+F6+F7+F8+F9+F10+F11+F12+F13+F14+F15+F16+F17+F18+F19+F20+F21+F22+F23+F24+F25+F26+F27+F28+F29+F30+F31;

EFMDT=FM1+FM2+FM3+FM4+FM5+FM6+FM7+FM8+FM9+FM10+FM11+FM12+FM13+FM14+FM15+FM16+FM17+FM18+FM19+FM20+FM21+FM22+FM23+FM24+FM25+FM26+FM27+FM28+FM29+FM30+FM31;

end

DJDF(1)=(EJ(2)-EJ(1))./(XTA(2)-XTA(1));

DJDF(N+1)=(EJ(N+1)-EJ(N))./(XTA(N+1)-XTA(N));

DJDF(N/2+1)=0;

DJDF(N/2+2)=0;

for i=2:(N/2);

DJDF(i)=(EJ(i+1)-EJ(i))./(XTA(i+1)-XTA(i))./2.0+(EJ(i)-EJ(i-1))./(XTA(i)-XTA(i-1))./2.0;

end

for i=(N/2+1):N;

    DJDF(i)=(EJ(i+1)-EJ(i))./(XTA(i+1)-XTA(i))./2.0+(EJ(i)-EJ(i-1))./(XTA(i)-XTA(i-1))./2.0;

end

%右端函数

ZZJ=84819.5e-6;      %主动组合链轮的等效转动惯量

CZJ=84819.5e-6;    %从动组合链轮的等效转动惯量

DJZJ=0.0;    %电机转子及连接轴的转动惯量

RG1J=7.01511e-6;  %齿轮1的转动惯量

RG2J=512.148e-6;   %齿轮2的转动惯量

RG3J=53.271e-6;    %齿轮3的转动惯量

RG4J=3889.12e-6;    %齿轮4的转动惯量

RG5J=720.0e-6;     %齿轮5的转动惯量

RG6J=11657.3e-6;     %齿轮6的转动惯量

RI1=38.0/13.0*38.0/13.0*41.0;   %264.8757  !第一级传动的速比

RI2=38.0/13.0*41.0;  %90.61538   !第二级传动的速比

RI3=41.0;       %31.0   !第三级传动的速比

E1=0.95;  %第一级传动的效率

E2=0.95;  %第二级传动的效率

E3=0.47;  %第三级传动的效率

%CDJ-传动部分的总的等效转动惯量

CDJ=ZZJ+CZJ+DJZJ*RI1^2*E1*E2*E3*2.0+RG1J*RI1^2*E1*E2*E3+RG2J*RI2^2*E2*E3+RG3J*RI2^2*E2*E3+RG4J*RI3^2*E3+RG5J*RI3^2*E3+RG6J;

%DLIJ-动力矩

x=6400;

DLIJ=(x-6722.0)/(5828.0-6722.0)*0.75+(x-5828.0)/(6722.0-5828.0)*0.5;

a=EJ+CDJ;

b=DJDF;

c=EF+EFMDT+DLIJ;

dy(1)=y(2);

dy(2)=-b*(y(2))^2/2/a+c/a;

调用四阶龙格-库塔法解微分方程：

%rk4解弹仓运动方程