Recently, I have been taking a course on ODEs and learning Runge-Kutta methods. To be specific, the 4th order Runge-Kutta method on solving initial value problems.
My instructor and the textbook told me the formula but didn't say anything about the thoughts behind the method. I wrote some code and found that the Runge-Kutta method does perform better than the Euler method, but I can't understand why.
Is anyone willing to give me a hand on how to get the formula of Runge-Kutta? Thanks!
Best Answer
On the history
See Butcher: A History of the Runge-Kutta method
In summary, people (Nystroem, Runge, Heun, Kutta,...) at the end of the 19th century experimented with success in generalizing the methods of numerical integration of functions in one variable $$\int_a^bf(x)dx,$$ like the Gauss, trapezoidal, midpoint and Simpson methods, to the solution of differential equations, which have an integral form $$y(x)=y_0+\int_{x_0}^x f(s,y(s))\,ds.$$
Carl Runge in 1895[1] came up with ("by some curious inductive process" - "auf einem eigentümlich induktiven Wege" wrote Heun 5 years later) the 4-stage 3rd order method \begin{align} k_1&=f(x,y)Δx,\\ k_2&=f(x+\tfrac12Δx,y+\tfrac12k_1)Δx\\ k_3&=f(x+Δx,y+k_1)Δx\\ k_4&=f(x+Δx,y+k_3)Δx\\ y_{+1}&=y+\tfrac16(k_1+4k_2+k_4) \end{align}
[1] "Über die numerische Auflösung von Differentialgleichungen", Math. Ann. 46, p. 167-178
Inspired by this Karl Heun in 1900[2] explored methods of the type $$ \left.\begin{aligned}k^i_m &= f(x+ε^i_mΔx,y+ε^i_mk^{i+1}_m)Δx,~~ i=1,..,s,\\ k^{s+1}_m&=f(x,y)Δx\end{aligned}\right\},~~ m=1,..,n\\ y_{+1}=y+\sum_{m=1}^n\alpha_mf(x+ε^0_mΔx,y+ε^0_mk^1_m)Δx $$ He computed the order conditions by Taylor expansion and constructed methods of this type up to order 4, however the only today recognizable Runge-Kutta methods among them are the order-2 Heun-trapezoidal method and the order 3 Heun method.
[2] "Neue Methode zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen", Z. f. Math. u. Phys. 45, p. 23-38
Wilhelm Kutta in his publication one year later in 1901[3] considered the scheme of Heun wasteful in the number of function evaluations and introduced what is today known as explicit Runge-Kutta methods, where each new function evaluation potentially contains all previous values in the $y$ update. \begin{align} k_1&=f(x,y)Δx,\\ k_m&=f(x+c_mΔx, y+a_{m,1}k_1+...+a_{m,m-1}k_{m-1})Δx,&& m=2,...,s\\[0.5em] y_{+1}&=y+b_1k_1+...+b_sk_s \end{align} He computed order conditions and presented methods up to order $5$ in parametrization and examples. He especially noted the 3/8 method for its symmetry and small error term and the "classical" RK4 method for its simplicity in using always only the last function value in the $y$ updates.
[3] "Beitrag zur näherungsweisen Lösung totaler Differentialgleichungen", Z. f. Math. u. Phys. 46, p. 435-453
On the order dependence of the performance
The Euler method has global error order 1. Which means that to get an error level of $10^{-8}$ (on well-behaved example problems) you will need a step size of $h=10^{-8}$. Over the interval $[0,1]$ this requires $10^8$ steps with $10^8$ function evaluations.
The classical RK4 method has error order 4. To get an error level of $10^{-8}$ you will thus need a step size of $h=0.01$. Over the interval $[0,1]$ this requires $100$ steps with $400$ function evaluations.
If you decrease the step by a factor of $10$ to $h=0.001$, the RK4 method will need $1000$ steps with $4000$ function evaluations to get an error level of $10^{-12}$. This is still much less effort than used in the Euler example above with a much better result.
Using
doubleprecision floating point numbers you will not get a much better result with any method using a fixed step size, as smaller step sizes result in an accumulating floating point noise that dominates the error of the method.