Comparing truncation errors from RK4

Consider the ODE $y'=f(x,y)$ with $x_0 = 0$, $y_0 = y(x_0) = 0$. We wish to approximate $y$ on the interval $[0,1]$. Let $h_1$ be some reasonable step size and $h_2 = \frac{1}{2} h_1$. These two choices create two partitions $P_i = \{ 0, h_i, 2h_i, ..., 1 \}$ of the interval $[0,1]$ ($i=1,2$). Denote by $\tilde{y}_j(h_i)$ the $j$:th RK4 approximation of $y(jh_i)$ ($j=0,1,...$).

I have an assignment to verify that RK4 has a global truncation error of $O(h^4)$ (for a specific ODE) but I'm not sure how to do that. My best guess would be to calculate $$\frac{y(jh_2)-\tilde{y}_j(h_2)}{y(jh_1)-\tilde{y}_j(h_1)}$$ and make the argument that this is close to $1/32 = 2^{-5}$ whence GTE is $O(h^{5-1})$, but I'm not sure at all and the approximations I have don't seem to support that assertion anyway.

• Have you resolved this question or is it still open/relevant? – Carl Christian May 10 '16 at 16:57
• @Carl Still open! :) – Erik Vesterlund May 11 '16 at 0:24
• There are some MATLAB codes here www8.cs.umu.se/kurser/5DV005/HT15/Programs/TOPICS which might be of use to you. In particular, /ODE/rode.m, /INT/rint.m, /DIF/rdif.m which does RE for ODE, integration and differentiation. All codes come with minimal working examples (MWE) and the numbers I quoted below were generated by a slight variation of the MWE for rode.m. The codes are only intended for teaching students some principles of numerical analysis. Use them at you own risk! – Carl Christian May 13 '16 at 8:13

Let $\Sigma_h$ denote the grid $$\Sigma_h = \{ jh \: : \: j \in \{0,1,2,\dotsc,\} \}$$ and let $v_h : \Sigma_h \rightarrow \mathbb{R}$ denote the grid function obtained by applying RK4 with stepsize $h$ to your initial value problem. Let $u : I \rightarrow \mathbb{R}$ denote the exact solution of your initial value problem. You intuition is to compute the fractions $$G_h(t) = \frac{u(t) - v_{2h}(t)}{u(t) - v_h(t)}, \quad t \in \Sigma_h \cap \Sigma_{2h} = \Sigma_{2h},$$ and one would hope to find that $G_h(t) \approx 16 = 2^4$. Unfortunately, this is not case and the numbers are in fact quite erratic. Instead, one should consider Richardson's fractions $F_h$ given by $$F_h(t) = \frac{v_{2h}(t) - v_{4h}(t)}{v_h(t) - v_{2h}(t)}, \quad t \in \Sigma_{4h}$$ Provided that the is sufficient differentiability available, you will find that these fractions tend (monotonically) to $2^p$, where $p$ is the order of the method. This nice behavior relies on the existence of an asymptotic error expansion of the form $$u(t) - v_h(t) = \alpha(t) h^p + \beta(t) h^{q} + O(h^r), \quad p < q < r.$$ Deriving the expansion is typically a bit complicated, but there is nothing stopping you from verifying numerically that Richardson's fractions behave as described.

As an illustration of this technique consider the ODE $$u'(t) = - \frac{1}{10} u(t), \quad t \in [0,1]$$ together with the initial condition $y(0) = 2$. The exact solution is of course $u(t) = 2 \exp(-t/10)$. This equation was integrated using a second order explicit Runge-Kutta method and several different stepsizes $h_j = 2^{j} h_0$. The results are given here


t          v(h)           F( 4h)          F( 2h)          F( 1h)          est              err
---------------------------------------------------------------------------------------------------------
0.000000e+00  2.00000000e+00             NaN             NaN             NaN  0.00000000e+00   0.00000000e+00
1.000000e-01  1.98009967e+00  4.00350138e+00  4.00174895e+00  4.00090553e+00 -5.15760027e-11  -5.15707477e-11
2.000000e-01  1.96039735e+00  4.00350145e+00  4.00174932e+00  4.00088998e+00 -1.02126011e-10  -1.02114761e-10
3.000000e-01  1.94089107e+00  4.00350134e+00  4.00174793e+00  4.00089746e+00 -1.51664533e-10  -1.51647805e-10
4.000000e-01  1.92157888e+00  4.00350125e+00  4.00174871e+00  4.00089391e+00 -2.00207406e-10  -2.00184758e-10
5.000000e-01  1.90245885e+00  4.00350138e+00  4.00174932e+00  4.00088602e+00 -2.47769583e-10  -2.47740273e-10
6.000000e-01  1.88352907e+00  4.00350137e+00  4.00174969e+00  4.00088205e+00 -2.94365347e-10  -2.94329894e-10
7.000000e-01  1.86478764e+00  4.00350144e+00  4.00174958e+00  4.00088337e+00 -3.40008984e-10  -3.39968276e-10
8.000000e-01  1.84623269e+00  4.00350147e+00  4.00174973e+00  4.00088345e+00 -3.84715223e-10  -3.84669407e-10
9.000000e-01  1.82786237e+00  4.00350140e+00  4.00175046e+00  4.00088093e+00 -4.28498348e-10  -4.28446834e-10
1.000000e+00  1.80967484e+00  4.00350139e+00  4.00175007e+00  4.00088449e+00 -4.71371534e-10  -4.71316319e-10


The last column contains the actual error, the second to last column contains Richardsons error estimates $$E_h(t) = \frac{v_h(t) - v_{2h}(t)}{2^p-1},$$ where $p=2$ is the order. You will notice that there is exceedingly good agreement between the actual error and the error estimate. However, here you should primarily be interested in Richardson's fraction. For each row, i.e. for each fixed value of $t$ you will observe that the fraction is tending monotonically to $4 = 2^2$ as you move to the right (decreasing the value of the stepsize). This is experimental verification of the fact that there exists an asymptotic error expansion of the given form and that the dominant error term, i.e $\alpha(t) h^p$ has $p = 2$.

Richardson's full name is Lewis Fry Richardson. Richardson extrapolation is discussed in almost any introduction to numerical analysis, but it primarily used to compute error estimates rather than probe a "new" method to experimentally determine the order as I have done here. It is possible to say substantially more about this. In particular, one needs to be careful when using high order methods and small stepsizes as there are numerical issues related to the limitations of floating point arithmetic (catastropic cancellation) which will affect the behaviour of the fractions.

• This is great stuff, I'll have to look into it more closely this weekend, but thank you for taking the time! Regarding Richardson extrapolation I've only seen it used with regards to integration, in neither of the two introductory books I have is it mentioned together with differential equations. Interesting... – Erik Vesterlund May 13 '16 at 0:01
• Thank you very much. In order to solve ODEs more rapidly, people want to use adaptive step size control and increase the step size when ever possible, i.e. when the solution is changing slowly. The mechanism for selecting the new step size hinges on the existence of suitable asymptotic error expansions. They get a good approximation rapidly, but lose the very accurate error estimate which requires a constant step size. This is probably the reason why most textbooks will discuss RE in the context of integration. – Carl Christian May 13 '16 at 7:15