请问用四阶龙格库塔法解二阶微分方程的思想是什么?

问题描述:

请问用四阶龙格库塔法解二阶微分方程的思想是什么?
最近我遇到了一个难题!就是求一个二阶微分方程,形式如:y''=f(y),初始条件是y(0)=0,y'(0)=0,y是t的函数.要求用四阶龙格库塔法求解!请问思路是什么?由于f(y)函数复杂,所以解析解是求不出来的,必须用四阶龙格库塔法解之,如果您有实例,希望能够和我分享一下!

令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;

while err

     X=[0 0];

     k=xishu*k;

     K=[K;k];

     Y=[];

    for i=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];

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

        break

    elseif ymax>92/rad2deg(1)

        xishu=0.98;

    else

        xishu=1.01;

    end

end

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

hold on

grid on

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

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

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

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

title(tstring,'fontsize',14);

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


function Y = 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);


function dy=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)));