Additive Runge-Kutta methods

Consider an ordinary differential equation (ODE) of the form

\[u'(t) = \sum_\nu^N f^\nu(t, u(t)).\]

An additive Runge-Kutta (ARK) method with $s$ stages is given by its Butcher coefficients

\[A^\nu = (a^\nu_{i,j})_{i,j} \in \mathbb{R}^{s \times s}, \quad b^\nu = (b^\nu_i)_i \in \mathbb{R}^{s}, \quad c^\nu = (c^\nu_i)_i \in \mathbb{R}^{s}.\]

Usually, the consistency condition

\[\forall i\colon \quad c^\nu_i = \sum_j a^\nu_{i,j}\]

is assumed, which reduces all analysis to autonomous problems.

The step from $u^{n}$ to $u^{n+1}$ is given by

\[\begin{aligned} y^i &= u^n + \Delta t \sum_\nu \sum_j a^\nu_{i,j} f^\nu(y^i), \\ u^{n+1} &= u^n + \Delta t \sum_\nu \sum_i b^\nu_{i} f^\nu(y^i), \end{aligned}\]

where $y^i$ are the stage values.

In RootedTrees.jl, ARK methods are represented as AdditiveRungeKuttaMethods.

Order conditions

The order conditions of RK methods can be derived using colored rooted trees. In RootedTrees.jl, this functionality is implemented in residual_order_condition. Thus, an AdditiveRungeKuttaMethod is of order $p$ if the residual_order_condition vanishes for all colored rooted trees with order up to $p$ and $N$ colors. The most important case is $N = 2$, i.e., BicoloredRootedTrees as special case of ColoredRootedTrees.

For example, the classical Störmer-Verlet method can be written as follows, see Table II.2.1 of Hairer, Lubich, Wanner (2002) Geometric numerical integration.

using RootedTrees

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     0
 1//2  1//2
b: 2-element Vector{Rational{Int64}}:
 1//2
 1//2
c: 2-element Vector{Rational{Int64}}:
 0
 1
2. RungeKuttaMethod{Rational{Int64}} with
A: 2×2 Matrix{Rational{Int64}}:
 1//2  0
 1//2  0
b: 2-element Vector{Rational{Int64}}:
 1//2
 1//2
c: 2-element Vector{Rational{Int64}}:
 1//2
 1//2

To verify that this method is at least second-order accurate, we can check the residual_order_conditions up to this order.

using Test

@testset "Störmer-Verlet, order 2" begin
  for o in 1:2
    for t in BicoloredRootedTreeIterator(o)
      @test iszero(residual_order_condition(t, ark))
    end
  end
end

Test Summary:           | Pass  Total  Time
Störmer-Verlet, order 2 |    6      6  0.1s

To verify that this method does not satisfy any of the order conditions for an order of accuracy of three, we can use the following code.

using Test

@testset "Störmer-Verlet, not order 3" begin
  for t in BicoloredRootedTreeIterator(3)
    @test !iszero(residual_order_condition(t, ark))
  end
end

Test Summary:               | Pass  Total  Time
Störmer-Verlet, not order 3 |   14     14  0.1s