Creating B-series
We have already seen some ways of creating B-series in the basic tutorial. In this tutorial, we look at ways to obtain the B-series of different time integration methods.
B-series for Runge-Kutta methods
BSeries.jl and RootedTrees.jl provide the type RungeKuttaMethod as wrapper of Butcher coefficients A, b, c of Runge-Kutta methods. For example, you can create the classical explicit, fourth-order Runge-Kutta method as follows.
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)RungeKuttaMethod{Rational{Int64}} with
A: 4×4 Matrix{Rational{Int64}}:
0//1 0//1 0//1 0//1
1//2 0//1 0//1 0//1
0//1 1//2 0//1 0//1
0//1 0//1 1//1 0//1
b: 4-element Vector{Rational{Int64}}:
1//6
1//3
1//3
1//6
c: 4-element Vector{Rational{Int64}}:
0//1
1//2
1//2
1//1
Instead of passing the Butcher coefficients explicitly to bseries, you can also pass the wrapper struct rk. In fact, this is the preferred method.
series = bseries(rk, 5)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
RootedTree{Int64}: Int64[] => 1//1
RootedTree{Int64}: [1] => 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//1
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//24We can check that the classical Runge-Kutta method is indeed fourth-order accurate.
series - ExactSolution(series)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
RootedTree{Int64}: Int64[] => 0//1
RootedTree{Int64}: [1] => 0//1
RootedTree{Int64}: [1, 2] => 0//1
RootedTree{Int64}: [1, 2, 3] => 0//1
RootedTree{Int64}: [1, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4] => 0//1
RootedTree{Int64}: [1, 2, 3, 3] => 0//1
RootedTree{Int64}: [1, 2, 3, 2] => 0//1
RootedTree{Int64}: [1, 2, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 5] => -1//120
RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//240
RootedTree{Int64}: [1, 2, 3, 4, 3] => -1//240
RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//120
RootedTree{Int64}: [1, 2, 3, 3, 3] => -1//120
RootedTree{Int64}: [1, 2, 3, 3, 2] => -1//240
RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//80
RootedTree{Int64}: [1, 2, 3, 2, 2] => 1//240
RootedTree{Int64}: [1, 2, 2, 2, 2] => 1//120order_of_accuracy(series)4B-series for additive Runge-Kutta methods
BSeries.jl and RootedTrees.jl also support additive Runge-Kutta methods via the wrapper AdditiveRungeKuttaMethod. For example, we can write the Störmer-Verlet method as additive Runge-Kutta method following Table II.2.1 of Hairer, Lubich, and Wanner (2002).
using BSeries
As = [
[0 0; 1//2 1//2],
[1//2 0; 1//2 0],
]
bs = [
[1 // 2, 1 // 2],
[1 // 2, 1 // 2],
]
ark = AdditiveRungeKuttaMethod(As, bs)AdditiveRungeKuttaMethod{Rational{Int64}} with methods
1. RungeKuttaMethod{Rational{Int64}} with
A: 2×2 Matrix{Rational{Int64}}:
0//1 0//1
1//2 1//2
b: 2-element Vector{Rational{Int64}}:
1//2
1//2
c: 2-element Vector{Rational{Int64}}:
0//1
1//1
2. RungeKuttaMethod{Rational{Int64}} with
A: 2×2 Matrix{Rational{Int64}}:
1//2 0//1
1//2 0//1
b: 2-element Vector{Rational{Int64}}:
1//2
1//2
c: 2-element Vector{Rational{Int64}}:
1//2
1//2
We can create the B-series as usual, truncated to order 3.
series = bseries(ark, 3)TruncatedBSeries{BicoloredRootedTree{Int64, Vector{Int64}, Vector{Bool}}, Rational{Int64}} with 21 entries:
ColoredRootedTree{Int64}: (Int64[], Bool[]) => 1//1
ColoredRootedTree{Int64}: ([1], Bool[0]) => 1//1
ColoredRootedTree{Int64}: ([1], Bool[1]) => 1//1
ColoredRootedTree{Int64}: ([1, 2], Bool[0, 0]) => 1//2
ColoredRootedTree{Int64}: ([1, 2], Bool[1, 0]) => 1//2
ColoredRootedTree{Int64}: ([1, 2], Bool[0, 1]) => 1//2
ColoredRootedTree{Int64}: ([1, 2], Bool[1, 1]) => 1//2
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 0, 0]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 0, 0]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 1, 0]) => 0//1
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 1, 0]) => 0//1
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 0, 1]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 0, 1]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 1, 1]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 1, 1]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[0, 0, 0]) => 1//2
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[1, 0, 0]) => 1//2
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[0, 1, 0]) => 1//4
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[1, 1, 0]) => 1//4
⋮ => ⋮For additive Runge-Kutta methods like this, we use colored rooted trees. Again, we can check the order of accuracy by comparing the coefficients to the exact solution.
series - ExactSolution(series)TruncatedBSeries{BicoloredRootedTree{Int64, Vector{Int64}, Vector{Bool}}, Rational{Int64}} with 21 entries:
ColoredRootedTree{Int64}: (Int64[], Bool[]) => 0//1
ColoredRootedTree{Int64}: ([1], Bool[0]) => 0//1
ColoredRootedTree{Int64}: ([1], Bool[1]) => 0//1
ColoredRootedTree{Int64}: ([1, 2], Bool[0, 0]) => 0//1
ColoredRootedTree{Int64}: ([1, 2], Bool[1, 0]) => 0//1
ColoredRootedTree{Int64}: ([1, 2], Bool[0, 1]) => 0//1
ColoredRootedTree{Int64}: ([1, 2], Bool[1, 1]) => 0//1
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 0, 0]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 0, 0]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 1, 0]) => -1//6
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 1, 0]) => -1//6
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 0, 1]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 0, 1]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[0, 1, 1]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 3], Bool[1, 1, 1]) => 1//12
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[0, 0, 0]) => 1//6
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[1, 0, 0]) => 1//6
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[0, 1, 0]) => -1//12
ColoredRootedTree{Int64}: ([1, 2, 2], Bool[1, 1, 0]) => -1//12
⋮ => ⋮order_of_accuracy(series)2We can also create LaTeX code for this B-series as follows.
using Latexify
latexify(series)\[F_{f}\mathopen{}\left( \varnothing \right)\mathclose{} + h \cdot F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + h \cdot F_{f}\mathopen{}\left( \rootedtree[o] \right)\mathclose{} + \frac{1}{2} \cdot h^{2} \cdot F_{f}\mathopen{}\left( \rootedtree[.[.]] \right)\mathclose{} + \frac{1}{2} \cdot h^{2} \cdot F_{f}\mathopen{}\left( \rootedtree[o[.]] \right)\mathclose{} + \frac{1}{2} \cdot h^{2} \cdot F_{f}\mathopen{}\left( \rootedtree[.[o]] \right)\mathclose{} + \frac{1}{2} \cdot h^{2} \cdot F_{f}\mathopen{}\left( \rootedtree[o[o]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[.[.]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[.[o]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[.[o]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[o[o]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[o[o]]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[.][.]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[o][.]] \right)\mathclose{} + \frac{1}{4} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[o][.]] \right)\mathclose{} + \frac{1}{8} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[.[o][o]] \right)\mathclose{} + \frac{1}{8} \cdot h^{3} \cdot F_{f}\mathopen{}\left( \rootedtree[o[o][o]] \right)\mathclose{}\]
When compiled using the preamble code shown in the docstring of RootedTrees.latexify, the output looks as follows.
References
- Ernst Hairer, Gerhard Wanner, Christian Lubich. Geometric numerical integration. Springer, 2002. DOI: 10.1007/3-540-30666-8
B-series for Rosenbrock methods
BSeries.jl and RootedTrees.jl also support Rosenbrock (Rosenbrock-Wanner, ROW) methods via the wrapper RosenbrockMethod. For example, a classical ROW method of Kaps and Rentrop (1979) can be parameterized as follows.
using BSeries
γ = [0.395 0 0 0;
-0.767672395484 0.395 0 0;
-0.851675323742 0.522967289188 0.395 0;
0.288463109545 0.880214273381e-1 -0.337389840627 0.395]
A = [0 0 0 0;
0.438 0 0 0;
0.796920457938 0.730795420615e-1 0 0;
0.796920457938 0.730795420615e-1 0 0]
b = [0.199293275701, 0.482645235674, 0.680614886256e-1, 0.25]
ros = RosenbrockMethod(γ, A, b)RosenbrockMethod{Float64} with
γ: 4×4 Matrix{Float64}:
0.395 0.0 0.0 0.0
-0.767672 0.395 0.0 0.0
-0.851675 0.522967 0.395 0.0
0.288463 0.0880214 -0.33739 0.395
A: 4×4 Matrix{Float64}:
0.0 0.0 0.0 0.0
0.438 0.0 0.0 0.0
0.79692 0.0730795 0.0 0.0
0.79692 0.0730795 0.0 0.0
b: 4-element Vector{Float64}:
0.199293275701
0.482645235674
0.0680614886256
0.25
c: 4-element Vector{Float64}:
0.0
0.438
0.8699999999995001
0.8699999999995001
We can create the B-series as usual, truncated to order 5.
series = bseries(ros, 5)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Float64} with 18 entries:
RootedTree{Int64}: Int64[] => 1.0
RootedTree{Int64}: [1] => 1.0
RootedTree{Int64}: [1, 2] => 0.5
RootedTree{Int64}: [1, 2, 3] => 0.166667
RootedTree{Int64}: [1, 2, 2] => 0.333333
RootedTree{Int64}: [1, 2, 3, 4] => 0.0416667
RootedTree{Int64}: [1, 2, 3, 3] => 0.0833333
RootedTree{Int64}: [1, 2, 3, 2] => 0.125
RootedTree{Int64}: [1, 2, 2, 2] => 0.25
RootedTree{Int64}: [1, 2, 3, 4, 5] => 0.00872428
RootedTree{Int64}: [1, 2, 3, 4, 4] => 0.00418004
RootedTree{Int64}: [1, 2, 3, 4, 3] => 0.0320513
RootedTree{Int64}: [1, 2, 3, 4, 2] => 0.0467417
RootedTree{Int64}: [1, 2, 3, 3, 3] => 0.05
RootedTree{Int64}: [1, 2, 3, 3, 2] => 0.00387949
RootedTree{Int64}: [1, 2, 3, 2, 3] => 0.0469263
RootedTree{Int64}: [1, 2, 3, 2, 2] => 0.09295
RootedTree{Int64}: [1, 2, 2, 2, 2] => 0.19998Again, we can check the order of accuracy by comparing the coefficients to the exact solution.
series - ExactSolution(series)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Float64} with 18 entries:
RootedTree{Int64}: Int64[] => 0.0
RootedTree{Int64}: [1] => 5.99965e-13
RootedTree{Int64}: [1, 2] => -1.52767e-13
RootedTree{Int64}: [1, 2, 3] => -2.12941e-13
RootedTree{Int64}: [1, 2, 2] => -2.50522e-13
RootedTree{Int64}: [1, 2, 3, 4] => -1.37459e-13
RootedTree{Int64}: [1, 2, 3, 3] => -8.44186e-14
RootedTree{Int64}: [1, 2, 3, 2] => -8.9137e-14
RootedTree{Int64}: [1, 2, 2, 2] => -3.5999e-13
RootedTree{Int64}: [1, 2, 3, 4, 5] => 0.000390949
RootedTree{Int64}: [1, 2, 3, 4, 4] => -0.0124866
RootedTree{Int64}: [1, 2, 3, 4, 3] => 0.00705128
RootedTree{Int64}: [1, 2, 3, 4, 2] => 0.0134083
RootedTree{Int64}: [1, 2, 3, 3, 3] => -9.4473e-14
RootedTree{Int64}: [1, 2, 3, 3, 2] => -0.0627872
RootedTree{Int64}: [1, 2, 3, 2, 3] => -0.00307372
RootedTree{Int64}: [1, 2, 3, 2, 2] => -0.00705
RootedTree{Int64}: [1, 2, 2, 2, 2] => -2.0e-5order_of_accuracy(series)4References
- Peter Kaps and Peter Rentrop. "Generalized Runge-Kutta methods of order four with stepsize control for stiff ordinary differential equations." Numerische Mathematik 33, no. 1 (1979): 55-68. DOI: 10.1007/BF01396495
B-series for the average vector field method
Consider the autonomous ODE
\[u'(t) = f\bigl( u(t) \bigr).\]
The average vector field method
\[u^{n+1} = u^{n} + \Delta t \int_0^1 f\bigl(\xi u^{n+1} + (1 - \xi) u^{n}\bigr) \mathrm{d} \xi\]
introduced by McLachlan, Quispel, and Robidoux (1999) is a second-order accurate method. Quispel and McLaren (2008) discovered that it is indeed a B-series method. Its coefficients are given explicitly as
\[\begin{align*} b(.) &= 1, \\ b([t_1, ..., t_n]) &= b(t_1)...b(t_n) / (n + 1). \end{align*}\]
by Celledoni, McLachlan, McLaren, Owren, Quispel, and Wright (2009). We can implement this up to order 5 in BSeries.jl as follows.
using BSeries
series = bseries(5) do t, series
if order(t) in (0, 1)
return 1 // 1
else
v = 1 // 1
n = 0
for subtree in SubtreeIterator(t)
v *= series[subtree]
n += 1
end
return v / (n + 1)
end
endTruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
RootedTree{Int64}: Int64[] => 1//1
RootedTree{Int64}: [1] => 1//1
RootedTree{Int64}: [1, 2] => 1//2
RootedTree{Int64}: [1, 2, 3] => 1//4
RootedTree{Int64}: [1, 2, 2] => 1//3
RootedTree{Int64}: [1, 2, 3, 4] => 1//8
RootedTree{Int64}: [1, 2, 3, 3] => 1//6
RootedTree{Int64}: [1, 2, 3, 2] => 1//6
RootedTree{Int64}: [1, 2, 2, 2] => 1//4
RootedTree{Int64}: [1, 2, 3, 4, 5] => 1//16
RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//12
RootedTree{Int64}: [1, 2, 3, 4, 3] => 1//12
RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//12
RootedTree{Int64}: [1, 2, 3, 3, 3] => 1//8
RootedTree{Int64}: [1, 2, 3, 3, 2] => 1//9
RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//12
RootedTree{Int64}: [1, 2, 3, 2, 2] => 1//8
RootedTree{Int64}: [1, 2, 2, 2, 2] => 1//5BSeries.jl also offers a convenience constructor using the type AverageVectorFieldMethod as follows.
series == bseries(AverageVectorFieldMethod(), 5)trueWe can check that this method is second-order accurate by comparing it to the B-series of the exact solution, truncated at the same order.
series - ExactSolution(series)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 18 entries:
RootedTree{Int64}: Int64[] => 0//1
RootedTree{Int64}: [1] => 0//1
RootedTree{Int64}: [1, 2] => 0//1
RootedTree{Int64}: [1, 2, 3] => 1//12
RootedTree{Int64}: [1, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4] => 1//12
RootedTree{Int64}: [1, 2, 3, 3] => 1//12
RootedTree{Int64}: [1, 2, 3, 2] => 1//24
RootedTree{Int64}: [1, 2, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 5] => 13//240
RootedTree{Int64}: [1, 2, 3, 4, 4] => 1//15
RootedTree{Int64}: [1, 2, 3, 4, 3] => 7//120
RootedTree{Int64}: [1, 2, 3, 4, 2] => 1//20
RootedTree{Int64}: [1, 2, 3, 3, 3] => 3//40
RootedTree{Int64}: [1, 2, 3, 3, 2] => 2//45
RootedTree{Int64}: [1, 2, 3, 2, 3] => 1//30
RootedTree{Int64}: [1, 2, 3, 2, 2] => 1//40
RootedTree{Int64}: [1, 2, 2, 2, 2] => 0//1order_of_accuracy(series)2References
- Robert I. McLachlan, G. Reinout W. Quispel, and Nicolas Robidoux. "Geometric integration using discrete gradients." Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 357, no. 1754 (1999): 1021-1045. DOI: 10.1098/rsta.1999.0363
- G. Reinout W. Quispel, and David Ian McLaren. "A new class of energy-preserving numerical integration methods." Journal of Physics A: Mathematical and Theoretical 41, no. 4 (2008): 045206. DOI: 10.1088/1751-8113/41/4/045206
- Elena Celledoni, Robert I. McLachlan, David I. McLaren, Brynjulf Owren, G. Reinout W. Quispel, and William M. Wright. "Energy-preserving Runge-Kutta methods." ESAIM: Mathematical Modelling and Numerical Analysis 43, no. 4 (2009): 645-649. DOI: 10.1051/m2an/2009020
Composing B-series
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//1
RootedTree{Int64}: [1] => 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//1
⋮ => ⋮Note that this method is fourth-order accurate:
order_of_accuracy(series_a)4Next, 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//1
RootedTree{Int64}: [1] => 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//1
⋮ => ⋮order_of_accuracy(series_b)3Note 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//1
RootedTree{Int64}: [1] => 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//1
⋮ => ⋮order_of_accuracy(series_c)3Finally, 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//1
RootedTree{Int64}: [1] => 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)5series - ExactSolution(series)TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 38 entries:
RootedTree{Int64}: Int64[] => 0//1
RootedTree{Int64}: [1] => 0//1
RootedTree{Int64}: [1, 2] => 0//1
RootedTree{Int64}: [1, 2, 3] => 0//1
RootedTree{Int64}: [1, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4] => 0//1
RootedTree{Int64}: [1, 2, 3, 3] => 0//1
RootedTree{Int64}: [1, 2, 3, 2] => 0//1
RootedTree{Int64}: [1, 2, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 5] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 4] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 3] => 0//1
RootedTree{Int64}: [1, 2, 3, 4, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 3, 3] => 0//1
RootedTree{Int64}: [1, 2, 3, 3, 2] => 0//1
RootedTree{Int64}: [1, 2, 3, 2, 3] => 0//1
RootedTree{Int64}: [1, 2, 3, 2, 2] => 0//1
RootedTree{Int64}: [1, 2, 2, 2, 2] => 0//1
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