B-series basics

In this tutorial we use BSeries.jl to investigate error expansions for Runge-Kutta (RK) methods, generically or when applied to a specific ordinary differential equation (ODE).

# Load the packages we will use.
# These must first be installed using: import Pkg; Pkg.add("package_name")
using BSeries
using Latexify  # Only needed for some pretty-printing cells below using `latexify`
import SymPyPythonCall; sp = SymPyPythonCall;

SymPyPythonCall

B-series for a generic ODE

First we specify the Butcher coefficients of the RK method. This can include symbolic expressions and parameterized families of methods. Here is a generic 2-stage, 2nd-order method:

α = sp.symbols("α", real = true)
A = [0 0; 1/(2*α) 0]; b = [1-α, α]; c = [0, 1/(2*α)]
coeffs2 = bseries(A, b, c, 3)
latexify(coeffs2, cdot=false)

\[F_{f}\mathopen{}\left( \varnothing \right)\mathclose{} + h F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{h^{2}}{2} F_{f}\mathopen{}\left( \rootedtree[.[.]] \right)\mathclose{} + \frac{h^{3}}{8 \alpha} F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}\]

The rooted trees are printed as nested lists, essentially in the form used in Butcher's book. The rooted trees written in this way can be rendered in LaTeX using the package forest; unfortunately, there is no easy way to render them in the browser. Nevertheless, you can render them using LaTeX with an appropriate preamble, see the docstring of RootedTrees.latexify. The rendered output looks like this:

We have generated the B-series up to terms of order $h^3$. The terms $F_f()$ represent elementary differentials, which are products of derivatives of the ODE right-hand side. Since we haven't specified an ODE, these are indicated simply by the associated rooted tree.

Here is a B-series for the classical 4th-order method, expanded up to 5th-order terms:

A = [0 0 0 0; 1//2 0 0 0; 0 1//2 0 0; 0 0 1 0];
b = [1//6, 1//3, 1//3, 1//6];
c = [0, 1//2, 1//2, 1];

coeffs4 = bseries(A, b, c, 5)
latexify(coeffs4, cdot=false)

\[F_{f}\mathopen{}\left( \varnothing \right)\mathclose{} + h F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{1}{2} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.]] \right)\mathclose{} + \frac{1}{6} h^{3} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{6} h^{3} F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]]]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.]]] \right)\mathclose{} + \frac{1}{8} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.]] \right)\mathclose{} + \frac{1}{96} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.][.]]]] \right)\mathclose{} + \frac{1}{48} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]][.]]] \right)\mathclose{} + \frac{1}{24} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]]][.]] \right)\mathclose{} + \frac{1}{144} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.][.]]] \right)\mathclose{} + \frac{1}{32} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.]][.]] \right)\mathclose{} + \frac{1}{32} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.[.]]] \right)\mathclose{} + \frac{5}{96} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.][.]] \right)\mathclose{} + \frac{5}{576} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.][.]] \right)\mathclose{}\]

We can also print out the B-series coefficients this way:

coeffs4

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
  RootedTree{Int64}: Int64[]         => 1
  RootedTree{Int64}: [1]             => 1
  RootedTree{Int64}: [1, 2]          => 1//2
  RootedTree{Int64}: [1, 2, 3]       => 1//6
  RootedTree{Int64}: [1, 2, 2]       => 1//3
  RootedTree{Int64}: [1, 2, 3, 4]    => 1//24
  RootedTree{Int64}: [1, 2, 3, 3]    => 1//12
  RootedTree{Int64}: [1, 2, 3, 2]    => 1//8
  RootedTree{Int64}: [1, 2, 2, 2]    => 1//4
  RootedTree{Int64}: [1, 2, 3, 4, 5] => 0
  RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//48
  RootedTree{Int64}: [1, 2, 3, 4, 3] => 1//48
  RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//24
  RootedTree{Int64}: [1, 2, 3, 3, 3] => 1//24
  RootedTree{Int64}: [1, 2, 3, 3, 2] => 1//16
  RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//16
  RootedTree{Int64}: [1, 2, 3, 2, 2] => 5//48
  RootedTree{Int64}: [1, 2, 2, 2, 2] => 5//24

In this form, the rooted trees are printed as level sequences. The corresponding coefficients are on the right.

You can use the function RootedTrees.set_printing_style to change the printing style globally. For example, you can use the notation of Butcher as follows.

RootedTrees.set_printing_style("butcher")
coeffs4

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
  ∅         => 1
  τ         => 1
  [τ]       => 1//2
  [[τ]]     => 1//6
  [τ²]      => 1//3
  [[[τ]]]   => 1//24
  [[τ²]]    => 1//12
  [[τ]τ]    => 1//8
  [τ³]      => 1//4
  [[[[τ]]]] => 0
  [[[τ²]]]  => 1//48
  [[[τ]τ]]  => 1//48
  [[[τ]]τ]  => 1//24
  [[τ³]]    => 1//24
  [[τ²]τ]   => 1//16
  [[τ][τ]]  => 1//16
  [[τ]τ²]   => 5//48
  [τ⁴]      => 5//24

To use the level sequence representation, you need to change the printing style again.

RootedTrees.set_printing_style("sequence")
coeffs4

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
  RootedTree{Int64}: Int64[]         => 1
  RootedTree{Int64}: [1]             => 1
  RootedTree{Int64}: [1, 2]          => 1//2
  RootedTree{Int64}: [1, 2, 3]       => 1//6
  RootedTree{Int64}: [1, 2, 2]       => 1//3
  RootedTree{Int64}: [1, 2, 3, 4]    => 1//24
  RootedTree{Int64}: [1, 2, 3, 3]    => 1//12
  RootedTree{Int64}: [1, 2, 3, 2]    => 1//8
  RootedTree{Int64}: [1, 2, 2, 2]    => 1//4
  RootedTree{Int64}: [1, 2, 3, 4, 5] => 0
  RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//48
  RootedTree{Int64}: [1, 2, 3, 4, 3] => 1//48
  RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//24
  RootedTree{Int64}: [1, 2, 3, 3, 3] => 1//24
  RootedTree{Int64}: [1, 2, 3, 3, 2] => 1//16
  RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//16
  RootedTree{Int64}: [1, 2, 3, 2, 2] => 5//48
  RootedTree{Int64}: [1, 2, 2, 2, 2] => 5//24

Exact series and local error

We can also get the B-series of the exact solution:

coeffs_ex = ExactSolution(coeffs4)

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
  RootedTree{Int64}: Int64[]         => 1
  RootedTree{Int64}: [1]             => 1
  RootedTree{Int64}: [1, 2]          => 1//2
  RootedTree{Int64}: [1, 2, 3]       => 1//6
  RootedTree{Int64}: [1, 2, 2]       => 1//3
  RootedTree{Int64}: [1, 2, 3, 4]    => 1//24
  RootedTree{Int64}: [1, 2, 3, 3]    => 1//12
  RootedTree{Int64}: [1, 2, 3, 2]    => 1//8
  RootedTree{Int64}: [1, 2, 2, 2]    => 1//4
  RootedTree{Int64}: [1, 2, 3, 4, 5] => 1//120
  RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//60
  RootedTree{Int64}: [1, 2, 3, 4, 3] => 1//40
  RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//30
  RootedTree{Int64}: [1, 2, 3, 3, 3] => 1//20
  RootedTree{Int64}: [1, 2, 3, 3, 2] => 1//15
  RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//20
  RootedTree{Int64}: [1, 2, 3, 2, 2] => 1//10
  RootedTree{Int64}: [1, 2, 2, 2, 2] => 1//5

latexify(coeffs_ex, cdot=false)

\[F_{f}\mathopen{}\left( \varnothing \right)\mathclose{} + h F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{1}{2} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.]] \right)\mathclose{} + \frac{1}{6} h^{3} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{6} h^{3} F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]]]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.]]] \right)\mathclose{} + \frac{1}{8} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.]] \right)\mathclose{} + \frac{1}{24} h^{4} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.]] \right)\mathclose{} + \frac{1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.[.]]]]] \right)\mathclose{} + \frac{1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.][.]]]] \right)\mathclose{} + \frac{1}{40} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]][.]]] \right)\mathclose{} + \frac{1}{30} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]]][.]] \right)\mathclose{} + \frac{1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.][.]]] \right)\mathclose{} + \frac{1}{30} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.]][.]] \right)\mathclose{} + \frac{1}{40} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.[.]]] \right)\mathclose{} + \frac{1}{20} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.][.]] \right)\mathclose{} + \frac{1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.][.]] \right)\mathclose{}\]

We can find the local error by subtracting the exact solution B-series from the RK method B-series:

latexify(coeffs4-coeffs_ex, cdot=false)

\[\frac{-1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.[.]]]]] \right)\mathclose{} + \frac{1}{480} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.][.]]]] \right)\mathclose{} + \frac{-1}{240} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]][.]]] \right)\mathclose{} + \frac{1}{120} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.[.]]][.]] \right)\mathclose{} + \frac{-1}{720} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.][.]]] \right)\mathclose{} + \frac{-1}{480} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.][.]][.]] \right)\mathclose{} + \frac{1}{160} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.[.]]] \right)\mathclose{} + \frac{1}{480} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.[.]][.][.]] \right)\mathclose{} + \frac{1}{2880} h^{5} F_{f}\mathopen{}\left( \rootedtree[.[.][.][.][.]] \right)\mathclose{}\]

This confirms that the method is of 4th order, since all terms involving smaller powers of $h$ vanish exactly. We don't see the $h^6$ and higher order terms since we only generated the truncated B-series up to 5th order. We can also obtain the order of accuracy awithout comparing the coefficients to the exact solution manually:

order_of_accuracy(coeffs4)

4

For the 2nd-order method, we get

order_of_accuracy(coeffs2)

2

with the following leading error terms:

latexify(coeffs2-coeffs_ex, cdot=false)

\[\frac{ - h^{3}}{6} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + h^{3} \left( \frac{-1}{6} + \frac{1}{8 \alpha} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}\]

This confirms again the accuracy of the method, and shows us that we can eliminate one of the leading error terms completely if we take $\alpha=3/4$ (this is known as Ralston's method, or sometimes as Heun's method).

B-series for a specific ODE

Next, let us define an ODE. We'll consider the Prothero-Robinson problem:

\[ y'(t) = \lambda(y-\sin(t)) + \cos(t).\]

For a non-autonomous ODE like this, it's convenient to rewrite the problem in autonomous form. We set $u=[y,t]^T$ and

\[\begin{align} u_1'(t) & = \lambda(u_1 - \sin(u_2)) + \cos(u_2) \\ u_2'(t) & = t. \end{align}\]

λ = sp.symbols("λ", real = true)
y, t = sp.symbols("y t", real = true)
h = sp.symbols("h", real = true)

u = [y, t]
ff = [λ*(u[1]-sin(t))+cos(t), 1]

\[\left[\begin{smallmatrix}λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\\1\end{smallmatrix}\right]\]

Finally, we get the B-Series for our RK method applied to our ODE:

evaluate(ff, u, h, coeffs4)[1]

\[\frac{h^{5} λ^{2} \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{96} + \frac{h^{5} λ \left(λ \cos{\left(t \right)} + \sin{\left(t \right)}\right)}{144} + \frac{5 h^{5} \left(- λ \sin{\left(t \right)} + \cos{\left(t \right)}\right)}{576} + \frac{h^{4} λ^{2} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{24} + \frac{h^{4} λ \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{24} + \frac{h^{4} \left(λ \cos{\left(t \right)} + \sin{\left(t \right)}\right)}{24} + \frac{h^{3} λ \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{6} + \frac{h^{3} \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{6} + \frac{h^{2} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{2} + h \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) + y\]

Notice that the series is truncated at the same order that we specified when we initially generated it from the RK coefficients.

Here's the B-Series for the exact solution of the same ODE:

evaluate(ff, u, h, coeffs_ex)[1]

\[\frac{h^{5} λ^{3} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{120} + \frac{h^{5} λ^{2} \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{120} + \frac{h^{5} λ \left(λ \cos{\left(t \right)} + \sin{\left(t \right)}\right)}{120} + \frac{h^{5} \left(- λ \sin{\left(t \right)} + \cos{\left(t \right)}\right)}{120} + \frac{h^{4} λ^{2} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{24} + \frac{h^{4} λ \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{24} + \frac{h^{4} \left(λ \cos{\left(t \right)} + \sin{\left(t \right)}\right)}{24} + \frac{h^{3} λ \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{6} + \frac{h^{3} \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{6} + \frac{h^{2} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{2} + h \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) + y\]

And their difference, which is the local error:

expr = sp.simplify(evaluate(ff, u, h, coeffs4) - evaluate(ff, u, h,coeffs_ex))[1]

\[- \frac{h^{5} λ^{3} \left(λ \left(λ \left(y - \sin{\left(t \right)}\right) + \cos{\left(t \right)}\right) - λ \cos{\left(t \right)} - \sin{\left(t \right)}\right)}{120} + \frac{h^{5} λ^{2} \left(λ \sin{\left(t \right)} - \cos{\left(t \right)}\right)}{480} - \frac{h^{5} λ \left(λ \cos{\left(t \right)} + \sin{\left(t \right)}\right)}{720} + \frac{h^{5} \left(- λ \sin{\left(t \right)} + \cos{\left(t \right)}\right)}{2880}\]

B-series for a generic RK method

We can also examine just the elementary differentials, without specifying a RK method:

elementary_differentials(ff, u, 5)

OrderedDict{RootedTree{Int64, Vector{Int64}}, Vector{SymPyCore.Sym{PythonCall.Core.Py}}} with 18 entries:
  RootedTree{Int64}: Int64[]         => [y, t]
  RootedTree{Int64}: [1]             => [λ*(y - sin(t)) + cos(t), 1]
  RootedTree{Int64}: [1, 2]          => [λ*(λ*(y - sin(t)) + cos(t)) - λ*cos(t)…
  RootedTree{Int64}: [1, 2, 3]       => [λ*(λ*(λ*(y - sin(t)) + cos(t)) - λ*cos…
  RootedTree{Int64}: [1, 2, 2]       => [λ*sin(t) - cos(t), 0]
  RootedTree{Int64}: [1, 2, 3, 4]    => [λ^2*(λ*(λ*(y - sin(t)) + cos(t)) - λ*c…
  RootedTree{Int64}: [1, 2, 3, 3]    => [λ*(λ*sin(t) - cos(t)), 0]
  RootedTree{Int64}: [1, 2, 3, 2]    => [0, 0]
  RootedTree{Int64}: [1, 2, 2, 2]    => [λ*cos(t) + sin(t), 0]
  RootedTree{Int64}: [1, 2, 3, 4, 5] => [λ^3*(λ*(λ*(y - sin(t)) + cos(t)) - λ*c…
  RootedTree{Int64}: [1, 2, 3, 4, 4] => [λ^2*(λ*sin(t) - cos(t)), 0]
  RootedTree{Int64}: [1, 2, 3, 4, 3] => [0, 0]
  RootedTree{Int64}: [1, 2, 3, 4, 2] => [0, 0]
  RootedTree{Int64}: [1, 2, 3, 3, 3] => [λ*(λ*cos(t) + sin(t)), 0]
  RootedTree{Int64}: [1, 2, 3, 3, 2] => [0, 0]
  RootedTree{Int64}: [1, 2, 3, 2, 3] => [0, 0]
  RootedTree{Int64}: [1, 2, 3, 2, 2] => [0, 0]
  RootedTree{Int64}: [1, 2, 2, 2, 2] => [-λ*sin(t) + cos(t), 0]