19. Runge-Kutta methods

Recall the integral formulation of the ODE on time interval \([t_i, t_{i+1}]\):

\[ y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} f(s, y(s)) ds \]

The goal of Runge-Kutta methods is to obtain higher accuracy by using more accurate numerical integration rule:

\[ y_{i+1} = y_i + \tau \sum_{j=0}^{s-1} b_j f(t_i+c_j \tau, y_i^j) \]

where

• \(s\) is called the number of stages

• \(c_j, b_j\) are points and weights of the quadrature rule

\[ \int_0^1 f(t) dt \approx \sum_{j=0}^{s-1} b_j f(c_j) \]

• \(y_i^j\) are approximations to \(y(t_i + c_j \tau)\), which are also unknown

The exact values at the stage-points are

\[ y(t_i+c_j \tau) = y(t_i) + \int_{t_i}^{t_i+c_j \tau} f(t, y(t)) dt \]

Again, we use numerical integration to define the numerical stage values (skipping the index \(i\)):

\[ y^j = y_i + \tau \sum_{l=0}^{s-1} a_{jl} f(t_i+c_l \tau, y^l) \qquad 0 \leq j < s \]

We use the same integration points as above, however the weigths are adjusted to the smaller interval such that

\[ \int_0^{c_j} f(t) dt \approx \sum_{l=0}^{s-1} a_{jl} f(c_l) \qquad 0 \leq j < s \]

A particular RK method is completely defined by providing the \(a\), \(b\), and \(c\) coefficients, which are usually given in the so called Butcher tableau:

\[\begin{split} \begin{array}{c|ccc} c_0 & a_{0,0} & \dots & a_{0,s-1} \\ \vdots & \vdots & & \vdots \\ c_{s-1} & a_{s-1,0} & \dots & a_{s-1,s-1} \\ \hline & b_0 & \dots & b_{s-1} \end{array} \end{split}\]

Instead of having the \(y^j\) as unknonws, one often solves for \(k\)-values (\(k\) like slope):

\[ k^j = f(y_i + c_j \tau, y_i + \tau \sum_{l=0}^{s-1} a_{jl} k^l) \]

and then

\[ y_{i+1} = y_i + \tau \sum_{l=0}^{s-1} b_l k^l \]

19.1. Explicit RK methods

Explicit RK methods can be seen as extension of the explicit Euler method. The \(a_{jk}\) coefficients form a strictly lower triangular matrix.

Examples are

• the \(RK_2\) - method, aka the explicit mid-point rule

  \[\begin{split} \begin{array}{c|cc} 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 \\ \hline & 0 & 1 \end{array} \end{split}\]

  It is evaluated as

  \[\begin{split} \begin{array}{rcl} y^0 & = & y_i \\ y^1 & = & y_i + \tfrac{\tau}{2} f(t_i, y^0) \end{array} \end{split}\]

  and then

  \[ y_{i+1} = y_i + \tau f(t_i+\tfrac{h}{2}, y^1) \]

  Evaluation in \(k\)-values:

  \[\begin{split} \begin{array}{rcl} k^0 & = & f(t_i, y_i) \\ k^1 & = & f(t_i+\frac{\tau}{2}, y_i + \frac{\tau}{2} k^0) \end{array} \end{split}\]

  and then

  \[ y_{i+1} = y_i + \tau k^1 \]

• the \(RK_4\) - method, aka THE classical Runge Kutta method

\[\begin{split} \begin{array}{c|cccc} 0 & \\ \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & 0 & \frac{1}{2} \\ 1 & 0 & 0 & 1 \\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array} \end{split}\]

\[\begin{split} \begin{array}{rcl} y^0 & = & y_i \\ y^1 & = & y_i + \tfrac{\tau}{2} f(t_i, y^0) \\ y^2 & = & y_i + \tfrac{\tau}{2} f(t_i+\tfrac{\tau}{2}, y^1) \\ y^3 & = & y_i + \tau f(t_i+\tfrac{\tau}{2}, y^2) \\ \end{array} \end{split}\]

and then

\[ y_{i+1} = y_i + \frac{\tau}{6} \big(f(t_i, y^0) + 2 f(t_i+\tfrac{h}{\tau}, y^1) + 2 f(t_i+\tfrac{\tau}{2}, y^2) + f(t_i+\tau, y^3) \big) \]

\[\begin{split} \begin{array}{rcl} k^0 & = & f(t_i, y_i) \\ k^1 & = & f(t_i+\frac{\tau}{2}, y_i + \frac{\tau}{2} k^0) \\ k^2 & = & f(t_i+\frac{\tau}{2}, y_i + \frac{\tau}{2} k^1) \\ k^3 & = & f(t_i + \tau, y_i + \tau k^2) \end{array} \end{split}\]

and then

\[ y_{i+1} = y_i + \frac{\tau}{6} \big( k^0 + 2 k^1 + 2 k^2 + k^3 \big) \]

19.2. Implicit Runge Kutta methods

Methods of optimal accuracy are obtained based on Gaussian quadrature:

Gauss-Legendre 2 method:

\[\begin{split} \begin{array}{c|cccc} \frac{1}{2} - \frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4} - \frac{\sqrt{3}}{6} \\ \frac{1}{2} + \frac{\sqrt{3}}{6} & \frac{1}{4} + \frac{\sqrt{3}}{6} & \frac{1}{4} \\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \end{split}\]

If not all time-scales are resolved, then \(c_s = 1\) leads to better stability (i.e. L-stability). The other points are chosen for optimal accuracy. These methods are called Radau IIA.

Two-point Radau IIA method:

\[\begin{split} \begin{array}{c|cc} \frac{1}{3} & \frac{5}{12} & \frac{-1}{12} \\ 1 & \frac{3}{4} & \frac{1}{4} \\ \hline & \frac{3}{4} & \frac{1}{4} \\ \end{array} \end{split}\]

Many examples can be found here: Butcher tableaus examples.

19.3. Implementation of a Runge-Kutta time-stepper for arbitrary Butcher tableaus.

They are implemented in ASC-ODE, branch RungeKutta. Merge that branch into your main branch:

git remote add upstream https://github.com/TUWien-ASC/ASC-ODE.git
git fetch upstream
git merge upstream/RungeKutta

Solve a non-linear system of equations for \(k^0, \ldots k^{s-1}\) such that

\[\begin{split} \left( \begin{array}{c} k^0 \\ \vdots \\ k^{s-1} \end{array} \right) - \left( \begin{array}{c} f( y_i + \tau (a_{0,0} k^0 + \cdots + a_{0,s-1} k^{s-1}) ) \\ \vdots \\ f( y_i + \tau (a_{s-1,0} k^0 + \cdots + a_{s-1,s-1} k^{s-1}) ) \end{array} \right) = 0 \end{split}\]

This system of equations is rewritten in compact form for the unknown \(k \in {\mathbb R}^{sn}\):

\[ k - \widetilde f ( \tilde y_i + \tau A \otimes k) = 0 \]

Where

• \(k = (k^0, \ldots k^{s-1})\), with \(k^j \in {\mathbb R}^n\)

• \(\tilde y_i = (\underbrace{y_i, \ldots , y_i}_{s-{\text times}}) \in {\mathbb R}^{sn}\)

• \(\tilde f : {\mathbb R}^{sn} \rightarrow {\mathbb R}^{sn} : (x_0, \ldots, x_{s_-1}) \mapsto (f(x_0), \ldots, f(x_{s-1}))\)

• for \(A \in {\mathbb R}^{h \times w}\) define \(A \otimes : {\mathbb R}^{wn} \rightarrow {\mathbb R}^{hn} : (x_0, \ldots, x_{w-1}) \mapsto (y_0, \ldots , y_{h-1})\) such that \(y_j = \sum_l a_{jl} x_l\).

We derive two classes from NonlinearFunction for that (named MultipleFunction and MatVecFunc).

Setup the Runge-Kutta time-stepper, and solve the ODE.

19.3.1. Runge-Kutta methods of arbitrary order

If we have fixed the integration points \(c_j\), we can compute the \(a_{jl}\) and \(b_j\) coefficients such that the integration rules are (at least) exact for polynomials up to order \(s-1\):

\[ \sum_{l=0}^{s-1} b_{l} {c_l}^k = \int_0^1 x^k \, dx \qquad 0 \leq k < s \]

and

\[ \sum_{l=0}^{s-1} a_{jl} {c_l}^k = \int_0^{c_j} x^k \, dx \qquad 0 \leq j,k < s \]

Setup and solve linear equations for computing the \(a\) and \(b\) coefficients.

Use the \(c\)-coefficients from Gauss-Legendre 2, and two-point Radau IIA, and check whether you recover the same \(a\) and \(b\) coefficients.

Try \(c\)-coefficients from Gauss-Legendre 3, and three-point Radau II A (from the web-page above).

Then, find \(c\)-coefficients by solving for the roots of Legendre-polynomials, and transform to \([0,1]\). A ready to use function can be taken from Numerical recipes, section 4.6

Similar for the Radau IIA: Find roots of the \(P^{(1,0)}\) Jacobi-polynomial, transform to \([0,1]\), and add \(c_j=1\). This is done by function gaujac from Numerical recipes. The logarithmic-gamma function is lgamma in C++.

19.4. Exercises

• Try out and compare various Runge-Kutta methods already implemented for the examples above.

• Implement and test an ExplicitRungeKutta time-stepper for arbitrary Butcher tableaus (with a lower triangular matrix for the \(a\)-coefficents).