Runge-Kutta order conditions

In this tutorial, we generate order conditions for a generic explicit Runge-Kutta method.

First, we create symbolic coefficient arrays with the appropriate structure. There are several symbolic packages you can use. Here, we will use SymPy.jl.

using BSeries, SymPy, Latexify

s = 4  # Stages
p = 4  # Desired order of accuracy

A = Array{Sym,2}(undef, s, s)
b = Array{Sym,1}(undef, s)
c = Array{Sym,1}(undef, s);
for i in 1:s
    b[i] = symbols("b$i", real = true)
    for j in 1:i-1
        A[i, j] = symbols("a$i$j", real = true)
    end
    for j in i:s
        A[i, j] = 0
    end
end

for i in 1:s
    c[i] = 0
    for j in 1:i-1
        c[i] += A[i, j]
    end
end

Next we generate the B-series for the RK method and the exact solution, and take their difference.

rk3 = bseries(A, b, c, p)
exact = ExactSolution(rk3)
error = rk3 - exact

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPy.Sym} with 9 entries:
  RootedTree{Int64}: Int64[]      => 0
  RootedTree{Int64}: [1]          => b1 + b2 + b3 + b4 - 1
  RootedTree{Int64}: [1, 2]       => a21*b2 + b3*(a31 + a32) + b4*(a41 + a42 + …
  RootedTree{Int64}: [1, 2, 3]    => a21*a32*b3 + b4*(a21*a42 + a43*(a31 + a32)…
  RootedTree{Int64}: [1, 2, 2]    => a21^2*b2 + b3*(a31 + a32)^2 + b4*(a41 + a4…
  RootedTree{Int64}: [1, 2, 3, 4] => a21*a32*a43*b4 - 1/24
  RootedTree{Int64}: [1, 2, 3, 3] => a21^2*a32*b3 + b4*(a21^2*a42 + a43*(a31 + …
  RootedTree{Int64}: [1, 2, 3, 2] => a21*a32*b3*(a31 + a32) + b4*(a21*a42 + a43…
  RootedTree{Int64}: [1, 2, 2, 2] => a21^3*b2 + b3*(a31 + a32)^3 + b4*(a41 + a4…

The output above is truncated. We can see the full expressions as follows:

for (tree, order_condition) in error
    println(order_condition)
end

0
b1 + b2 + b3 + b4 - 1
a21*b2 + b3*(a31 + a32) + b4*(a41 + a42 + a43) - 1/2
a21*a32*b3 + b4*(a21*a42 + a43*(a31 + a32)) - 1/6
a21^2*b2 + b3*(a31 + a32)^2 + b4*(a41 + a42 + a43)^2 - 1/3
a21*a32*a43*b4 - 1/24
a21^2*a32*b3 + b4*(a21^2*a42 + a43*(a31 + a32)^2) - 1/12
a21*a32*b3*(a31 + a32) + b4*(a21*a42 + a43*(a31 + a32))*(a41 + a42 + a43) - 1/8
a21^3*b2 + b3*(a31 + a32)^3 + b4*(a41 + a42 + a43)^3 - 1/4