标准四阶Runge-Kutta公式在一般的数值计算方法中都有，其一元常微分方程初值问题对应形式为

$\begin{align*}\label{eq:odeipv1rk44} y_{n+1} &= y_n+\dfrac{h}{6} (k_1+2k_2+2k_3+k_4),\\ k_1 &= f(x_n,y_n),\\ k_2 &= f(x_n+\dfrac{1}{2}h,y_n+\dfrac{1}{2}h\,k_1)\\ k_3 &= f(x_n+\dfrac{1}{2}h,y_n+\dfrac{1}{2}h\,k_2)\\ k_4 &= f(x_n+h,y_n+ h\,k_3)\\ \end{align}$

下面结合例子给出MATLAB程序和结果。

先看一个一元常微分方程初值问题数值解。

$\dfrac{dy}{dx}=x(1-y),0\leq x\leq1,y(0)=1.$

MATLAB程序如下

odefun = @(x,y)(10*x*(1-y)); %微分方程匿名函数，好处是不用再另写一个函数文件

h = 0.1;  x = 0:h:1; y = x; y(1)=0;%参数初值

for n = 1:length(x)-1

       k1 = odefun(x(n),y(n));

       k2 = odefun(x(n)+h/2,y(n)+h/2*k1);

       k3 = odefun(x(n)+h/2,y(n)+h/2*k2);

       k4 = odefun(x(n)+h,y(n)+h*k3);

       y(n+1) = y(n)+h/6*(k1+2*k2+2*k3+k4);%四级四阶Runge-Kutta公式

end

plot(x,y,'.',x,1-exp(-5*x.^2),'r-.d') %绘制数值解与解析解

legend('数值解','解析解','Location','northwest')

运行结果为

再看一个复杂点的例子,先看一个二阶常微分方程

$y''-2y'+2y = e^{2t}\sin t,0\leqslant t\leqslant1,y(0)=-0.4,y'(0)=-0.6,$

先化为方程组

$y_1'=y_2,y_2' = 2y_2-2y_1+e^{2t}\sin t,$

标准四阶Runge-Kutta公式在一般的数值计算方法中都有，其向量形式为

$\begin{align*}\label{eq:odeipv1rk44v} \textbfsymbol{y}_{n+1} &= \textbfsymbol{y}_n+\dfrac{h}{6} (\textbfsymbol{k}_1+2\textbfsymbol{k}_2+2\textbfsymbol{k}_3+\textbfsymbol{k}_4),\\ \textbfsymbol{k}_1 &= \textbfsymbol{F}(x_n,\textbfsymbol{y}_n),\\ \textbfsymbol{k}_2 &= \textbfsymbol{F}(x_n+\dfrac{1}{2}h,\textbfsymbol{y}_n+\dfrac{1}{2}h\,k_1)\\ \textbfsymbol{k}_3 &= \textbfsymbol{F}(x_n+\dfrac{1}{2}h,\textbfsymbol{y}_n+\dfrac{1}{2}h\,k_2)\\ \textbfsymbol{k}_4 &= \textbfsymbol{F}(x_n+h,\textbfsymbol{y}_n+ h\,k_3)\\ \end{align}$

上面常微分方程组问题的求解程序为

%%%%%%%%%%%%%%%%% 二阶ODE %%%%%%%%%%

%% 二阶常微分方程 y''-2*y'+2*y = exp(2*t)*sin(t),0<=t<=1;y(0) = -0.4,y'(0) = -0.6

odefun2 = @(t,y)[y(2); -2*y(2)+2*y(1)+exp(2*t)*sin(t)]; % 化二阶常微分方程为二阶常微分方程组

t0 = 0; tf = 1; h = 0.1;  t = t0:h:tf; y = zeros(length(t),2); y(1,:)=[-0.4,-0.6];%参数初值

for n = 1:length(t)-1

       k1 = odefun2(t(n),y(n,:)');

       k2 = odefun2(t(n)+h/2,y(n,:)'+h/2*k1); % 这里要特别注意同型向量相加

       k3 = odefun2(t(n)+h/2,y(n,:)'+h/2*k2);

       k4 = odefun2(t(n)+h,y(n,:)'+h*k3);

       y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';%四级四阶Runge-Kutta公式

end

plot(t,y(:,1),'.',t,y(:,2),'r-.') %绘制数值解

legend('y(t)','Dy(t)','Location','northwest')

hold on

[t,y] = ode45(odefun2,[t0:h:tf],[-0.4,-0.6]);

plot(t',y(:,1),'k-',t,y(:,2),'k-.') %绘制ode45得到的数值解

结果为

再看个复杂点的滑动摆例子

这个滑动摆满足微分方程组

$\dfrac{d^2\theta}{d t^2} = \dfrac{-M\cos\theta \sin\theta\left( \dfrac{d \theta}{d t}\right)^2-\dfrac{g}{l}\sin \theta}{1-M\cos^2\theta}, \dfrac{d^2x}{d t^2} = \dfrac{ Mg\cos\theta \sin\theta +Ml\sin\theta\left( \dfrac{d \theta}{d t}\right)^2}{1-M\cos^2\theta}.$

先降阶，令 $y_1 = \theta,y_2 = \dfrac{d \theta}{d t},y_3 = x,y_4 = \dfrac{d x}{d t}$, 则上式可化为四个一阶微分方程：

$\dfrac{d y_1}{d t} = y_2, \dfrac{d y_2}{d t} = \dfrac{-My_2^2\cos y_1 \sin y_1 -\dfrac{g}{l}\sin y_1}{1-M\cos^2y_1},$

$\dfrac{d y_3}{d t} = y_4, \dfrac{d y_4}{d t} = \dfrac{ Mg\cos y_1 \sin y_1 +Mly_2^2\sin y_1 }{1-M\cos^2 y_1}.$

编写MATLAB程序为

g = 9.8; m1 = 3; m2 = 3; l =1; M = m1/(m1+m2);%%参数设置

odefunhdb = @(t,y)[y(2); (-M*sin(y(1))*cos(y(1))*y(2)^2-g/l*sin(y(1)))/(1-M*cos(y(1))^2);...

y(4); (M*g*sin(y(1))*cos(y(1))+M*l*sin(y(1))*y(2)^2)/(1-M*cos(y(1))^2)];

%%%% 用四级四阶Runge-Kutta公式求解 %%%%%%%

h = 0.05;  t = (0:h:5.5)'; y = zeros(length(t),4); y(1,:)=[pi/4,0,-cos(pi/4)/2,0];%参数设置与初值

for n = 1:length(t)-1

       k1 = odefunhdb(t(n),y(n,:)');

       k2 = odefunhdb(t(n)+h/2,y(n,:)'+h/2*k1);

       k3 = odefunhdb(t(n)+h/2,y(n,:)'+h/2*k2);

       k4 = odefunhdb(t(n)+h,y(n,:)'+h*k3);

       y(n+1,:) = y(n,:)+h/6*(k1+2*k2+2*k3+k4)';

end

[t1, y1] = ode45(odefunhdb, [0:0.05:5.5], [pi/4,0,-cos(pi/4)/2,0]); %微分方程组

% 画曲线

figure('color',[1 1 1],'unit','normalized','position',[0.1 0.1 0.45 0.5])

[ax,h1,h2] = plotyy(t,y(:,3),t,180*y(:,1)/pi);      % 双坐标系

xlabel('时间{it t}','FontName','Times New Roman','FontSize',20)

ylabel(ax(1),'滑块位置{it x}','FontName','Times New Roman','FontSize',20)

ylabel(ax(2),'滑块位置{it theta}','FontName','Times New Roman','FontSize',20)

set(ax(1),'FontName','宋体','FontSize',20,'xlim',[0, 5.5])

set(ax(2),'FontName','宋体','FontSize',20,'xlim',[0, 5.5],'ylim',[-60,60])

set(h1,'LineStyle','-','Marker','*','LineWidth',2)

set(h2,'LineStyle','-.','Marker','.','LineWidth',2)

  这种方法在进一步的数值测验中表现并不太好，所以一般采用其他算法，比如MATLAB自带的 ode45( )。