一、计算公式

对于形如以下的常微分方程：

采用四阶龙格库塔法的计算公式：

四阶 龙格库塔法精度为4，属于单步递推法。

单步递推法的基本思想是从（x(i),y(i)）点出发，以某一斜率沿直线达到（x(i+1),y(i+1)）点。

二、实例计算

从上述定义可以看出，龙格库塔实质上是求一阶微分方程，但是如果将一阶导看作变量，则二阶导也不过是这个变量的一阶导而已。

接下来就看实例吧：

对于下述二阶方程：

1、基本思想：

令位移为

q的一阶导，即位移的一阶导（速度）为

q的二阶导

求解位移u(1)的龙格库塔计算方法如下：

KK1=u(2);KK2=u(2)+h/2*KK1;KK3=u(2)+h/2*KK2;KK4=u(2)+h*KK3;u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);

求解速度u(2)的龙格库塔计算方法如下：

K1=-2*eptheton*u(2)-u(1)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);K2=-2*eptheton*(u(2)+h/2*K1)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);K3=-2*eptheton*(u(2)+h/2*K2)-(u(1)+h/2)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);K4=-2*eptheton*(u(2)+h*K3)-(u(1)+h)+Fm+Fah*wh*wh*sin(wh*tao+fav_h);u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);

2、编程实现

参数设置：

h=0.001; % 步长t0=0:h:400; % 总时长w=5;ep=0.02;Fm=0.1;Fah=0.05;u(1)=0;u(2)=0; % 初值

总的程序实现

h=0.001;t0=0:h:400;w=5;ep=0.02;Fm=0.1;Fah=0.05;u(1)=0;u(2)=0;for i=1:length(t0) % 进行多次迭代tao=t0(i);u=RK(u,tao,h,ep,w,Fm,Fah);Result(i,:)=u; % 将每次迭代的位移和速度保存endfigure(1)subplot(2,1,1)plot(t0,Result(:,1))% 绘制位移图xlabel('Time')ylabel('displacement')subplot(2,1,2)plot(t0,Result(:,2))% 绘制速度图xlabel('Time')ylabel('velocity')function u=RK(u,tao,h,ep,w,Fm,Fah)KK1=u(2);KK2=u(2)+h/2*KK1;KK3=u(2)+h/2*KK2;KK4=u(2)+h*KK3;u(1)=u(1)+h/6*(KK1+2*KK2+2*KK3+KK4);K1=-2*ep*u(2)-u(1)+Fm+Fah*cos(w*tao);K2=-2*ep*(u(2)+h/2*K1)-u(1)-h/2+Fm+Fah*cos(w*tao);K3=-2*ep*(u(2)+h/2*K2)-u(1)-h/2+Fm+Fah*cos(w*tao);K4=-2*ep*(u(2)+h*K3)-u(1)-h+Fm+Fah*cos(w*tao);u(2)=u(2)+h/6*(K1+2*K2+2*K3+K4);end

结果图如下所示

值得特别注意的是

u=RK(u,tao,h,ep,w,Fm,Fah);

输入的u与输出的u一定要符号一致，从而确保下一次迭代的初始值是上一次的值。确保迭代能一直进行下去。

三、辅助验证

接下来用MATLAB自带的ode45函数来进行验证。

之前已经写过ode45函数的用法，在此不再进行介绍。

源程序如下：

t0=0:0.001:400;w=5;ep=0.02;Fm=0.1;Fah=0.05;y0=[0 0];[t,u]=ode45(@(t,u) odefun(t,u,w,ep,Fm,Fah),t0,y0);figure(1)subplot(2,1,1)plot(t,u(:,1))xlabel('Time')ylabel('displacement')subplot(2,1,2)plot(t,u(:,2))xlabel('Time')ylabel('velocity')function du=odefun(t,u,w,ep,Fm,Fah)du=[u(2);-2*ep*u(2)-u(1)+Fm+Fah*cos(w*t)]; end

运算结果图如下所示

两中方法计算的结果是一样的。

创作不易，如有帮助，记得点个赞呐。