令y1=y';y2=y''=(y1)';得到了二阶微分方程的一阶形式,然后可以按照普通的一节微分方程组的四阶龙格库塔法求解了；这样说的有点抽象,你可以把方程复制出来我有空帮你编写一下

　　好的非常感谢呢

　　292200y''=k*(90-y)-8390.7*sin(y);y(0)=0,y'(0)=0;

　　其中y是时间t的函数，k是自己可以调的一个系数；步长取0.001吧，最少计算一千次，因为这是折叠机翼的转动微分方程，要求在1秒内展开，所以最少计算一千次。解这个微分方程的目的是选取k，从而实现1秒内展开，即y（t=1）>=90度

　　非常感谢哈

　　默认y的单位是弧度

　　k=1000;

　　t=0:0.001:1;

　　Y=[];

　　err=1

　　K=[];

　　Ymax=[];

　　xishu=1.01;

　　whileerr

　　X=[00];

　　k=xishu*k;

　　K=[K;k];

　　Y=[];

　　fori=1:1001

　　Y_1=Runge_Kutta41(t(i),X,@folded_wing,0.001,k);

　　Y=[Y;t(i),Y_1];

　　X=Y_1;

　　end

　　ymax=max(Y(:,2));

　　Ymax=[Ymax;ymax];

　　ifymax>=pi/2&&ymax<=92/rad2deg(1)

　　break

　　elseifymax>92/rad2deg(1)

　　xishu=0.98;

　　else

　　xishu=1.01;

　　end

　　end

　　plot(Y(:,1),Y(:,2),'-b','linewidth',2);

　　holdon

　　gridon

　　xlabel('t-time','fontsize',14)

　　ylabel('y(rad)','fontsize',14)

　　plot([0,1],[pi/2pi/2],'-g','linewidth',2);

　　tstring=['292200y"=',num2str(k),'*(90-y)-8390.7*sin(y);'];

　　title(tstring,'fontsize',14);

　　legend('y"','y=pi/2')

　　functionY=Runge_Kutta41(t,X,f,h,k)

　　K1=f(t,X,k);

　　K2=f(t+h/2,X+h/2*K1,k);

　　K3=f(t+h/2,X+h/2*K2,k);

　　K4=f(t+h,X+h*K3,k);

　　Y=X+h/6*(K1+2*K2+2*K3+K4);

　　functiondy=folded_wing(t,y,k)

　　%y1=y;y2=y'

　　dy=zeros(1,2);

　　dy(1)=y(2);

　　dy(2)=1/292200*(k*(90-y(1))-8390.7*sin(y(1)));