33.2 Runge Kutta Integration

The RK integration scheme starts from an initial set of first order ODE dependent variable values ${\bf y}_0$ and advances the ODE solution ${\bf y}$ through well defined steps in the independent variable $x$. Control of the dependent variable errors allows for optimum step size control. We describe here the embedded RK schemes, which use the same intermediate vectors ${\bf k}$ to get an estimate of the error vector ${\bf e}$. An embedded RK step can be represented by the following equations:

$$\frac{d}{dx}{\bf y}(x) = f({\bf y}(x),x)$$ (33.1)

$${\bf k}_i = hf({\bf y}(x)+\sum_{j=1}^sa_{ij}{\bf k}_j,x+c_ih)$$ (33.2)

$${\bf y}(x+h) = {\bf y}(x)+\sum_{i=1}^sb_i{\bf k}_i + O[h^{p+2}]$$ (33.3)

$${\bf y}(x+h) = {\bf y}(x)+\sum_{i=1}^sb_i^*{\bf k}_i + O[h^{p+1}]$$ (33.4)

$${\bf e} = \sum_{i=1}^s(b_i-b_i^*){\bf k}_i + O[h^{p+1}]$$ (33.5)

Here $h$ denotes the step size and $s$ is the size of the embedded RK scheme. Each particular embedded RK scheme is characterized by the four sets of coefficients $a_{ij},b_i,b_i^{*},c_i$, which can be arranged into a table known as a Butcher tableau:

$c_1$	$a_{11}$	$a_{12}$	$\dotsm$	$a_{1s}$
$c_2$	$a_{21}$	$a_{22}$	$\dotsm$	$a_{2s}$
$\vdots$	$\vdots$	$\vdots$	$\ddots$	$\vdots$
$c_s$	$a_{s1}$	$a_{s2}$	$\dotsm$	$a_{ss}$
	$b_1$	$b_2$	$\dotsm$	$b_s$
	$b_1^*$	$b_2^*$	$\dotsm$	$b_s^*$

The coefficients $b$ and $b^{*}$ correspond to RK methods of order $p+1$ and $p$, from which the error vector is constructed via equation 33.5. If only the strict lower triangle of the $a_{ij}$ values are non-zero, then the embedded RK scheme is explicit. Otherwise we have an implicit scheme. Explicit RK schemes are much easier to solve, because all intermediate ${\bf k}$ vectors can be evaluated one after another using previously evaluated ${\bf k}$ vectors. Implicit RK schemes on the other hand have to be solved in an iteratively fashion and are much more expensive. FLASH4 currently has only explicit RK schemes implemented. Error control on each RK step is done by passing a maximum allowed absolute error vector ${\bf e}_{max}$ corresponding to each of the dependent variables in ${\bf y}$. Step size control is currently implemented using information from the obtained error vector ${\bf e}$ after the RK step and the maximum error vector ${\bf e}_{max}$. Using the fact that the error scales as ${\bf O}[h^{p+1}]$ due to the lower $p$ order embedded RK method, the new step size $h_{new}$ is adjusted from the old step size $h_{old}$ as follows

$$h_{new} = h_{old}S\left(max\left\vert\frac{{\bf e}_{max}}{{\bf e}}\right\vert\right)^{1/(p+1)}$$ (33.6)

and used as the next trial step size. A safety factor (typically $S=0.9$) is used to account for the fact that only the leading $h^{p+1}$ term from ${\bf O}[h^{p+1}]$ was used in deriving this formula.

In some situations one needs an initial guess for the very first initial trial step size. This can be done from the supplied ODE vector ${\bf y}$ and the maximum allowed error vector ${\bf e}_{max}$ as follows. Consider the Taylor expansion of ${\bf y}(x+h)$ around $x$:

$${\bf y}(x+h) = {\bf y}(x) + h{\bf y}^{(1)}(x) + h^2{\bf y}^{(2)}(x)/2 + \cdots + h^p{\bf y}^{(p)}(x)/p! + {\bf O}[h^{p+1}],$$ (33.7)

where ${\bf y}^{(p)}$ denotes the $p$-th derivative of ${\bf y}$ with respect to $x$. The expansion error ${\bf O}[h^{p+1}]$ can be directly related to the ${\bf O}[h^{p+1}]$ error of a RK method of order $p$ (cf. equation 33.4). As an approximation one can equate the absolute value of the Taylor remainder term with the error goal of the RK method (i.e. plugging in formulas for both ${\bf O}[h^{p+1}]$ terms):

$$h^{p+1}\vert{\bf y}^{(p+1)}(x^*)\vert/(p+1)! = \vert{\bf e}_{max}\vert,$$ (33.8)

where $x\leq x^*\leq x+h$. Since we are dealing with an approximation, we further set $x^*=x$ (this in effect approximates ${\bf O}[h^{p+1}]$ by its leading term and ignores the $>p+1$ terms), leading to

$$h = min\left((p+1)!\left\vert\frac{{\bf e}_{max}}{{\bf y}^{(p+1)}(x)}\right\vert\right)^{1/p+1}.$$ (33.9)

We then solve for $h$ and take the minimum to satisfy all components. To use this formula we need to evaluate the derivative. The following are second order central difference formulas based on equal spacing $h$ (not to be confused with the RK step size) between functions:

$${\bf y}^{(m+2k+\delta)}(x) = \frac{1}{h^{2k+\delta}}\sum_{n=-(k+\delta)}^{k+\delta} (-1)^{k+n-... ...lta}(2k+\delta)!}{(k+n+\delta)!(k-n+\delta)!}{\bf y}^{(m)}(x+nh) +{\bf R}[h^2],$$ (33.10)

where the remainder term is equal to

$${\bf R}[h^2] = -\frac{k+2\delta}{12}h^2{\bf y}^{(m+2k+2+\delta)}(x^*)$$ (33.11)

with $x-(k+\delta)h\leq x^*\leq x+(k+\delta)h$ and $\delta=0,1$, depending on if one wants even or odd higher order derivatives. For $m = 0$ the derivatives are given in terms of the original functions ${\bf y}(x+ih)$, whereas for $m>0$ we get higher order derivatives from lower order derivatives. These formulas can be derived by setting up a system of derivative linear equations from the Taylor expansions of all the (derivative) functions ${\bf y}^{(m)}(x\pm nh)$ for a certain range of $n$ integers and solving these equations for the lowest derivatives (the technique is discribed in J. Mathews, Computer Derivations of Numerical Differentiation Formulae, International Journal of Mathematics Education in Science and Technology, Volume 34 No 2, pp.280-287). The functions ${\bf y}(x\pm nh)$ could in principle be obtained using the RK equation 33.3, but a much more economical approach would be to use the 1st derivatives (i.e. using $m=1$ in equation 33.10). We have from equation 33.1:

$${\bf y}^{(1)}(x+nh) = f({\bf y}(x+nh),x+nh)$$ (33.12)

$$\approx f({\bf y}(x)+nh{\bf y}^{(1)}(x),x+nh)$$ (33.13)

$$\approx f({\bf y}(x)+nhf({\bf y}(x),x),x+nh),$$ (33.14)

where, since the $nh$ are small, the first order expansion of ${\bf y}(x+nh)$ has been used. Note the nested use of the derivative ODE function $f$ in equation 33.14. Since the goal is to estimate ${\bf y}^{(p+1)}(x)$ via equation 33.10, we set $p=2k+\delta$, from which $k$ and $\delta$ are uniquely determined.

In practice it turns out that equation 33.10 is numerically very unstable for higher order derivatives. We will use it only for estimating low order derivatives. This leaves us with the problem of estimating higher order derivatives in equation 33.9. We follow ideas presented in J. Comp. Appl. Math. 18, 176-192 (1987). To this end we first split the maximum error vector into an error base vector and an error fraction scalar

$${\bf e}_{max} = e_{frac}\cdot {\bf e}_{base}$$ (33.15)

and have the user supplied both quantities. The purpose of this splitting is that we separate the maximum error vector into an accuracy part (the $e_{frac}$, typically of the order of $10^{-6}$ to $10^{-12}$) and a part which relates to the magnitude of the vectors involved in the ODE (the ${\bf e}_{base}$). We further assume that the user has sufficient knowledge about his ODE problem, such that he can construct these two quantities obeying roughly the following inequalities:

$$e_{frac} < 1$$ (33.16)

$${\bf e}_{base} \neq {\bf0}$$ (33.17)

$$\vert{\bf e}_{base}\vert \geq max\;\left\{\vert h^n{\bf y}^{(n)}(x)/n!\vert;n=0,1,\ldots \infty \right\}.$$ (33.18)

With this choice of ${\bf e}_{base}$, all terms in equation 33.7 can be divided by ${\bf e}_{base}$. If we now assume a smooth convergence of the Taylor series and further postulate an exponential type decay between the individual ${\bf e}_{base}$-rescaled Taylor terms (which all must be $\leq 1$ due to the assumptions in equation 33.18)

$$\left\vert\frac{h^n{\bf y}^{(n)}(x)}{n!\;{\bf e}_{base}}\right\vert = \left\vert\frac{h^m{\bf y}^{(m)}(x)}{m!\;{\bf e}_{base}}\right\vert^{n/m}\;\;\;n\geq m,$$ (33.19)

we can replace the higher order derivatives in equation 33.9 by lower order ones and obtain:

$$h = e_{frac}^{1/p+1}\cdot min\left(m!\left\vert\frac{{\bf e}_{base}}{{\bf y}^{(m)}(x)}\right\vert\right)^{1/m}.$$ (33.20)

Note the similar structure to equation 33.9. The higher order contribution has been relegated to $e_{frac}$, which, due to the condition in equation 33.16, leads to larger initial step size estimates for higher order RK schemes. In practice we calculate step size estimates using only $m=2$ and $m=1$. In rare cases when all 1st and 2nd order derivatives of ${\bf y}(x)$ vanish, the (user supplied) default maximum step size estimate is used. This maximum step size estimate can be given for example as a cell/domain size measure, when it is clear that the RK stepper should not step beyond a cell or outside of the domain.

A separate RK stepping routine exists in FLASH4, which performs confined RK steps, which will be explained next. Looking at equation 33.2, we see that for an RK step the ODE function $f$ has to be evaluated repeatedly for several intermediate ${\bf y}_i$ vectors

$${\bf y}_i = {\bf y}(x)+\sum_{j=1}^sa_{ij}{\bf k}_j.$$ (33.21)

For certain applications it may not be possible to evaluate the ODE function at some of these intermediate vectors due to lack of sufficient data. An example would be tracing the (classical) path of a particle through a cell. Only the acceleration vector field of that particular cell is available. If one of the position variables in ${\bf y}_i$ steps outside of the cell boundaries, the acceleration vector at that point cannot be evaluated and the step size $h$ must be reduced. The confined RK stepping routine needs extra boundary vectors corresponding to the variables in the ${\bf y}$ vectors. The boundary vectors can be only a subset of all the ${\bf y}$ vector variables. Those variables which have no corresponding boundary vector entries are assumed to be unbound. Since the boundaries might itself depend on the variables in ${\bf y}$, the boundary vectors are transmitted to the confined RK stepper as functions, which have to be designed by the user.

Currently there are five embedded explicit RK schemes implemented in FLASH, each corresponding to a particular order: 2nd order Euler-Heun, 3rd order Bogacki-Shampine, 4th order Fehlberg, 5th order Fehlberg and 5th order Cash-Karp. Addition of other embedded RK schemes to the code is easy and consists only in extending the existing Butcher tableau repository and adding the appropriate RK scheme keyword to the Runge Kutta unit.