GitHub

Composing B-series

B-series can be composed, i.e., the result of one B-series can be inserted into another B-series. We saw this operation already in a previous tutorial creating the B-series of a Runge-Kutta method step by step. Here, we present additional examples of composing B-series.

Richardson extrapolation

Richardson extrapolation is a well-known technique to increase the order of accuracy of a rime integration method. Assume we are given a Runge-Kutta method with step size $h$. Let $u^{n+1}_h$ be the solution after one step with step size $h$, and let $u^{n+1}_{2 \times \frac{h}{2}}$ be the solution after two steps with step size $\frac{h}{2}$. The difference between these two approximations can be used to estimate the leading-order error term and adapt the step size accordingly.

Moreover, we can combine the two approximations to obtain a new approximation $u^{n+1}_R$ by extrapolating to step size $h = 0$ as follows:

\[u^{n+1}_R = u^{n+1}_{2 \times \frac{h}{2}} + \frac{ u^{n+1}_{2 \times \frac{h}{2}} - u^{n+1}_h }{ 2^p - 1 },\]

where $p$ is the order of accuracy of the original method. First, we use the classical fourth-order Runge-Kutta method as our base method.

using BSeries

# classical RK4
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]
rk = RungeKuttaMethod(A, b)
series_rk = bseries(rk, 6)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 0
  ⋮                                     => ⋮
order_of_accuracy(series_rk)
4

Next, we compute the B-series of two steps with half the step size.

series_rk_2 = compose(series_rk, series_rk, normalize_stepsize = true)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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//128
  RootedTree{Int64}: [1, 2, 3, 4, 4]    => 13//768
  RootedTree{Int64}: [1, 2, 3, 4, 3]    => 19//768
  RootedTree{Int64}: [1, 2, 3, 4, 2]    => 13//384
  RootedTree{Int64}: [1, 2, 3, 3, 3]    => 19//384
  RootedTree{Int64}: [1, 2, 3, 3, 2]    => 17//256
  RootedTree{Int64}: [1, 2, 3, 2, 3]    => 13//256
  RootedTree{Int64}: [1, 2, 3, 2, 2]    => 77//768
  RootedTree{Int64}: [1, 2, 2, 2, 2]    => 77//384
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 5//4608
  ⋮                                     => ⋮
order_of_accuracy(series_rk_2)
4

Finally, we can combine the two B-series to obtain the Richardson-extrapolated B-series.

p = order_of_accuracy(series_rk)
series_richardson = series_rk_2 + (series_rk_2 - series_rk) / (2^p - 1)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 1//864
  ⋮                                     => ⋮
order_of_accuracy(series_richardson)
5

We can also perform this composition of the methods at the level of the Butcher tableau.

AA = [A zero(A) zero(A)
      zero(A) A/2 zero(A)
      zero(A) repeat(b'/2, length(b)) A/2]
bb = vcat(-b / (2^p - 1),
          b / 2 * (1 + 1 // (2^p - 1)),
          b / 2 * (1 + 1 // (2^p - 1)))
rk_richardson = RungeKuttaMethod(AA, bb)
series_richardson_2 = bseries(rk_richardson, 6)
all(iszero, values(series_richardson - series_richardson_2))
true

Effective order and composition

B-series representing a mapping such as a time integration method can be composed. This operation is equivalent to performing a step with one method after another. Here, we present the example of Butcher's method of effective order 5. This is a fourth-order Runge-Kutta method resulting in a fifth-order method when composed with a special starting and finishing procedure.

First, we set up the B-series of the main method (method "a" in Butcher's paper):

using BSeries

A = [0 0 0 0 0;
     1//5 0 0 0 0;
     0 2//5 0 0 0;
     3//16 0 5//16 0 0;
     1//4 0 -5//4 2 0]
b = [1 // 6, 0, 0, 2 // 3, 1 // 6]
rk_a = RungeKuttaMethod(A, b)
series_a = bseries(rk_a, 6)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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//48
  RootedTree{Int64}: [1, 2, 3, 4, 2]    => 1//30
  RootedTree{Int64}: [1, 2, 3, 3, 3]    => 1//24
  RootedTree{Int64}: [1, 2, 3, 3, 2]    => 1//15
  RootedTree{Int64}: [1, 2, 3, 2, 3]    => 5//96
  RootedTree{Int64}: [1, 2, 3, 2, 2]    => 5//48
  RootedTree{Int64}: [1, 2, 2, 2, 2]    => 5//24
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 0
  ⋮                                     => ⋮

Note that this method is fourth-order accurate:

order_of_accuracy(series_a)
4

Next, we set up the starting procedure (method "b" in Butcher's paper):

A = [0 0 0 0 0;
     1//5 0 0 0 0;
     0 2//5 0 0 0;
     75//64 -9//4 117//64 0 0;
     -37//36 7//3 -3//4 4//9 0]
b = [19 // 144, 0, 25 // 48, 2 // 9, 1 // 8]
rk_b = RungeKuttaMethod(A, b)
series_b = bseries(rk_b, 6)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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]       => 13//320
  RootedTree{Int64}: [1, 2, 3, 3]       => 13//160
  RootedTree{Int64}: [1, 2, 3, 2]       => 121//960
  RootedTree{Int64}: [1, 2, 2, 2]       => 121//480
  RootedTree{Int64}: [1, 2, 3, 4, 5]    => 13//1600
  RootedTree{Int64}: [1, 2, 3, 4, 4]    => 13//800
  RootedTree{Int64}: [1, 2, 3, 4, 3]    => 139//6400
  RootedTree{Int64}: [1, 2, 3, 4, 2]    => 13//400
  RootedTree{Int64}: [1, 2, 3, 3, 3]    => 139//3200
  RootedTree{Int64}: [1, 2, 3, 3, 2]    => 13//200
  RootedTree{Int64}: [1, 2, 3, 2, 3]    => 2003//38400
  RootedTree{Int64}: [1, 2, 3, 2, 2]    => 2003//19200
  RootedTree{Int64}: [1, 2, 2, 2, 2]    => 2003//9600
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 0
  ⋮                                     => ⋮
order_of_accuracy(series_b)
3

Note that this method is only third-order accurate - as is the finishing procedure given by

A = [0 0 0 0 0;
     1//5 0 0 0 0;
     0 2//5 0 0 0;
     161//192 -19//12 287//192 0 0;
     -27//28 19//7 -291//196 36//49 0]
b = [7 // 48, 0, 475 // 1008, 2 // 7, 7 // 72]
rk_c = RungeKuttaMethod(A, b)
series_c = bseries(rk_c, 6)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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]       => 41//960
  RootedTree{Int64}: [1, 2, 3, 3]       => 41//480
  RootedTree{Int64}: [1, 2, 3, 2]       => 119//960
  RootedTree{Int64}: [1, 2, 2, 2]       => 119//480
  RootedTree{Int64}: [1, 2, 3, 4, 5]    => 41//4800
  RootedTree{Int64}: [1, 2, 3, 4, 4]    => 41//2400
  RootedTree{Int64}: [1, 2, 3, 4, 3]    => 463//19200
  RootedTree{Int64}: [1, 2, 3, 4, 2]    => 41//1200
  RootedTree{Int64}: [1, 2, 3, 3, 3]    => 463//9600
  RootedTree{Int64}: [1, 2, 3, 3, 2]    => 41//600
  RootedTree{Int64}: [1, 2, 3, 2, 3]    => 639//12800
  RootedTree{Int64}: [1, 2, 3, 2, 2]    => 639//6400
  RootedTree{Int64}: [1, 2, 2, 2, 2]    => 639//3200
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 0
  ⋮                                     => ⋮
order_of_accuracy(series_c)
3

Finally, we can compose the three methods to obtain

series = compose(series_b, series_a, series_c, normalize_stepsize = true)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 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
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => 121//87480
  ⋮                                     => ⋮

Note that this composition has to be read from left to right. Finally, we check that the resulting series is indeed fifth-order accurate:

order_of_accuracy(series)
5
series - ExactSolution(series)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 entries:
  RootedTree{Int64}: Int64[]            => 0
  RootedTree{Int64}: [1]                => 0
  RootedTree{Int64}: [1, 2]             => 0
  RootedTree{Int64}: [1, 2, 3]          => 0
  RootedTree{Int64}: [1, 2, 2]          => 0
  RootedTree{Int64}: [1, 2, 3, 4]       => 0
  RootedTree{Int64}: [1, 2, 3, 3]       => 0
  RootedTree{Int64}: [1, 2, 3, 2]       => 0
  RootedTree{Int64}: [1, 2, 2, 2]       => 0
  RootedTree{Int64}: [1, 2, 3, 4, 5]    => 0
  RootedTree{Int64}: [1, 2, 3, 4, 4]    => 0
  RootedTree{Int64}: [1, 2, 3, 4, 3]    => 0
  RootedTree{Int64}: [1, 2, 3, 4, 2]    => 0
  RootedTree{Int64}: [1, 2, 3, 3, 3]    => 0
  RootedTree{Int64}: [1, 2, 3, 3, 2]    => 0
  RootedTree{Int64}: [1, 2, 3, 2, 3]    => 0
  RootedTree{Int64}: [1, 2, 3, 2, 2]    => 0
  RootedTree{Int64}: [1, 2, 2, 2, 2]    => 0
  RootedTree{Int64}: [1, 2, 3, 4, 5, 6] => -1//174960
  ⋮                                     => ⋮

References

  • Butcher, J. C. "The effective order of Runge-Kutta methods." In Conference on the numerical solution of differential equations, pp. 133-139. Springer, Berlin, Heidelberg, 1969. DOI: 10.1007/BFb0060019