Original post is here: eklausmeier.goip.de

Below text is based on the results in

Bibliography:

Peter Albrecht.
Michael E. Hosea and LinkedIn.
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:

Let $h$ be our step-size, $t_j = t_0 + jh$, with $j\in\mathbb{N}$.
Let $p\in\mathbb{N}$ be the order of our Runge-Kutta method, see below.
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.
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$
Let $y_{j+c_i}$ be the approximation according to below Runge-Kutta formula.
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} $$

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

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

[markmap]

Runge-Kutta methods #

explicit #

implicit #

semi-implicit #

diagonally implicit #

[/markmap]

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. $$

. 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} $$

. 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:

the solution $Q_{t+c}$ exists and is unique,
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.

. 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$.

. 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} } $$

. 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 #

. 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 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 } $$

. 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.

#Runge-Kutta formula #Peter Albrecht #IVP #ordinary differential equations #Michael Hosea