经典四阶 Runge-Kutta 方法

本文介绍如何计算经典四阶 Runge-Kutta 数值解及其收敛阶，并用 MATLAB 编程，最后展示了输出结果。

经典四阶 Runge-Kutta 方法

$$\left\{ \begin{aligned} &u_{n+1}=u_{n}+\frac{h}{6}\left(k_{1}+2 k_{2}+2 k_{3}+k_{4}\right), \\ &k_{1} =f\left(t_{n}, u_{n}\right), \\ &k_{2} =f\left(t_{n}+\frac{h}{2}, u_{n}+\frac{h}{2}k_{1}\right), \\ &k_{3} =f\left(t_{n}+\frac{h}{2}, u_{n}+\frac{h}{2}k_{2}\right), \\ &k_{4} =f\left(t_{n}+h, u_{n}+hk_{3}\right). \end{aligned} \right.$$

线性 ODE 例子

$$\left\{\begin{aligned} &\frac{\mathrm{d}u}{\mathrm{d}t}=t^2+t-u ,\quad 0 < t \leqslant 1, \\ &u(0)=0. \end{aligned} \right.$$

方程的真解：$u(x)=−e^{−t}+t^2−t+1$.

计算数值解

MATLAB 程序

% RungeKutta.m
% Runge-Kutta method for the ODE model
% u'(x)=x^2+x-u, x in [0,1] 
% Initial condition: u(0)=0 ;
% Exact solution: u(x)=-exp(-x)+x^2-x+1.
clear all; close all;
h=0.1;
x=0:h:1;                     % interval partition
N=length(x)-1;
u(1)=0;                      % initial value
fun=@(x,u) x.^2+x-u;         % RHS
for n=1:N
    k1=fun(x(n),u(n));
    k2=fun(x(n)+h./2,u(n)+h.*k1/2);
    k3=fun(x(n)+h./2,u(n)+h.*k2/2);
    k4=fun(x(n)+h,u(n)+h.*k3);
    u(n+1)=u(n)+h.*(k1+2.*k2+2.*k3+k4)./6;
end
ue=-exp(-x)+x.^2-x+1;        % exact solution
plot(x,ue,'b-',x,u,'r+','LineWidth',1)
legend('Exact','Numerical','location','northwest')
% title('Runge-Kutta Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('x','fontsize',16), ylabel('u','fontsize',16)

% print -dpng -r600 RungeKutta.png
% print -depsc2 RungeKutta.eps

输出结果

计算误差阶

MATLAB 程序

% RungeKutta_error.m
% Runge-Kutta method for the ODE model
% u'(t)=t^2+t-u, t \in [0,2] 
% Initial value : u(0)=0 ;
% Exact solution : u(t)=-exp(-t)+t^2-t+1.
clear all; close all;
Nvec=[10 50 100 500 1000];        % Number of partitions
Error=[];
fun=@(t,u) t.^2+t-u;              % RHS
for k=1:length(Nvec)
    N=Nvec(k);
    h=1/N;
    x=0:h:1;                      % interval partition
    u(1)=0;                       % initial value
    for n=1:N
        k1=fun(x(n),u(n));
        k2=fun(x(n)+h./2,u(n)+h.*k1/2);
        k3=fun(x(n)+h./2,u(n)+h.*k2/2);
        k4=fun(x(n)+h,u(n)+h.*k3);
        u(n+1)=u(n)+h.*(k1+2.*k2+2.*k3+k4)./6;
    end
    ue=-exp(-x)+x.^2-x+1;         % exact solution
    error=max(abs(u-ue));
    Error=[Error,error];
end
plot(log10(Nvec),log10(Error),'ro-','MarkerFaceColor','w','LineWidth',1)
hold on
plot(log10(Nvec),log10(Nvec.^(-4)), '--')
grid on
%title('Convergence of Runge-Kutta Method','fontsize',12)
set(gca,'fontsize',12)
xlabel('log_{10}N','fontsize',16), ylabel('log_{10}Error','fontsize',16)

% add annotation of slope
ax = [0.57 0.53];
ay = [0.68 0.63];
annotation('textarrow',ax,ay,'String','slope = -4 ','fontsize',14)

% computating convergence order
for n=1:length(Nvec)-1
    order(n)=-log(Error(n)/Error(n+1))/(log(Nvec(n)/Nvec(n+1)));
end
Error
order

% print -dpng -r600 RungeKutta_error.png
% print -depsc2 RungeKutta_error.eps

输出结果

Error =
    1.0e-05 *
        0.1051    0.0002    0.0000    0.0000    0.0000

order =
    4.0193    4.0057    4.0025    3.9532