# How to estimate the stage error for Runge kutta method

Consider an ordinary differential equation (ODE) in the form $$u_t=g(t,u(t))$$ and apply the explicit Runge-Kutta method, as defined by the following Butcher tableau: $$\mathrm{RK}(s,p):\begin{array}{c|cccc} 0 & 0 & 0 & 0 \\ c_1 & a_{1,0} & 0 & 0 &0\\ c_2 & a_{2,0} & a_{2,1} & 0 &0\\ \vdots & \cdots & \cdots & \cdots &0\\ \hline c_s & a_{s,0} & a_{s,1} & \cdots & 0 \end{array},$$

Applying the taylor expansions to obtain $$u(t_{n,i}) = u(t_n) + \sum\limits_{k=1}^{p} \frac{(c_i \tau)^k}{k!}u^{(k)}(t_n) + \mathcal{O}(\tau^{p+1}),$$ $$g(t_{n,j},u(t_{n,j})) = \sum\limits_{k=1}^{p} \frac{(c_j \tau)^{k-1}}{(k-1)!}u^{(k)}(t_n) + \mathcal{O}(\tau^{p})$$ where $$\tau = t_{n+1}- t_n$$ and $$u^{(k)}(t_n) = \frac{d^k u}{dt^{k}}|_{t_n}.$$

The numerical scheme is $$u_{n,i} = u^n + \tau \sum\limits_{j=0}^{i-1}a_{i,j} g(t_{n,j},u_{n,j}), \quad i=1,...,s.$$ Assuming $$u^n = u(t_n)$$, inserting the exact solutions in the numerical scheme $$u(t_n+c_i \tau) = u(t_{n,i}) = u(t_n) + \tau \sum\limits_{j=0}^{i-1} a_{i,j}g(t_{n,j},u(t_{n,j})) + \delta_{n,i}, \quad i=1,...,s.$$ where $$\delta_{n,i} = \sum\limits_{k=1}^{p} \tau^k (\frac{c_i^k}{k!} - \frac{1}{(k-1)!} \sum\limits_{j=0}^{i-1}a_{i,j} c_j^{k-1})u^{(k)}(t_n) + \mathcal{O}(\tau^{p+1}),\quad i=1,...,s.$$

Let $$e^{n,i} = u(t_{n,i}) - u_{n,i}$$. What I want to obtain is the $$e^{n+1}-e^{n} = e^{n,s}-e^{n} = \mathcal{O}(\tau^{p+1})$$, and the error equation can be written as $$e_{n,i} = e^n + \tau \sum\limits_{j=0}^{i-1} a_{i,j} \eta_{n,j} -\delta_{n,i}.$$ where $$\eta_{n,j} = g(t_{n,j},u(t_{n,j}) - g(t_{n,j},u_{n,j}))$$. Due to the coefficients of the last step satisfying the condition $$(\frac{c_s^k}{k!} - \frac{1}{(k-1)!} \sum\limits_{j=0}^{s-1}a_{s,j} c_j^{k-1})=0, k=1,...,p.$$ I arrive at the expression $$e_{n,s} = e^n + \tau \sum\limits_{j=0}^{s-1} a_{i,j} \eta_{n,j} + \mathcal{O}(\tau^{p+1}).$$ How to prove the second term on the right hand $$\tau \sum\limits_{j=0}^{s-1} a_{i,j} \eta_{n,j} = \mathcal{O}(\tau^{p+1})$$. It seems more complex to estimate the stage error $$e^{n,i}$$,

It seems to be related to the order condition of Runge Kutta method. Is there any literature on this, namely stage error(maybe)? Thank you in advance.

• Yes, in the design of methods the order of the stages is used as one of the tools to simplify the system of order conditions. This is discussed in the design papers for higher-order methods, like the 10 (or 12?) stage (7)8 methods. // I'm uncertain if our approach has any merit, the stages are carefully constructed to lie aside the exact path, so that their error terms cancel in the linear combinations. Using the derivatives of the exact path may miss some higher-order cancellations. Commented Dec 13, 2023 at 10:41
• There are a number of textbooks that provide what you're looking for. The most widely used are either vol. I of Hairer, Norsett, & Wanner, or the text of Butcher. Commented Dec 14, 2023 at 4:33
• To me it's not 100% clear what you're asking, probably just due to language issues (particularly, I don't understand the sentence containing your last equation). Commented Dec 14, 2023 at 4:35
• I don't understand using the specific properties of $g$. In order to find the order of a scheme, it needs to have this order regardless of what $g$ you choose. I will also say that I don't understand what exactly this condition is expressing, among other reasons because I don't understand what $e^{n,j}$ is supposed to mean. Commented Dec 14, 2023 at 6:09
• @DavidKetcheson Thx u, I have modified my question and it now expresses my problem clearly. Commented Dec 14, 2023 at 8:28

The extension of Taylor expansions to Runge-Kutta methods are the B-series based on single-rooted trees as formalized by Butcher. Such a series has the general form $$\Delta y=\Psi(y,h)=\sum_{\tau\in trees}h^{|\tau|}\frac{\psi(\tau)}{\sigma(\tau)}F(\tau).$$ This all for some fixed right side in $$y'=f(y)$$, autonomous for simplicity. The coefficients of the series are $$\psi(\tau)$$, having the trees as arguments. $$\sigma(\tau)$$ is a constant of the tree representing its combinatorial structure, and $$F(\tau)$$ is a nested derivative construct according to the tree structure.

If one removes the root of the tree, it splits into sub-trees with multiplicity, written as $$\tau = [\tau_1^{m_1},...,\tau_k^{m_k}].$$ The recursive parts of the formula are then defined by $$m(\tau)=m_1+...+m_k\\ |\tau|=1+m_1|\tau_1|+...m_k|\tau_k|\\ \sigma(\tau)=\prod m_i!\sigma(\tau_i)^{m_i}\\ F(\tau)=f^{(m(\tau))}(y)[F(\tau_1)^{\odot m_1},...,F(\tau_k)^{\odot m_k}].$$ The last is the derivative of order $$m(\tau)$$ evaluated at $$y$$ and applied in the directions of the derivative expressions for the sub-trees, with copies according to their multiplicity. Note that the derivative is a vector-valued symmetric multi-linear form (or a tensor), thus the symmetric tensor product sign to denote the type of power.

The exact solution has a B-series with $$\phi(\tau)=\frac1{\gamma(\tau)}.$$ The new tree-constant collects the factors due to the integration of the powers of $$h$$, that is $$\gamma(\tau)=|\tau|\gamma(\tau_1)^{m_1}...\gamma(\tau_k)^{m_k}$$

The point I wanted to demonstrate was that you get multiple term for the same order, and the coefficients in each term have to be equal in exact solution and numerical approximation. Some of the order conditions are redundant, so that the order conditions in higher orders can be satisfied with less stages than pure counting suggests.

For your question, you can get multiple measures for the intermediate orders. You get the error order of the midpoints $$y_{n,i}$$ as approximations of $$y(x_n+c_ih)$$, or you could consider the order of $$k_i=hf(y_{n,i})$$ as when used in an Euler step.

For RK4 the leading error terms turn out to be

y + k1 = y(x+h) - 1/2*F([2,0])*h^2 - 1/6*F([3,0])*h^3 - 1/6*F([3,1])*h^3

y_n1 = y(x+h/2) - 1/8*F([2,0])*h^2 - 1/48*F([3,0])*h^3 - 1/48*F([3,1])*h^3

y + k2 = y(x+h) - 1/6*F([3,0])*h^3 - 1/24*F([3,1])*h^3

y_n2 = y(x+h/2) + 1/8*F([2,0])*h^2 - 1/48*F([3,0])*h^3 + 1/24*F([3,1])*h^3

y + k3 = y(x+h) + 1/12*F([3,0])*h^3 - 1/24*F([3,1])*h^3

y_n3 = y(x+h)   + 1/12*F([3,0])*h^3 - 1/24*F([3,1])*h^3

y + k4 = y(x+h) + 1/2*F([2,0])*h^2 + 1/3*F([3,0])*h^3 + 1/3*F([3,1])*h^3 + 5/24*F([4,0])*h^4 + 1/12*F([4,1])*h^4 + 3/16*F([4,2])*h^4 + 1/8*F([4,3])*h^4

========================

y_next = y + F([1,0])*h^1 + 1/2*F([2,0])*h^2 + 1/6*F([3,0])*h^3 + 1/6*F([3,1])*h^3 + 1/24*F([4,0])*h^4 + 1/24*F([4,1])*h^4 + 1/16*F([4,2])*h^4 + 1/24*F([4,3])*h^4 + 1/96*F([5,1])*h^5 + 1/96*F([5,2])*h^5 + 1/144*F([5,3])*h^5 + 1/48*F([5,4])*h^5 + 1/64*F([5,5])*h^5 + 1/32*F([5,6])*h^5 + 5/288*F([5,7])*h^5 + 5/576*F([5,8])*h^5

= y(x+h) - 1/120*F([5,0])*h^5 + 1/480*F([5,1])*h^5 - 1/480*F([5,2])*h^5 - 1/720*F([5,3])*h^5 + 1/240*F([5,4])*h^5 - 1/960*F([5,5])*h^5 + 1/160*F([5,6])*h^5 + 1/1440*F([5,7])*h^5 + 1/2880*F([5,8])*h^5


As one can see, the individual orders of the stages give little insight into the final order of the method.

The trees are

done
T[ 1, 0] = []
T[ 2, 0] = [[]] = [T[1, 0]]
T[ 3, 0] = [[[]]] = [T[2, 0]]
T[ 3, 1] = [[][]] = [2*T[1, 0]]
T[ 4, 0] = [[[[]]]] = [T[3, 0]]
T[ 4, 1] = [[[][]]] = [T[3, 1]]
T[ 4, 2] = [[[]][]] =  [T[2, 0], T[1, 0]]
T[ 4, 3] = [[][][]] = [3*T[1, 0]]
T[ 5, 0] = [[[[[]]]]] = [T[4, 0]]
T[ 5, 1] = [[[[][]]]] = [T[4, 1]]
T[ 5, 2] = [[[[]][]]] = [T[4, 2]]
T[ 5, 3] = [[[][][]]] = [T[4, 3]]
T[ 5, 4] = [[[[]]][]] = [T[3, 0]], T[1, 0]]
T[ 5, 5] = [[[][]][]] = [T[3, 1], T[1, 0]]
T[ 5, 6] = [[[]][[]]] = [2*T[2, 0]]
T[ 5, 7] = [[[]][][]] = [T[2, 0], 2*T[1, 0]]
T[ 5, 8] = [[][][][]] = [[4*T[1, 0]]]


[] denotes the rooted tree that is only the root node, or the leafs in trees of higher order.

• Thank you for your detailed answer. It is very inspiring to me. I will learn more about the Runge-Kutta tree root theorem and order condition. Commented Dec 18, 2023 at 2:07

What you're hoping to prove is (in general) not true.

In words, you want to show that for a method of order $$p$$, each stage of the Runge-Kutta method also approximates the solution (at $$t_n + c_i \Delta t$$) to order $$p$$. But this is not the case; for instance, in the classical 4th-order method of Kutta, none of the intermediate stages are 4th-order accurate. Rather, the leading errors in the stages are made to cancel at the end of the step, in order to obtain a 4th-order approximation.

The accuracy of the stages does play a crucial role in some situations -- particularly in application to stiff ODEs, or PDEs with time-dependent boundary conditions. For these problems, the observed rate of convergence can be less than the classical order and at worst equal to the lowest order of the stages; this is known as order reduction.