, 15 min read

Recursive Generation of Runge-Kutta Formulas

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:

$$ \dot y(t) = f(t,y(t)), \qquad y(t_0) = y_0 \in\mathbb{R}. $$

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\ge1$. They are not necessarily sorted, and they can even repeat.
  4. Let $Y_{j+c_i} = y(t_j+c_ih)$ 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

$$ k_i := f\left(t_j+c_ih, y_j + h \sum_{\ell=1}^s a_{i\ell}k_\ell\right), \qquad i=1,\ldots,s; $$

with

$$ y_{j+1} = y_j + h \sum_{\ell=1}^s b_i k_\ell, \qquad j=0,1,\ldots,m-1. $$

$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

$$ y_{j+c_i} := y_j + h \sum_{\ell=1}^s a_{i\ell} k_\ell. $$

We have

$$ k_i = f(t_j+c_ih,y_{j+c_i}), $$

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

$$ \eqalign{ y_{j+c_i} &= y_j + h \sum_{\ell=1}^s a_{i\ell} f(t_j+c_ih,y_{j+c_i}), \qquad i=1,\ldots,s,\cr y_{j+1} &= y_j + h \sum_{\ell=1}^s b_\ell f(t_j+c_\ell h,y_{j+c_i}), \qquad j=0,\ldots,m-1. } $$

In matrix notation this is

$$ \eqalign{ y_{j+c} &= y_j e + h A f(t_j+c_ih,y_{j+c}), \qquad\hbox{"internal" stages}\cr y_{j+1} &= y_j + hb f(t_{j+c},y_{j+c}), \qquad\hbox{"last" stage.} } $$

Using the matrix

$$ E = \pmatrix{ 0 & \ldots & 0 & 1\cr \vdots & \ddots & \vdots & \vdots\cr 0 & \ldots & 0 & 1 } \in \mathbb{R}^{s\times(s+1)}. $$

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:

$$ y_{j+c} := \pmatrix{ y_{j+c_1}\cr \vdots\cr y_{j+c_s} }, \quad e := \pmatrix{ 1\cr \vdots\cr 1 }, \quad f(t_{j+c},y_{j+c}) := \pmatrix{ f(t_tj+c_1h,y_{j+c_1})\cr \vdots\cr f(t_j+c_sh,y_{j+c_s}) }, $$

and

$$ A = \pmatrix{ a_{11} & \ldots & a_{1s}\cr \vdots & \ddots & \vdots\cr a_{s1} & \ldots & a_{ss} }, \qquad b := \pmatrix{ b_1\cr \vdots\cr b_s }, \qquad c = \pmatrix{ c_1\cr \vdots\cr c_s } $$

This corresponds to the classical Runge-Kutta Butcher tableau:

$$ \begin{array}{c|c} c & A\\ \hline & b^T \end{array} $$

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\gt i$ but $a_{ii}$ for at least one index $i$.
  4. It is called diagonally implicit, or DIRK, if $a_{ii}\ne0$ and $a_{ij}=0,$ for $j\gt i$.
# Runge-Kutta methods ## explicit ## implicit ### semi-implicit ### diagonally implicit

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

$$ c^2 = \pmatrix{c_1^2\cr \vdots\cr c_s^2}, \qquad\ldots\qquad, c^\ell = \pmatrix{c_1^\ell\cr \vdots\cr c_s^\ell}. $$

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.

$$ \begin{array}{c|cccc} {1\over2} & {1\over2}\\ {1\over2} & 0 & {1\over2}\\ 1 & 0 & 0 & 1\\ \hline & {1\over6} & {1\over3} & {1\over3} & {1\over6} \end{array} $$

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

$$ \begin{array}{c|cccc} {1\over3} & {1\over3}\\ {2\over3} & -{1\over3} & 1\\ 1 & 1 & -1 & 1\\ \hline & {1\over8} & {3\over8} & {3\over8} & {1\over8} \end{array} $$

Gill's method of order 4 with 4 stages.

$$ \begin{array}{c|cccc} {1\over2} & {1\over2}\\ {1\over2} & {\sqrt{2}-1\over2} & {2-\sqrt{2}\over2}\\ 1 & 1 & -{\sqrt{2}\over2} & 1+{\sqrt{2}\over2}\\ \hline & {1\over6} & {2-\sqrt{2}\over6} & {2+\sqrt{2}\over6} & {1\over6} \end{array} $$

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

Butcher method of order 5 with 6 stages.

$$ \begin{array}{c|ccccc} {1\over4} & {1\over4}\\ {1\over4} & {1\over8} & {1\over8}\\ {1\over2} & 0 & -{1\over2} & 1\\ {3\over4} & {3\over16} & 0 & 0 & {9\over16}\\ 1 & -{3\over7} & {2\over7} & {12\over7} & -{12\over7} & {8\over7}\\ \hline & {7\over90} & 0 & {32\over90} & {12\over90} & {32\over90} & {7\over90} \end{array} $$

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.

$$ \begin{array}{c|cccccc} {1\over4} & {1\over4}\\ {3\over8} & {3\over32} & {9\over32}\\ {12\over13} & {1932\over2197} & -{7200\over2197} & {7296\over2197}\\ 1 & {439\over216} & -8 & {3680\over513} & -{845\over4104}\\ {1\over2} & -{8\over27} & 2 & -{3544\over2565} & {1859\over4104} & -{11\over40}\\ \hline p=5 & {16\over135} & 0 & {6656\over12825} & {28561\over56430} & -{9\over50} & {2\over55}\\ p=4 & {25/216} & 0 & {1408\over2465} & {2197\over4104} & -{1\over5} & 0 \end{array} $$

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

$$ \begin{array}{c|ccccccc} {1\over5} & {1\over5}\\ {3\over10} & {3\over40} & {9\over40}\\ {4\over5} & {44\over45} & -{56\over15} & {32\over9}\\ {8\over9} & {19372\over6561} & -{25360\over2187} & {64448\over6561} & -{212\over729}\\ 1 & {9017\over3168} & -{355\over33} & {46732\over5247} & {49\over176} & -{5103\over18656}\\ 1 & {35\over384} & 0 & {500\over1113} & {125\over192} & -{2187\over6784} & {11\over84}\\ \hline p=5 & {35\over384} & 0 & {500\over1113} & {125\over192} & -{2187\over6784} & {11\over84}\\ p=4 & {5179\over57600} & 0 & {7571\over16695} & {393\over640} & -{92097\over339200} & {187\over 2100} & {1\over40} \end{array} $$

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

$$ \begin{array}{c|cc} {3-\sqrt{3}\over6} & {1\over4} & {3-2\sqrt{3}\over12}\\ {3+\sqrt{3}\over6} & {3+2\sqrt{3}\over12} & {1\over4}\\ \hline & {1\over2} & {1\over2} \end{array} $$

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

$$ \begin{array}{c|ccc} {5-\sqrt{15}\over10} & {5\over36} & {10-3\sqrt{15}\over45} & {25-6\sqrt{15}\over180}\\ {1\over2} & {10+3\sqrt{15}\over72} & {2\over9} & {10-3\sqrt{15}\over72}\\ {5+\sqrt{15}\over10} & {25+6\sqrt{15}\over180} & {10+3\sqrt{15}\over45} & {5\over36}\\ \hline & {5\over18} & {4\over9} & {5\over18} \end{array} $$

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 $d_{j+1}\in\mathbb{R}$: one for the "internal" stages, and one for the "last" stage.

$$ \eqalign{ Y_{j+c} &= Y_j e + h A f(t_{j+c},Y_{j+c}) + h d_{j+c},\cr Y_{j+1} &= Y_j + h b^T f(t_{j+c},Y_{j+c}) + h \hat d_{j+1}, \qquad j=0,\ldots,m-1. }\tag{*} $$

2. Definition. The $d_{j+c}$ and $d_j$ are called local discretization errors.

Using

$$ Y_j^{(i)} = {d^i y(t_j)\over dt^i}, \qquad Y_j^{(1)} = f(t_j,y(t_j)), $$

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

$$ \eqalign{ Y_{j+c_1} &= y(t_j+c_1h) = Y_j + {1\over1!}Y_j^{(1)} c_1h + {1\over2!}Y_j^{(2)}(c_1h)^2 + \cdots\cr f(t_{j+c_1},y(t_{j+c_1})) &= \dot y(t_j+c_1h) = Y_j^{(1)} + {1\over1!}Y_j^{(2)}c_1h + {1\over2!}Y_j^{(3)}(c_1h)^2 + \cdots } $$

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

$$ \eqalign{ d_{j+c} &= \pmatrix{ \gamma_{11} Y_j^{(1)} h + \gamma_{12} Y_j^{(2)} h^2 + \gamma_{13} Y_j^{(3)} h^3 + \cdots\cr \qquad\vdots \qquad\qquad\qquad \ddots\cr \gamma_{s1} Y_j^{(1)} h + \gamma_{s2} Y_j^{(2)} h^2 + \gamma_{s3} Y_j^{(3)} h^3 + \cdots } = \sum_{\ell=1}^{p+1} \gamma_\ell Y_j^{(\ell)} h^{\ell-1} + \cal{O}(h^{p+1}),\cr \hat d_{j+1} &= \sum_{\ell=1}^{p+1} \hat\gamma_\ell Y_j^{(\ell)} h^{\ell-1} + \cal{O}(h^{p+1}). } $$

Using

$$ \Gamma = \left( \gamma_1, \ldots, \gamma_{p+1} \right) = \pmatrix{ \gamma_{11} & \ldots & \gamma_{1,p+1}\cr \vdots & \ddots & \vdots\cr \gamma_{s1} & \ldots & \gamma_{s,p+1} } \in \mathbb{R}^{s\times(p+1)}, $$

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

$$ \eqalign{ \gamma_\ell &= {1\over\ell!} c^\ell - {1\over(\ell-1)!} A c^{\ell-1}, \qquad \ell=1,\ldots,p+1;\cr \hat\gamma_\ell &= {1\over\ell!} - {1\over(\ell-1)!} b^T c^{\ell-1}, \qquad c^\ell := \left(c_1^\ell,\cdots,c_s^\ell\right)^T } $$

The "internal" stages are consistent iff $\gamma_\ell=0\in\mathbb{R}^s$ for $\ell=1,\ldots,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$:

$$ q_{j+c} := Y_{j+c} - y_{j+c} \in\mathbb{R}^s, \qquad \hat q_{j+1} := Y_{j+1} - y_{j+1} \in\mathbb{R} $$

and for $f(t,y)$:

$$ Q_{j+c} := f(t_{j+c},Y_{t+c}) - f(t_{j+c},y_{j+c}) \in\mathbb{R}^s. $$

2. Theorem. (General order condition.) Assume that the local discretization errors $\hat d_j$ is $\cal O(h^p)$, $\forall j$.

(a) The Runge-Kutta method then converges with order $p$, i.e., $\hat q_j = \cal O(h^p)$, iff

$$ b^T Q_{t+c} = {\cal O}(h^p), \quad\forall j. $$

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

$$ \eqalign{ q_{j+c} &= {\cal O}(h^p) + h d_{j+c} + h A Q_{j+c},\cr h d_{j+c} &= \sum_{i=2}^{p+1} \gamma_i Y^{(i)}_j h^i + {\cal O}(h^{p+2}). } $$

Proof: Use the fact that

$$ \hat d_{j+1} = {\cal O}(h^p), \qquad h \sum_{\ell=1}^j \hat d_{\ell+1} = \cal O(h^p). $$

Then

$$ \eqalign{ q_{j+c} &= \hat q_j e + h d_{j+c} + h A Q_{t+j},\cr \hat q_0 &= 0,\cr \hat q_{j+1} &= \hat q_j + h \hat d_{j+1} + h b^T Q_{j+c},\cr &= h \sum_{\ell=0}^j \hat d_{\ell+1} + h \sum_{\ell=0}^j b^T Q_{\ell+c}. } $$

    ☐

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

$$ g_\ell(t) := {(-1)^{\ell+1}\over\ell!} {\partial^\ell\over\partial y^\ell}f(t,y(t)), \qquad D := \mathop{\rm diag}(c_1,\ldots,c_s) \in\mathbb{R}^{s\times s}. $$

Further

$$ G_\ell(t) := \mathop{\rm diag}\left( g_\ell(t_j+c_1h),\ldots,g_\ell(t_j+c_sh)\right). $$

By Taylor expansion at $t_j$:

$$ G_\ell(t) = g_\ell(t_j) + h D g_l'(t_j) + {1\over2!}h^2 D^2 g_\ell''(t_j) + \cdots \tag{G} $$

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

$$ Q_{j+c} = G_1(t_j) q_{j+c} + G_2(t_j) q_{j+c}^2 + G_3(t_j) q_{j+c}^3 + \cdots. $$

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

$$ Q_{j+c_i} = f(t_j+c_ih,Y_{j+c_i}) - f(t_j+c_ih,y_{j+c_i}) \in\mathbb{R}. $$

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

$$ -f(t_j+c_ih,y_{j+c_i}) = -f(t_j+c_ih,Y_{j+c_i}) + \sum_{\ell=1}^p g_\ell(t_j+c_ih)(Y_{j+c_i}-y_{j+c_i})^\ell + \cdots. $$

Hence

$$ Q_{j+c_i} = \sum_{\ell=1}^p g_\ell(t_j+c_ih) q_{j+c_i}^\ell + \cdots. $$

    ☐

The following two nonlinear equations vanish:

$$ \pmatrix{0\cr 0} = \pmatrix{u\cr v} := \pmatrix{ q_{j+c} - h d_{j+c} - h A Q_{t+c} - {\cal O}(h^p)\cr Q_{t+c} - G_1(t_j) q_{t+c} - G_2(t_j) q_{j+c}^2 - G_3(t_j) q_{j+c}^3 - \cdots } \in \mathbb{R}^{s+1}. $$

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

$$ \pmatrix{q_{t+c}(h)\cr Q_{t+c}(h)} = \pmatrix{ r_2(t_j) h^2 + r_3(t_j) h^3 + \cdots + r_{p-1} h^{p-1} + {\cal O}(h^p)\cr w_2(t_j) h^2 + w_3(t_j) h^3 + \cdots + w_{p-1} h^{p-1} + {\cal O}(h^p) } \tag{T} $$

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:

$$ b^T w_i(t_j) = 0 \in\mathbb{R}, \qquad i=2,\ldots,p-1. $$

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:

$$ \prod_{k=1}^n \left( \sum_{\nu\ge0} a_{\nu k} z_k^\nu \right) = \sum_{|\alpha|\ge0} a_\alpha z^\alpha, $$

with the multiindex notation:

$$ \eqalign{ \alpha &= (\alpha_1,\ldots,\alpha_n)\cr |\alpha| &= \alpha_1 + \cdots + \alpha_n\cr a_\alpha &= a_{1\alpha_1} a_{2\alpha_2} \cdots a_{n\alpha_n}\cr z &= (z_1,\ldots,z_n)\cr z^\alpha &= z_1^{\alpha_1} \cdots z_n^{\alpha_n} } $$

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:

$$ \prod_{k=1}^n \left( \sum_{\nu\ge\beta_k} a_{\nu k} z_k^\nu \right) = \sum_{|\alpha|\ge0} a_{\alpha+\beta} z^{\alpha+\beta}, $$

with the multiindex notation:

$$ \eqalign{ \alpha &= (\alpha_1,\ldots,\alpha_n)\cr \beta &= (\beta_1,\ldots,\beta_n)\cr \alpha+\beta &= (\alpha_1+\beta_1,\ldots,\alpha_n+\beta_n)\cr |\alpha| &= \alpha_1 + \cdots + \alpha_n\cr a_\alpha &= a_{1\alpha_1} a_{2\alpha_2} \cdots a_{n\alpha_n}\cr z &= (z_1,\ldots,z_n)\cr z^\alpha &= z_1^{\alpha_1} \cdots z_n^{\alpha_n} } $$

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

$$ \prod_{k\ge1}^n \left( \sum_{\nu\ge\nu_0} a_{\nu k} z_k^\nu \right) = \sum_{|\alpha|\ge\nu_0} a_\alpha z^{\alpha+\nu_0}, $$

with the infinite multiindex notation:

$$ \eqalign{ \alpha &= (\alpha_1,\alpha_2,\ldots)\cr |\alpha| &= \alpha_1 + \alpha_2 + \cdots\cr a_\alpha &= a_{1\alpha_1} a_{2\alpha_2} \cdots\cr z &= (z_1,z_2,\ldots)\cr z^\alpha &= z_1^{\alpha_1} z_2^{\alpha_2} \cdots } $$

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

$$ z^{\alpha+\nu_0} = z^{\alpha_1+\nu_0} z^{\alpha_2+\nu_0} \cdots $$

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

$$ \prod_{k\ge1} \left( \sum_{\nu\ge0} a_{k\nu} z^\nu \right) = \sum_{|\alpha|\ge0} a_\alpha z^{|\alpha|} $$

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,\ldots,p$:

$$ r_i(t_j) = \gamma_i Y_j^{(i)} + A w_{i-1}(t_j), \qquad w_1(t_j) := 0 \in\mathbb{R}^s, $$

and using the infinite multiindex $\alpha=(\alpha_1,\alpha_2,\ldots)$

$$ \eqalign{ w_i(t_j) &= \sum_{\ell=0}^{i-2} D^\ell {1\over\ell!} g_1^{(\ell)}(t_j) r_{i-\ell}(t_j)\cr &+ \sum_{\ell=0}^{i-4} D^\ell \sum_{\substack{\alpha_1,\alpha_2\ge2\\ \alpha_1+\alpha_2+\ell=i}} {1\over\ell!} g_2^{(\ell)}(t_j) r_{\alpha_1}(t_j) r_{\alpha_2}(t_j)\cr &+ \sum_{\ell=0}^{i-6} D^\ell \sum_{\substack{\alpha_1,\alpha_2,\alpha_3\ge2\\ \alpha_1+\alpha_2+\alpha_3+\ell=i}} {1\over\ell!} g_2^{(\ell)}(t_j) r_{\alpha_1}(t_j) r_{\alpha_2}(t_j) r_{\alpha_3}(t_j)\cr &+ \cdots\cr &= \sum_{k\ge1} \sum_{\ell=0}^{i-2k} D^\ell \sum_{\substack{\alpha\ge2\\ |\alpha|+\ell=i}} {1\over\ell!} g_k^{(\ell)}(t_j) r_\alpha(t_j) . } $$

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

$$ Q_{j+c} = G_1 \left(\sum_\ell r_\ell h^\ell\right) + G_2 \left(\sum_\ell r_\ell h^\ell\right)^2 + G_3 \left(\sum_\ell r_\ell h^\ell\right)^3 + \cdots $$

Now use (G) and Cauchy product formula:

$$ \eqalign{ Q_{j+c} &= \sum_{\nu=1}^p G_\nu q_{j+c}^\nu + {\cal O}(h^p)\cr &= \sum_{\nu=1}^p G_\nu \sum_{|\alpha|\ge0} r_{\alpha_\nu+2} h^{|\alpha|+2} + {\cal O}(h^p)\cr &= \sum_{\nu=1}^p \sum_{|\alpha|\ge0} \sum_{\ell\ge0} {1\over\ell!} D^\ell g_\nu^{(\ell)} r_{\alpha_\nu+2} h^{|\alpha|+\ell+2} + {\cal O}(h^p) } $$

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

$$ \eqalign{ w_i(t_j) &= \sum_{\ell=1}^{m_i} \rho_{i\ell} \alpha_{i\ell} e_{i\ell}(t_j),\cr r_i(t_j) &= \sum_{\ell=0}^{m_i-1} \sigma_{i\ell} \beta_{i\ell} e_{i-1,\ell}(t_j), \qquad e_{i-1,0} := Y^{(i)}. } $$

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

$$ b^T \alpha_{i\ell} = 0, \qquad i=2,\ldots,p, \quad \ell=1,\ldots,m_i. $$

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

Recursion 1:

$$ \eqalign{ R_i^\alpha := &\{\gamma_i\} \cup \left\{Aw \mid w\in W_{i-1}^\alpha\right\}, \qquad W_1^\alpha:= \emptyset, \quad i=2,3,\ldots\cr W_i^\alpha := &\{D^\ell r_1 \mid r_1\in R_{i-\ell}^\alpha, \quad\ell=0,\ldots,i-2\}\cr & \cup \left\{D^\ell r_1\cdot r_2 \mid r_1\in R_{\alpha_1}^\alpha, \quad r_2\in R_{\alpha_2}^\alpha, \quad\ell=0,\ldots,i-4, \quad \alpha_1+\alpha_2+\ell=i \right\}\cr & \cup \left\{D^\ell r_1\cdot r_2\cdot r_3 \mid r_1\in R_{\alpha_1}^\alpha, \: r_2\in R_{\alpha_2}^\alpha, \: r_3\in R_{\alpha_3}^\alpha, \: \ell=0,\ldots,i-6, \: |\alpha|+\ell=i\right\}\cr & \cup \quad \cdots } $$

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\ge3$ if the last stage has order of consistency of $p$, i.e.,

$$ b^T c^{\ell-1} = {1\over\ell}, \qquad\ell=1,\ldots,p, $$

and if

$$ b^T \alpha_{i\ell} = 0, \qquad i=2,\ldots,p-1, \quad \ell=1,\ldots,m_\ell. $$

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.