I've been learning about butcher tables and am having some difficulty understanding how to translate them when it comes to implicit methods. Specifically, I'm looking at backwards Euler:

\begin{array} {c|c} 1&1 \\ \hline & 1 \end{array}

Now generally, I am aware of how implicit methods work. A to-be-determined value lies on both sides of the equal sign and must be solved for algebraically:

$y_{n+1} = y_n + hf(t_{n+1}, y_{n+1})$

But what I am having a hard time figuring out is the notation when it comes to finding $k_i$ values from a butcher table in the implicit case.

What I got is this:

$$k_1 = hf(t_n + h, y + hk_1)\\ y_{n+1} = y_n + k_1$$

Am I to assume that $y + hk_1 = y_{n+1}$? If that is the case, how am I to understand the value of $k_1$ algorithmic-ally? It is an unknown, and appears twice. I am trying to learn how to turn Butcher tableaus into code, and this one is confusing me. The higher order explicit methods make sense to me though!


To fix notation, denote the Butcher tableau by $$ \begin{array}{c|c} c & A\\ \hline & b^T \end{array} $$ where $b$ and $c$ are vectors of length $s$ (the number of stages) and $A$ is a $s \times s$ matrix.

Consider the ODE $$ y' = f(t, y) $$ and suppose that $y_n$ is a given approximation to $y$ at time $t = t_n$.

A general Runge-Kutta method (explicit or implicit) can be written as \begin{align} k_i &= f\left( t_n + \Delta t c_i, y_n + \Delta t \sum_{j=1}^s a_{ij} k_j \right), \tag{1} \\ y_{n+1} &= y_n + \Delta t \sum_{i=1}^s b_i k_i. \tag{2} \end{align} Advancing from $y_n$ to $y_{n+1}$ first requires finding the values $k_i$ for $1 \leq i \leq s$ using equation (1), and then computing $y_{n+1}$ using (2).

Note that if the method is explicit, then the matrix $A$ is strictly lower triangular, and each $k_{i}$ can be computed using $k_1, \ldots, k_{i-1}$ without solving a system of equations. On the other hand, if the method is implicit, then finding each $k_i$ will require solving a system of equations.

For concreteness, you give the example of backward Euler, which has the Butcher tableau $$ \begin{array}{c|c} 1 & 1\\ \hline & 1 \end{array} $$ There is only one stage, so substituting the values $c_1 = 1$ and $a_{11} = 1$ into equation (1), we obtain the following equation for $k_1$: $$ k_1 = f(t_n + \Delta t, y_n + \Delta t k_1). $$ As you mention, because $k_1$ also appears on the right-hand side of this equation, it will require solving an algebraic equation.

Once you have found $k_1$, you can compute $y_{n+1}$ using (2), which, after substituting $b_1 = 1$ gives $$ y_{n+1} = y_n + \Delta t k_1. $$

  • $\begingroup$ How do you solve for it algebraically in an algorithmic sort of way? Isn't it going to depend entirely on the ODE and the initial conditions? For example, I tried $y' = y$ and got $k_n = \frac{1}{y_{n-1} - h}$, where h is my stepsize. Is it the job of the programmer to figure out this relationship for each new ODE he is trying to solve, or is there a way to automate this process? $\endgroup$ Nov 1 '19 at 0:05
  • $\begingroup$ For $y'= y$ you will get $k_1 = \frac{y_n}{1-\Delta t}$. Algorithmically, you can try to use a general-purpose solver to compute the solution. For example, you could use Matlab's fsolve where the objective function is $k_1 - f(t_n + \Delta t, y_n + \Delta t k_1)$. For more complicated examples it will require more work on the part of the programmer. $\endgroup$
    – Will P.
    Nov 1 '19 at 0:17
  • $\begingroup$ The whole question is because I am trying to write my own code rather than using a generic solver. $\endgroup$ Nov 1 '19 at 0:22
  • $\begingroup$ It will depend on your ODE. Is it a scalar equation or is it a system? Is it linear or nonlinear? You can try to implement methods like Newton's method or bisection to solve the nonlinear equations. $\endgroup$
    – Will P.
    Nov 1 '19 at 0:42

Your Answer

