13th August 2024, 15 min read

Recursive Generation of Runge-Kutta Formulas

Original post is here eklausmeier.goip.de/blog/2024/08-13-recursive-generation-of-runge-kutta-formulas.

• 1. Notation
• 2. Example Runge-Kutta methods
• 3. Local discretization error
• 4. Order condition
• 5. Power series expansion for global error
• 6. Cauchy product formula
• 7. Recursive calculation of the order condition

Below text is based on the results in

1. Peter Albrecht: "A New Theoretical Approach to Runge–Kutta Methods"
2. Peter Albrecht: "The Runge–Kutta Theory in a Nutshell"
3. Michael E. Hosea: "A New Recurrence for Computing Runge–Kutta Truncation Error Coefficients"

Bibliography:

1. Peter Albrecht.
2. Michael E. Hosea and LinkedIn.
3. rktec.c computes the coefficients of Albrecht's expansion of the local truncation error of Runge-Kutta formulas.

1. Notation

Consider the ordinary differential value problem with initial condition:

Nomenclature and assumptions:

1. Let $h$ be our step-size, $t_j=t_0+jh$, with $j\in \mathbb{N}$.
2. Let $p\in \mathbb{N}$ be the order of our Runge-Kutta method, see below.
3. The constants $c_i$ are between zero and one, i.e., $0<c_i\ge 1$. They are not necessarily sorted, and they can even repeat.
4. Let $Y_{j+c_i}=y(t_j+c_{i}h)$ be the exact solution of above initial value problem at point $t=t_j+c_{i}h$
5. Let $y_{j+c_i}$ be the approximation according to below Runge-Kutta formula.
6. We assume $y(t)$ is $p+2$ times continously differentiable.

Runge-Kutta methods are written in their nonlinear form

with

$s$ is the number of internal stages of the Runge-Kutta method. $s$ is usually in the range of 2 to 9. The higher $s$ is, the more work you have to do for each step $j$.

Substituting $k_i$ with $y_{j+c_i}$ we obtain

We have

and thus get the $(s+1)$-stage linear representation

In matrix notation this is

Using the matrix

we could write the whole process as ${\tilde{y}}_{j+1}=E{\tilde{y}}_j+\cdots$.

Here we use $c$ as vector and multiindex simultaneously:

and

This corresponds to the classical Runge-Kutta Butcher tableau:

1. Definition. Above method is called an s-stage Runge-Kutta method.

1. It is an explicit method, if $a_{ij}=0$ for $j\ge i$, i.e., $A$ is a lower triangular matrix.
2. It is implicit otherwise.
3. It is called semi-implicit, or SIRK, if $a_{ij}=0$ for $j>i$ but $a_{ii}$ for at least one index $i$.
4. It is called diagonally implicit, or DIRK, if $a_{ii}\ne 0$ and $a_{ij}=0,$ for $j>i$.

In the following we use componentwise multiplications for the vector $c$:

2. Example Runge-Kutta methods

See the book Peter Albrecht, "Die numerische Behandlung gewöhnlicher Differentialgleichungen: Eine Einführung unter besonderer Berücksichtigung zyklischer Verfahren", 1979.

The classical Runge-Kutta method of order 4 with 4 stages.

Kutta's method or 3/8-method of order 4 with 4 stages.

Gill's method of order 4 with 4 stages.

Above examples show that order $p$ could be obtained with $s=p$ stages. Butcher showed that for $p\ge 5$ this is no longer possible.

Butcher method of order 5 with 6 stages.

Runge-Kutta-Fehlberg method of order 4 with 5 internal stages. Also called RKF45. The embedded 5-th order method is only used for step-size control.

Dormand-Prince method of order 4 with internal 5 stages, called DOPRI45.

Implicit Gauß-method of order 4 with 2 internal stages.

Implicit Gauß-method of order 6 with 3 internal stages.

Erwin Fehlberg (1911-1990), John C. Butcher (1933), Carl Runge (1856-1927), Martin Wilhelm Kutta (1867-1944).

3. Local discretization error

The local discretization error is when you insert the exact solution into the numerical formula and look at the error that ensues.

There are two local discretization errors $d_{j+c}\in {\mathbb{R}}^s$ and ${\hat{d}}_{j+1}\in \mathbb{R}$: one for the "internal" stages, and one for the "last" stage.

2. Definition. The $d_{j+c}$ and ${\hat{d}}_j$ are called local discretization errors.

Using

and by Taylor expansion at $t_j$ for $Y_j$ we get for index $i=1$:

The other indexes $i=2,\dots ,s$ are similar. We thus get for all the stages

Using

the "error vectors" ${\gamma }_{\ell }\in {\mathbb{R}}^s$ and the error factor ${\hat{\gamma }}_{\ell }\in \mathbb{R}$ are

The "internal" stages are consistent iff ${\gamma }_{\ell }=0\in {\mathbb{R}}^s$ for $\ell =1,\dots ,p$. Now comes the kicker:

the last stage of the method may furnish approximations of order $p$ even if ${\gamma }_{\ell }=0$ does not hold.

4. Order condition

Define the global error for $y_j$:

and for $f(t,y)$:

2. Theorem. (General order condition.) Assume that the local discretization errors ${\hat{d}}_j$ is $\mathcal{O}\mathcal{(}{\mathcal{h}}^{\mathcal{p}}\mathcal{)}$, $\mathrm{\forall }j$.

(a) The Runge-Kutta method then converges with order $p$, i.e., ${\hat{q}}_j=\mathcal{O}\mathcal{(}{\mathcal{h}}^{\mathcal{p}}\mathcal{)}$, iff

(b) This happens iff for the global error $q_{j+c}$ of the internal stages the following holds:

Proof: Use the fact that

Then

    ☐

Above theorem gives the general order condition for Runge-Kutta methods. However, in this form it is not practical to find the parameters $c,$ $b,$ and $A$.

Further outline:

By one-dimensional Tayler expansion the $Q_{j+c}$ are expressed by $q_{j+c}$.

5. Power series expansion for global error

Let

Further

By Taylor expansion at $t_j$:

3. Theorem. With above definitions the $Q_{j+c}$ can be expressed by powers of $q_{j+c}$:

Proof: The $i$-th component of $Q_{j+c}$ is

Taylor expansion at $y=Y_{j+c_i}$ gives

Hence

    ☐

The following two nonlinear equations vanish:

The right side is analytical at $h=0$, $Q_{t+c}=0$, and $q_{j+c}=0$. The theorem of implicit functions says:

1. the solution $Q_{t+c}$ exists and is unique,
2. the solution is itself analytical.

Its Taylor expansion is

Once can see that the term with $h^0$ and $h^1$ are equal to zero.

The general order condition $b^T{Q}_{t+c}=0$ reduces to below orthogonal relations:

6. Cauchy product formula

In below recap we are interested in the formula, not the convergence conditions, therefore skip these conditions.

4. Theorem. (Cauchy product formula for finitely many power series.) Multiply $n$ power series:

with the multiindex notation:

Similarly, when the power series start at some index ${\beta }_i$.

5. Theorem. (Cauchy product formula for finitely many power series.) Multiply $n$ power series:

with the multiindex notation:

6. Theorem. (Cauchy product formula for infinite many power series.) Multiply infinite many power series:

with the infinite multiindex notation:

For the special case $z_1=z_2=\cdots$

you can substitute the right power of $z$ with $z^{|\alpha |+{\nu }_0}$, therefore

7. Recursive calculation of the order condition

7. Theorem. (Recursion 0.) The $r_i$ and $w_i$ from (T) can be obtained from below recursion formula for $i=2,3,\dots ,p$:

and using the infinite multiindex $\alpha =({\alpha }_1,{\alpha }_2,\dots )$

Proof: The original proof by Albrecht is using induction. This proof is direct and uses the Cauchy product.

We omit $t_j$ in the following. According theorem 3 and (T) we have

Now use (G) and Cauchy product formula:

Now group by common powers of $h$ and you get the formula for $w_i$.     ☐

Recursion 0 is the basis of all further considerations. With the Ansatz

Hosea (1995) in his program rktec.c computes these ${\rho }_{i\ell }$ and ${\sigma }_{i\ell }$.

Main result: The condition $b^T{w}_i(t_j)=0$ must be satisfied for any $f$, and thus for various $e_{i\ell }$. This is the case if

Thus the order condition becomes independent of the $e_{i\ell }$, i.e., independent of the initial value problem at hand.

Recursion 1:

8. Theorem. (Main result) Let the ${\alpha }_{i\ell }\in {\mathbb{R}}^s$ be obtained from Recursion 1. The Runge-Kutta method then converges with order $p\ge 3$ if the last stage has order of consistency of $p$, i.e.,

and if

Hosea (1995) notes that

Through order sixteen, for instance, ... implies a total of 1,296,666 terms (including the quadrature error coefficients) while there are "only" 376,464 distinct order conditions.