Prove that Runge Kutta Method (RK4) is of Order 4

Question

Please somebody help me, recently we have been studying numerical methods for solving ODEs and we went over proofs for the Euler method being order 1 and Huen’s method being order 2. But our lecturer didn’t cover why the Runge Kutta method (RK4) is of fourth order, he set it as an extension exercise but all I get out is that the error is order 3 which is not right. Could some one please show me a correct proof. The method I am defining is: $$k_1 = f(t_0,x(t_0))$$ $$k_2 = f(t_0+h/2,x(t_0)+hk_1/2)$$ $$k_3 = f(t_0+h/2,x(t_0)+hk_2/2)$$ $$k_4 = f(t_0+h,x(t_0)+hk_3)$$ Then $$x(t_0+h)= x(t_0) + h\frac{k_1+2k_2+2k_3+k_4}{6}$$

Thank you in advance!

Answer 1

You can find the order conditions and what they mean in the first two presentations of https://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/, along with various examples up to explicit 4-stage order 4 methods.

Answer 2

The proof is basically taken from section II.2 of Solving Ordinary Differential Equations I by E. Hairer S. P. Nørsett and G. Wanner.

Step 1. Convert the ODE $$ \frac{dx(t)}{dt} = f(t, x(t)) $$ to an autonomous form (the one that does not explicitly depend on the independent variable). This basically is done to simplify the following series expansions.

I am assuming that $x(t) \in \mathbb R^n$. The case $n = 1$ does not really simplify things. Let's introduce another independent variable $s$, equivalent to the variable $t$. The ODE can be now rewritten as $$ \begin{aligned} \frac{dt(s)}{ds} &= 1\\ \frac{dx(s)}{ds} &= f(t(s), x(s)) \end{aligned} $$ Now the right hand side depends only on $t(s)$ and $x(s)$, but not on the $s$ itself. To simplify it further, let's introduce $$ X(s) = \begin{pmatrix} t(s)\\ x(s) \end{pmatrix} \in \mathbb R^{n+1}, \qquad F(X) = \begin{pmatrix} 1\\ f(X) \end{pmatrix} \in \mathbb R^{n+1}, \qquad $$

$$ \frac{dX(s)}{ds} = F(X(s)). $$

It's not hard to verify that RK4 formula can also be rewritten as $$ \begin{aligned} K_1 &= F(X(s_0))\\ K_2 &= F(X(s_0) + h K_1 / 2)\\ K_3 &= F(X(s_0) + h K_2 / 2)\\ K_4 &= F(X(s_0) + h K_3)\\ X(s_0 + h) &= X(s_0) + h \frac{K_1 + 2 K_2 + 2 K_3 + K_4}{6} \end{aligned} $$

Step 2. To show that the method has the fourth order we need to verify that $$ Y(s_0 + h) = Y(s_0) + h \frac{K_1 + 2 K_2 + 2 K_3 + K_4}{6} + O(h^5). $$ Here $X(s)$ is the numerical solution, obtained from the RK4 procedure, and $Y(s)$ is the true solution of the ODE and $X(s_0) = Y(s_0)$. The truncation error is one order of magnitude smaller, since it is a local (one step) error, and the global error is $O(h^{-1})$ (number of steps) times bigger.

All is need to be done now is to simply plug the Taylor's expansion to the every value in the right and the left hand side.

The left side is straightforward: $$ Y(s_0 + h) = Y(s_0) + h Y'(s_0) + \frac{h^2}{2} Y''(s_0) + \frac{h^3}{6} Y'''(s_0) + \frac{h^4}{24} Y^{IV}(s_0) + O(h^5). $$ Since $Y(s)$ is the solution of the ODE, we can plug $$ Y'(s_0) = F(Y(s_0)). $$ Moreover, $$ Y''(s_0) = \frac{\partial F}{\partial Y}(s_0) Y'(s_0) = \frac{\partial F}{\partial Y}(s_0) F(Y(s_0)). $$ Recall that $F : \mathbb R^{n+1} \to \mathbb R^{n+1}$ and care should be taken while differentiation, since $\frac{\partial F}{\partial Y}$ is a matrix and cannot be easily swapped with $F$, which is vector. Differentiating further will introduce $\frac{\partial^2 F}{\partial Y^2}$ derivative which is an $(1,2)$ tensor. Personally I recommend using the Einstein's summation rule to minimize errors while using the chain rule to differentiate: $$ \begin{aligned} (Y')_i &= F_i\\ (Y'')_i &= \frac{\partial F_i}{\partial Y_j} F_j\\ (Y''')_i &= \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j F_k + \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} F_k \\ (Y^{IV})_i &= \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} \frac{\partial F_j}{\partial Y_l} F_l F_k + \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j \frac{\partial F_k}{\partial Y_l} F_l + \frac{\partial^2 F_i}{\partial Y_j \partial Y_l} F_l \frac{\partial F_j}{\partial Y_k} F_k + \frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_l F_k + \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} \frac{\partial F_k}{\partial Y_l} F_l = \\ &= \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + 3\frac{\partial^2 F_i}{\partial Y_j \partial Y_k} \frac{\partial F_j}{\partial Y_l} F_l F_k + \frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_l F_k + \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} \frac{\partial F_k}{\partial Y_l} F_l. \end{aligned} $$ Those differentials look really horrible, so the tree notation was invented to represent each elementary differential like $\frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_l F_k$ as a labeled tree with indices. The main purpose is again to simplify the notation.

So far the left hand side is done. It's time to work with the right hand side. Let's start with $K_1$. That's pretty easy, $$ K_1 = F(Y(s_0)). $$ Now $K_2$: $$ \begin{aligned} (K_2)_i &= F_i(Y_j(s_0) + h/2 (K_1)_j) = F_i + \frac{h}{2} \frac{\partial F_i}{\partial Y_j} (K_1)_j + \frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} (K_1)_j (K_1)_k + \frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} (K_1)_j (K_1)_k (K_1)_l + O(h^4) = \\ &= F_i + \frac{h}{2} \frac{\partial F_i}{\partial Y_j} F_j + \frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j F_k + \frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + O(h^4) \end{aligned} $$ Note the same kind of differentials in the expansion. Expanding $K_3$ is a bit trickier. Basically it's the same as $K_2$: $$ (K_3)_i = F_i(Y_j(s_0) + h/2 (K_2)_j) = F_i + \frac{h}{2} \frac{\partial F_i}{\partial Y_j} (K_2)_j + \frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} (K_2)_j (K_2)_k + \frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} (K_2)_j (K_2)_k (K_2)_l + O(h^4), $$ but we cannot just plug $F$ instead of $K_2$, like we've done for the expansion of $K_2$. We need to plug the whole expression for the $K_2$ now! This work can be simplified by an observation: for the subexpression $\frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} (K_2)_j (K_2)_k (K_2)_l$ we need only $O(h)$ correct expression for $K_2$, any extra terms will anyway go into the last $O(h^4)$ term, since the subexpression has leading $\frac{h^3}{48}$ multiplier. Similarly, for the $\frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} (K_2)_j (K_2)_k$ we need $O(h^2)$ correct expression for $K_2$ and for the $\frac{h}{2} \frac{\partial F_i}{\partial Y_j} (K_2)_j$ — only $O(h^3)$ correct expression.

$$ \begin{aligned} (K_3)_i = F_i + \frac{h}{2} \frac{\partial F_i}{\partial Y_j} \left( F_j + \frac{h}{2} \frac{\partial F_j}{\partial Y_k} F_k + \frac{h^2}{8} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_k F_l \right) + \frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} \left( F_j + \frac{h}{2} \frac{\partial F_j}{\partial Y_l} F_l \right) \left( F_k + \frac{h}{2} \frac{\partial F_k}{\partial Y_l} F_l \right) + \frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + O(h^4) = \\ = F_i + \frac{h}{2} \frac{\partial F_i}{\partial Y_j} F_j + \frac{h^2}{4} \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} F_k + \frac{h^2}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j F_k + \frac{h^3}{16} \frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_k F_l + \frac{h^3}{8} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j \frac{\partial F_k}{\partial Y_l} F_l + \frac{h^3}{48} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + O(h^4). \end{aligned} $$

Finally, expanding $K_4$ gives $$ K_4 = F_i(Y_j(s_0) + h (K_3)_j) = F_i + h \frac{\partial F_i}{\partial Y_j} (K_3)_j + \frac{h^2}{2} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} (K_3)_j (K_3)_k + \frac{h^3}{6} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} (K_3)_j (K_3)_k (K_3)_l + O(h^4). $$

$$ \begin{aligned} K_4 &= F_i + h \frac{\partial F_i}{\partial Y_j} \left( F_j + \frac{h}{2} \frac{\partial F_j}{\partial Y_k} F_k + \frac{h^2}{4} \frac{\partial F_j}{\partial Y_k} \frac{\partial F_k}{\partial Y_l} F_l + \frac{h^2}{8} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_k F_l \right) + \frac{h^2}{2} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} \left(F_j + \frac{h}{2} \frac{\partial F_j}{\partial Y_l} F_l\right) \left(F_k + \frac{h}{2} \frac{\partial F_k}{\partial Y_l} F_l\right) + \frac{h^3}{6} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + O(h^4) = \\ &= F_i + h \frac{\partial F_i}{\partial Y_j} F_j + \frac{h^2}{2} \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} F_k + \frac{h^2}{2} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j F_k + \frac{h^3}{4} \frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} \frac{\partial F_k}{\partial Y_l} F_l + \frac{h^3}{8} \frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_k F_l + \frac{h^3}{2} \frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j \frac{\partial F_k}{\partial Y_l} F_l + \frac{h^3}{6} \frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l + O(h^4). \end{aligned} $$

Step 3. I'm going to make a table with all the expansion coefficients we have:

$$ \begin{array}{c|ccccc} &\frac{Y_i(s_0 + h) - Y_i(s_0)}{h} & (K_1)_i & (K_2)_i & (K_3)_i & (K_4)_i\\\hline F_i & 1 & 1 & 1 & 1 & 1\\ h\frac{\partial F_i}{\partial Y_j} F_j & 1/2 & 0 & 1/2 & 1/2 & 1\\ h^2\frac{\partial^2 F_i}{\partial Y_j \partial Y_k} F_j F_k & 1/6 & 0 & 1/8 & 1/8 & 1/2\\ h^2\frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} F_k & 1/6 & 0 & 0 & 1/4 & 1/2\\ h^3\frac{\partial^3 F_i}{\partial Y_j \partial Y_k \partial Y_l} F_j F_k F_l & 1/24 & 0 & 1/48 & 1/48 & 1/6\\ h^3\frac{\partial^2 F_i}{\partial Y_j \partial Y_k} \frac{\partial F_j}{\partial Y_l} F_l F_k & 3/24 & 0 & 0 & 1/8 & 1/2\\ h^3\frac{\partial F_i}{\partial Y_j} \frac{\partial^2 F_j}{\partial Y_k \partial Y_l} F_l F_k & 1/24 & 0 & 0 & 1/16 & 1/8\\ h^3\frac{\partial F_i}{\partial Y_j} \frac{\partial F_j}{\partial Y_k} \frac{\partial F_k}{\partial Y_l} F_l & 1/24 & 0 & 0 & 0 & 1/4 \end{array} $$

It's now easy to see that each term of $\frac{Y_i(s_0 + h) - Y_i(s_0)}{h}$ cancels with a corresponding term of $\frac{K_1 + 2K_2 + 2K_3 + K_4}{6}$, except for the $O(h^4)$, which means that the method is indeed of the fourth order.

PS. Personally I find this manual proof extremely error prone and useless. For RK methods it has been developed a framework of representation notations including elementary differential trees along with corresponding manipulation rules (c.f. Faa di Bruno's Formula). Using this technique and modern computer algebra systems the process of deriving order conditions for RK methods can be greatly simplified.