Runge-Kutta order of accuracy

Question

A generic RK method could be written in the following form: \begin{equation} \begin{cases} u^{[n+1]}=u^{[n]}+\Delta t F(t_n, u^{[n]})\\ F(t_n, u^{[n]})=\sum\limits_{i=1}^s{b_iK_i}\\ K_i=f\left(t_n+c_i\Delta t, u^{[n]}+\Delta t\sum\limits_{j=1}^s{a_{ij}K_j}\right) \end{cases}. \end{equation} The method is said to have order $p$ if the local truncation error is of order $p+1$: \begin{equation} E(\Delta t)=\big{\vert} y_{exact}(t+\Delta t)-y_{numeric}(t+\Delta t)\big{\vert}\sim O(\Delta t^{p+1}). \end{equation} To achieve a second order method the following conditions on the coefficients of the method must be verified: \begin{cases} \sum\limits_i{b_i}=1\\ \sum\limits_{j}{a_{ij}}=c_i \\ \sum\limits_{i}{b_ic_i} =\frac{1}{2}\ \end{cases} Suppose the coefficients of a particular method are given by the following Butcher's tableau: \begin{array}{c|ccc} \renewcommand{\arraystretch}{1.5} \gamma & \gamma & 0 \\ 1 & 1-\gamma & \gamma \\ \hline & 1-b_2 & b_2 \\ \end{array} The order conditions are satisfied if \begin{equation} b_2=\frac{1}{2}\frac{2\gamma -1}{\gamma -1} \end{equation} I tried to verify this by using a Matlab routine that solves a test equation for $t\in [0.t_{end}]$ using the method given above. By doubling the number of grid points (i.e. dividing the step size by 2) at each iteration, one can estimate the order of accuracy by using the fact that \begin{equation} \frac{E(2\Delta t)}{E\left(\Delta t\right)}= 2^{p+1}\ \ \ \ \Rightarrow\ \ \ \ \frac{1}{\log(2)}\log\left(\frac{E(2\Delta t)}{E\left(\Delta t\right)}\right)=p+1 \end{equation}. The code I used is the following:

clear 
close all

%Set ending time of the evaluation
tend=2.5;

%Initial condition
CI=1;

%Values for gamma, b2
g=0.2; %gamma
b2=0.5.*(2*g-1)./(g-1);

%Set initial number of grid points
ngrid_init=10;

%Evaluate exact solution at tend
exact_solution_tend=exp(-tend^2/2).*cos(4*tend);

%Set total number of grid points multiplications
nmax=8;

%% Accuracy evaluation

%Set error, accuracy and ngrid_points vectors to zero each loop
%iteration
err=zeros(1,nmax);
npoints=zeros(1,nmax);
ngrid=ngrid_init;

%Loop on different number of grid points to evaluate the error at the
%last point 
for j=1:nmax
    t=linspace(0,tend,ngrid);
    [~,u]=SDIRK_SolverND(@test_solution,g,b2, t, CI);
    err(j)=abs(u(end)-exact_solution_tend);
    npoints(j)=ngrid;
    ngrid=ngrid*2;
end

accuracy=log(err(1:end-1)./err(2:end))./log(2)

This is the test Cauchy problem used:

function f=test_solution(t,u)
     f=-t.*u-4*sin(4*t).*exp(-t.^2/2);

The output of the script is the following vector:

accuracy=[2.2769, 2.0867, 2.0338, 2.0148, 2.0070, 2.0034, 2.0016]

This means that the order o the method is $p=1$ and not $p=2$ as expected. I have no idea why this is happening; can someone give me a hint?

Answer 1

You have to differentiate between the local truncation error and the global error that is composed of all those local errors in the fashion of a compounding saving scheme with unit interest rate $L$ (Lipschitz constant), $e_{n+1}\le e^{Lh}e_n+Ch^{p+1}$ so that in first order $e_{n+1}+\frac CLh^p\le e^{Lh}(e_n+\frac CLh^p)$.

Thus the global error is of the form $e_N\le\frac{C}{L}(e^{Lt_{end}}-1)h^p$, $t_{end}=t_0+Nh$, which for $N=1,~t_{end}=h$ gives again the local error order $O(h^{p+1})$ and for $t_{end}\gg h$ the global error order $O(h^p)$.

Thus the order $2$ that you computed is $p$, as is correct for the method.