7.5 Runge–Kutta methods

In the method-of-lines approach one ends up effectively integrating a system of ordinary differential equations, which we now generically denote by
-d y(t) = f (t,y). (7.127 ) dt
The majority of approaches in numerical relativity use one-step, multi-stage, explicit Runge–Kutta (RK) methods, which take the following form.

Definition 16. An explicit, s-stage, Runge–Kutta method is of the form

k = f (t ,y ) , 1 n n k2 = f (tn + c2Δt, yn + Δta21k1 ), k = f (t + c Δt, y + Δt (a k + a k )), 3 n 3 n 31 1 32 2 ... ks = f (tn + csΔt, yn + Δt (as1k1 + ...+ as,s−1ks−1)), yn+1 = yn + Δt(b1k1 + ...+ bsks).
with bl(1 ≤ l ≤ s ), cj(2 ≤ j ≤ s), and aji(1 ≤ i ≤ j − 1) real numbers.

Next we present a few examples, in increasing order of accuracy. The simplest one is a forward finite-difference scheme [cf. Eq. (7.70View Equation)].

Example 42. The Euler method (first-order, one-stage):

k1 = f (tn,yn), yn+1 = yn + Δtk1.
It corresponds to all the coefficients set to zero except b1 = 1.

Example 43. Second-order, two-stages RK :

k = f (t ,y ), 1 (n n ) Δt- Δt- k2 = f tn + 2 , 2 k1 , y = y + Δtk . n+1 n 2
The non-vanishing coefficients are c2 = 1∕2 = a21,b2 = 1.

Example 44. Third-order, four-stages RK:

k1 = f (tn,yn ), ( ) k2 = f tn + Δt-,yn + Δt-k1 , 2 2 k3 = f (tn + Δt, yn + Δtk2 ), k4 = f (tn + Δt(, yn + Δtk3 ), ) k1- 2k2- k4- yn+1 = yn + Δt 6 + 3 + 6 .
In the previous case, the number of stages is larger than the order of the scheme. It turns out that it is possible to have a third-order RK scheme with three stages.

Example 45. Third-order, three-stages Heun–RK method

k1 = f ((tn,yn), ) Δt- Δt- k2 = f tn + 3 ,yn + 3 k1 , ( ) k = f t + 2Δt-,y + 2Δt-k , 3 n 3 n 3 2 ( ) yn+1 = yn + Δt k1-+ 3k3- . 4 4
It is also possible to find a fourth-order, four-stages RK scheme. Actually, there are multiple ones, with one and two-dimensional free parameters. A popular choice is

Example 46. The standard fourth-order, four-stages RK method:

k = f (t ,y ) , 1 (n n ) k = f t + Δt-,y + Δt-k , 2 n 2 n 2 1 ( Δt Δt ) k3 = f tn + ---,yn + ---k2 , 2 2 k4 = f (tn + Δt, yn + Δtk3) , ( ) yn+1 = yn + Δt k1+ 2k2-+ 2k3-+ k4- . 6 6 6 6
At this point it is clear that since the RK methods have the same structure, it is not very efficient to explicitly write down all the stages and the final step. Butcher’s tables are a common way to represent them; they have the following structure:


Table 1: The structure of a Butcher table.
0
c2 a21
c 3 a 31 a 32
.. .
c s a s1 a s2 a s,s− 1






b 1 b 2 b s−1 b s

Table 2 shows representations of Examples 42, 43, 44, 45, while Table 3 shows the classical fourth-order Runge–Kutta method of Example 46.


Table 2: From left to right: The Euler method, second and third-order RK, third-order Heun.
0


1
 
0
0 1/2



0 1
 
0
1/2 1/2
1 0 1
1 0 0 1





1/6 2/3 0 1/6
 
0
1/3 1/3
2/3 0 2/3




1/4 0 3/4


Table 3: The standard fourth-order Runge–Kutta method.
0
1/2 1/2
1/2 0 1/2
1 0 0 1





1/6 1/3 1/3 1/6

The above examples explicitly show that up to, and including, fourth-order accuracy there are Runge–Kutta methods of order p and s stages with s = p. It is interesting that even though the first RK methods date back to the end of the 19h century, the question of whether there are higher-order (than four) RK methods remained open until the following result was shown by Butcher in 1963 [93]: s = p cannot be achieved anymore starting with fifth-order accurate schemes, and there are a number of barriers.

Theorem 13. For p ≥ 5 there are no Runge–Kutta methods with s = p stages.

However, there are fifth and sixth-order RK methods with six and seven stages, respectively. Butcher in 1965 [94] and 1985 [95] respectively showed the following barriers.

Theorem 14. For p ≥ 7 there are no Runge–Kutta methods with s = p + 1 stages.

Theorem 15. For p ≥ 8 there are no Runge–Kutta methods with s = p + 2 stages.

Seventh and eighth-order methods with s = 9 and s = 11 stages, respectively, have been constructed, as well as a tenth-order one with s = 17 stages.

7.5.1 Embedded methods

In practice many approaches in numerical relativity use an adaptive timestep method. One way of doing so is to evolve the system of equations two steps with timestep Δt and one with 2Δt. The difference in both solutions at t + 2Δt can be used, along with Richardson extrapolation, to estimate the new timestep needed to achieve any given tolerance error.

In more detail: if we call y2 the solution at t + 2Δt evolved from time t in two steps of size Δt, and &tidle;y1 the solution at the same time advanced from time t in one step of size 2Δt, then the following holds

y2 −-&tidle;y1 ( p+1) y (t + 2Δt ) − y2 = 2p − 1 + 𝒪 (Δt) , (7.128 )
where y denotes the exact solution. Therefore, the term
y2 −-y&tidle;1 (7.129 ) 2p − 1
can be used as an estimate of the error and to choose the next timestep.

Embedded methods also compute two solutions and use their difference to estimate the error and adapt the timestep. However, this is done by reusing the stages. Two Runge–Kutta methods, of order p and ′ p (in most cases – but not always – p′ = p + 1 or p = p′ + 1) are constructed, which share the intermediate function values, so that there is no overhead cost. Therefore, their Butcher table looks as follows:


Table 4: The structure of embedded methods.
0
c2 a21
c3 a31 a32 . ..
c s a s1 a s2 a s,s− 1






b 1 b 2 b s−1 b s






b′ 1 b′ 2 b′ s−1 b′ s

Embedded methods are denoted by ′ p(p ), where p is the order of the scheme, which advances the solution. For example, a 5 (4 ) method would be of fifth order, with a fourth-order scheme, which shares its function calls used to estimate the error.

Table 5 shows the coefficients for a popular embedded method, the seven stages Dormand–Prince 5(4) [146]. Dormand–Prince methods are embedded methods, which minimize a quantification of the truncation error for the highest-order component, which is the one used for the evolution.


Table 5: The 5(4) Dormand–Prince method.
0
1- 5 1- 5
3-- 10 3-- 40 -9- 40
4 5- 44 45- 56 − 15- 32 9--
8- 9 19372- 6561 − 25360- 2187 64448- 6561 − 212- 729
1 9017 ----- 3168 355 − ---- 33 46732 ------ 5247 49 ---- 176 5103 − ------ 18656
1 -35- 384 0 500-- 1113 125- 192 2187- − 6784 11- 84








y1 35 ---- 384 0 500 ----- 1113 125 ---- 192 2187 − ----- 6784 11 --- 84 0








y′1 5179 ------ 57600 0 7571 ------ 16695 393 ---- 640 92097 − ------- 339200 187 ----- 2100 1 --- 40


  Go to previous page Go up Go to next page