The Runge-Kutta Family (RK)

Now we revert to a ``Taylor series'' interpretation of problem (see 0.1.3). We consider in some detail the EXPLICIT CASE and will make only cursory comments on the IMPLICIT Runge-Kutta methods. RK are single-step methods and can be either explicit or implicit. Later we'll consider multi-step methods.

Perhaps the most popular ODE integrator around, because it is explicit, easily programmable, of high order, and applicable to even mildly stiff problems: $4^{th}$-order ERK (Explicit Runge Kutta) or ERK4. The $2^{nd}$ order $\rightarrow ERK2$ is known as Heun's Method and is also popular but used less often.

$$\begin{array}{ll} \mbox {\underline {Pro's:}} & \mbox {simple... ...ion (invariably one uses small-order ERK method).} \end{array}$$

All ERK's are written as

$$y_{n+1}=y_n+h\sum^\nu_{j=1}b_j f\left(x_n+c_jh, y(x_n +c_j h))\equiv y_n+h F(x_n, y_n, h, f)\right) n\ge 0,$$

where $c_j's$ are between 0 and 1. The whole point is to specify $y$ at the locations $x_n +c_1h, x_n+c_2h, \cdots, x_n+c_\nu h$ and find the corresponding $b_j$. These $b_j's$ must sum to 1 so that we get a weighted average. What's the criteria? The choice of the entries in the vectors ${\bf b}$ and ${\bf c}$ make the leading terms in the truncation error equal to zero. Additionally, we want

$$F(x, Y(x), h; f)\approx Y'(x)=f(x, Y(x)),$$$$\mbox {for small } h$$

Example: Take trapezoidal with an Euler predictor step

$$y_{n+1}=y_n+\frac{h}{2}\Big[f(x_n, y_n)+f(x_{n+1}, y_n+hf(x_n, y_n))\Big]$$ (34)

or $y_{n+1}=y_{n}+hF$    What's $F$? See Figure 6.

Figure 6: Geometrical interpretation of average slope

$\therefore F$ is average slope on $[x, x+h]!$

Example) Another method based on average slope is

$$y_{n+1}=y_n+hf\left(x_n+\frac12 h, y_n+\frac12hf(x_n, y_n)\right)$$ (35)

$$F=f\left(x+\frac12 h, y_n+\frac12 hf\right)$$

both (34) and (35) are $2^{nd}$ order.

To illustrate general procedure: ERK2 ``Heun's Method''

$$Y(x+h)=Y(x)+hY'(x)+\frac{h^2}{2!}Y''(x)+\frac{h^3}{3!}Y'''(x)\cdots$$ (36)

and

$$\left\{ \begin{array}{l} Y'=f\\ Y''=f_x+f_Yy'=f_x+f_Yf\\ Y'... ...x+f_Yf)f_Y+f(f_{xY}+f_{YY}f)\\ \mbox {etc.} \end{array}\right.$$

so (36) can be written as

$$Y(x+h) = Y+hf+\frac12 h^2(f_x+f_Yf)+{\cal O }(h^3)$$ (37)

$$= Y+\frac{h}{2} f+\frac12h[f+hf_x+hff_Y]+{\cal O}(h^3)$$ (38)

but note that $f(x+h, Y+hf)=f+h f_x+hff_Y+{\cal O}(h^2)$

so substitute in (36)

$$Y(x+h)=Y+\frac12 hf+\frac{h}{2}f(x+h, Y+hf)+{\cal O}(h^3)$$

$$\left[\begin{array}{c} \therefore y(... ... ERK2 method}\\ \mbox{known as \underline {Heun's Method}}. \end{array}\right.$$

ERK4 Classical

$$y(x+h)=y(x)+\frac{1}{6}(F_1+2F_2+2F_3+F_4)$$

$$\left\{\begin{array}{l} F_1=hf(x,y)\\ F_2=hf(x+\frac12 h, y+\dis... ...12 h, y+\displaystyle \frac12 F_2)\\ \\ F_4=hf(x+h, y+F_3) \end{array}\right.$$

General Method:

$$y_{n+1}=y_n+h\sum^\nu_{j=1}b_jf(x_n+c_jh, y(x_n+c_jh))\quad n=0, 1\ldots$$

let ``approximant''

$$y(x_n+c_jh)=\xi_j\quad j=1, 2, \cdots \nu$$

let

$$c_1=0\quad {\mbox so }\quad\xi_1=y_n$$

The idea is to express each $\xi_j$ with $j=2, 3, \cdots \nu$ by updating $y_n$ with a linear combination of

$f(x_n,\xi), f(x_n, +hc_2, \xi_2)\cdots f(x_n+h c_{j-1}, \xi_{j-1}),$i.e.

$$\xi_1 = y_n$$

$$\xi_2 = y_n+ha_{2,1}f(x_n,\xi_1)$$

$$\xi_3 = y_n+ha_{3,1}f(x_n, \xi_1)+ha_{3,2}f(x_n+c_2h, \xi_2)$$

$$\vdots\quad\vdots\quad \nu-1$$

$$\xi_{\nu} = y_n+h\sum_{i=1}a_{\nu,i}f(x_i+c_ih,\xi_i)$$

$$\Rightarrow y_{n+1} = y_n+h\sum^{\nu}_{i=1}b_jf(x_n+c_jh, \xi_j)$$

The matrix $A \equiv A_{j,i}$     $j,i=1,2\cdots \nu\quad\equiv RK$ matrix

$${\bf b} = [b_1 b_2\cdots b_{\nu}]^T\equiv RK$$

$${\bf c} = [c_1 c_2\cdots c_{\nu}]^T\equiv RK$$

and we say that method has ``$\nu$stages''

Take simple case, with $\nu=2$. Assume $f$ smooth, for simplicity,

$$\xi_1=y_n$$

$$f(x_n+c_2h,\xi_2)=f(x_n+c_2h,y_n+a_{21}, hf(x_n, y_n))$$

$$=f(x_n,y_n)+h\left[c_2\frac{\partial f}{\partial x}(x_n,y_n)+a_{2,1}\frac{\partial f}{\partial y} (x_n, y_n)f(x_n, y_n)\right]+{\cal O}(h^2)$$

$$\therefore y_{n+1}=y_n+h(b_1+b_2)f(x... ...rtial f} {\partial x}+a_{21}\frac{\partial f}{\partial y}f\right]+{\cal O}(h^3)$$ (39)

but we note that $Y''=\displaystyle \frac{\partial f}{\partial x}+\frac{\partial f}{\partial Y}f$    from $Y'=f(x,y),$ the IVP, and exact solution

$$Y_{n+1}=Y_n+hf(Y_n,x_n)+\frac12h^2\Big[\frac{\partial f}{\partial x}(x_n,Y_n)+ \frac{\partial f}{\partial Y}(x_n,Y_n)f\Big]+ \mathcal{O} (h^3)$$ (40)

$\therefore$ comparing (39) and (40) gives us

$$b_1+b_2=1\quad b_2c_2=\frac12\quad a_{2,1}=c_2$$ (41)

So we see that a 2-stage is not uniquely defined. Popular choices are

0
$\frac12$	$\frac12$
	0	$1$

        

0
$\frac23$	$\frac23$
	$\frac14$	$\frac34$

        

0
$1$	$1$
	$\frac12$	$\frac12$

``RK Tableaux''

${\bf c}$	$A$
	${\bf b}^T$

3-Stage Examples (both are 3$^{rd}$ order)

``Classical''

0
$\frac12$	$\frac12$
$1$	$-1$	$2$
	$\frac16$	$\frac23$	$\frac16$

         System

0
$\frac23$	$\frac23$
$\frac23$	0	$\frac23$
	$\frac14$	$\frac38$	$\frac38$

4-Stage (fourth-order)

0
$\frac12$	$\frac12$
$\frac12$	0	$\frac12$
$1$	0	0	$1$
	$\frac16$	$\frac13$	$\frac13$	$\frac16$

Compare to $Y_{n+1}=y_n+\displaystyle \frac16(F_1+2F_2+2F_3+F_4)$

$$\left\{\begin{array}{l} F_1=hf(x_n,y_n)\\ F_2=hf\left(x_n+\frac1... ..., y_n + \displaystyle \frac12 F_2)\\ F_4=hf(x_n+h, y_n+F_3) \end{array}\right.$$

Relation between order and number of stages:

Max stages	1	2	3	4	5	6	7	8
Max order	1	2	3	4	4	5	6	6

Not worth it usually if number of stages is much greater than max order.

IMPLICIT RK(IRK): We won't study in detail. Complicated. Used in stiff equation solutions because they exhibit superior stability properties.

Generally: $\xi_j=y_n+h\displaystyle \sum^{\nu}_{i=1}a_{j,i}f(x_n+c_ih,\xi_i)\quad j=1,2, \cdots \nu$

$$y_{n+1}=y_{n}+h\sum^{\nu}_{j=1}b_jf(x_1 +c_jh, \xi_j)$$

Here $A\equiv(a_{j,i})$ is no longer lower triangular.

$$\sum^{\nu}_{i=1}a_{j,i}=c_j\quad j=1,2,\cdots \nu$$   by convention$$.$$

So for $y\in \mathbb{R}^d$ we get $\nu$ coupled algebraic equations.

For references see J.C. Butcher ``The Numerical Anaylsis of ODE'S,'' John Wiley Publ. (He's one of the world's experts).

$\Box$

RK-FEHLBERG (RKF) This is an adaptive step size method which illustrates how error control can be incorporated into RK. The technique is not only applicable to RK, though. It is based on an idea similar to Richardson extrapolation The idea is to adapt the step size to control the error and ensure error is kept within a (reasonable) specified bound, $epsilon$. It is an example of a scheme that incorporates ``error control''.

Most popular $\left\{\begin{array}{ll} 5 \mbox { stage} & 4^{th} \mbox{ order}\\ 6 \mbox { stage} & 5^{th} \mbox { order} \end{array}\right.$

General strategy:

Take $\left\{\begin{array}{ll} Y'=f(x,Y)\\ Y(x_0)=Y_0 \end{array}\right.$

For presentation purposes we will assume that the schemes under consideration are explicit, single step and that the scheme for $w_{n+1}$ has a truncation error $\tau_{n+1}(h)=O(h^m)$ Assume the scheme is of the form $\left\{\begin{array}{ll} w_{n+1}=w_n+hf(x_n, w_n, h)\\ w_0=Y_0 \end{array}\right.$

Use another scheme with $\tilde{\tau}_{n+1}(h)={\cal O}(h^{m+1})$

Take $\left\{\begin{array}{ll} z_{n+1}=z_n+hf(x_n, z_n, h)\\ z_0=Y_0 \end{array}\right.$

Assume that $w_n\approx z_n\approx Y(x_n)$, i.e. assume that the schemes $w$ and $z$ are convergent. The if we subtract

$$Y_{n+1}-w_{n+1} = Y_{n+1}-w_n-hf(x_n, w_n,h)\approx Y_{n+1}-Y_n-hf(x_n, y_n, h)$$

$$= h\tau_{n+1}(h)$$

so     $\tau_{n+1}\approx\displaystyle \frac{1}{h}\bigg[Y_{n+1}-w_{n+1}\bigg] =\frac{1}{h} [Y_{n+1}-z_{n+1}]+\frac{1}{h}[z_{n+1}-w_{n+1}]$

$\therefore \tau_{n+1}\approx \tilde{\tau}_{n+1}\quad +\frac{1}{h}[z_{n+1}-w_{n+1}]$

but $\tau_{n+1}=O(h^m)$ and $\tau_{n+1}=O(h^{m+1})$ hence the major error contribution

$$\tau_{n+1}\approx \frac{1}{h}[z_{n+1}-w_{n+1}]$$

The idea is to adjust the step size to stay within a certain error bound.

Since $\tau_{n+1}(h)={\cal O}(h^m) \therefore \tau_{n+1}(h)=Kh^m$

$$\tau_{n+1}(qh)\approx K(qh)^m=q^m(Kh^m)\approx q^m\tau_{n+1}(h)= \frac{q^m}{h}(z_{n+1}-w_{n+1})$$

So choose $q$ such that

$$\frac{q^m}{h}\vert z_{n+1}-w_{n+1}\vert\approx \vert\tau_{n+1}(qh)\vert\le \varepsilon$$

or

$$q\le\left(\frac{\varepsilon h}{\vert z_{n+1}-w_{n+1}\vert}\right)^{\frac{1}{m}}$$

$\Box$ Hence, $q$ is a multiplicative factor that scales the time step $h$, $epsilon$ is the error tolerance and is a specified input to the code.

You might be wondering about the fact that in the above expression $z_{n+1}$ and $w_{n+1}$ appear and these are quantities that are sought after. What's done algorithmically is not unique: one possibility is to take a step and produce some proxi $z_{n+1}$ and $w_{n+1}$. Make the estimate as per above equation, and if it is not satisfied, reduce $q$ till it is. Then do the real time step.

Note, also, that on implementation, it is possible for the above condition not to be satisfied (either because you chose an $epsilon$ that is unreasonable small, or because the method should not be applied to the IVP in the first place). Hence, an escape sequence should be supplied in the algorithm so that the user gets a warning of the method's failure.

An implementation of this is ERKF4: use a 5$^{th}$ order RK to estimate the 4$^{th}$ order RK (there's one for 5$^{th}$-order equations as well and it is easily derivable.)

ALGORITHM

INPUT         $x_0,b$, initial condition $\alpha$, TOL, $h_{max}$, $h_{min}$
OUTPUT     $x, y,h$

Step 1	$x=x_0$     ;     $y=\alpha, h=h_{max}$, FLAG=1; output($x,y$);
Step 2	While (FLAG=1) do Step 3-11
Step 3

$$K_1 = hf(x,y)$$

$$K_2 = hf(x+\frac14h, y+\frac14K_1)$$

$$K_3 = hf\left(x+\frac{3}{8}h, y+\frac{3}{32}K_1+\frac{9}{32}K_2\right)$$

$$K_4 = hf\left(x+\frac{12}{13}h, y+\frac{1932}{2197}K_1-\frac{7200}{2197}K_2+\frac{7296}{2197}K_3\right)$$

$$K_5 = hf\left(x+h, y+\frac{439}{216}K_1-8K_2+\frac{3680}{513}K_3 -\frac{845}{4104}K_4\right)$$

$$K_6 = hf\left(x+\frac{h}{2}, y-\frac{8}{27}K_1+2K_2-\frac{3544}{2565}K_3\right.$$

$$+ \left.\frac{1859}{4104}K_4-\frac{11}{40}K_5\right)$$

$$R=\frac{1}{h}\Big\vert\frac{1}{360}K_1-\frac{128}{4275}K_3-\frac{2197}{75240}K_4+\frac{1}{50} K_5+\frac{2}{55}K_6\Big\vert$$

$\%$ Note: $R=\Big\vert\tilde{y}_{n+1}-y_{n+1}\Big\vert\displaystyle \frac{h}{2}$

Step 5	$\delta=0.84$ (TOL/R) $^{\frac14}$
Step 6	if $R\le$ TOL do steps $7\& 8$
Step 7	$x=x+h$
	$y=y+\displaystyle \frac{25}{216}K_1+\frac{1408}{2565}K_3+ \frac{2197}{4104}K_4-\frac15 K_5$
Step 8	output $(x,y, \delta)$
Step 9	if $\delta\le 0.1$ then $h=0.1h$
	else if $\delta\ge 4$ then set $h=4h$
	else set $h=\delta h \%$ (Calc new $h$).
Step 10	if $h>h_{max}$ then $h=h_{max}$
Step 11	If $x\ge b$ the FLAG $=0$
	else if $x+h>b$ then $h=b-x$
	else if $h<h_{min}$ then
	set FLAG $=0$
	output ('minimum $x$ step is exceeded').
Step 12	STOP, END

$\Box$

Two final remarks are in order:

1) the formal introduction of the method in a schematic way is meant to convey that the trick is general and can be constructed using other schemes, in addition to RK. The overall choice of methods to combine should be dictated by the overall goal: to exploit adaptivity to make the computation more efficient, and to make codes that are more robust and less prone to users' bad choice of time stepping parameters.

2) All too often scientists will try to use adaptivity in step size to try to circumvent a problematic nature of an IVP. Adaptivity is useful when the problem exhibits solutions that have periods of high activity, followed by periods of low activity. It is important that you learn to recognize the difference between an IVP that is STIFF and one that merely has periods of heightened activity. It is common for people to think that one can still use an explicit method with adaptivity to counteract the stiffness of an IVP. This is not true in general, although it might work in some circumstances. Adaptivity might force a method to stay within its stability region, by making $h$ small enough. It works properly if the stability range gets significantly smaller and larger as the code steps through the approximate solution. It does not work, when the code is forced to use a small step and then is forced to keep such small step for the rest of the integration interval: you're better off coding something up that is simpler and more robust and forego adaptivity, since it brings no benefit. And of course, it goes without saying, that if you choose a method for which no $h$ value could produce a stable approximation, then adaptivity is not going to help things...