Rosenbrock methods

Consider an ordinary differential equation (ODE) of the form

\[u'(t) = f(t, u(t)).\]

A Rosenbrock (or Rosenbrock-Wanner, ROW) method with $s$ stages is given by its coefficients

\[\gamma = (\gamma_{i,j})_{i,j} \in \mathbb{R}^{s \times s}, \quad A = (a_{i,j})_{i,j} \in \mathbb{R}^{s \times s}, \quad b = (b_i)_i \in \mathbb{R}^{s}, \quad c = (c_i)_i \in \mathbb{R}^{s}.\]

Usually, the consistency condition

\[\forall i\colon \quad c_i = \sum_j a_{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} k^i &= \Delta t f\bigl(u^n + \sum_j a_{i,j} k^j \bigr) + \Delta t J \sum_j \gamma_{ij} k_j, \\ u^{n+1} &= u^n + \sum_i b_{i} k^i. \end{aligned}\]

In RootedTrees.jl, ROW methods are represented as RosenbrockMethods.

Order conditions

The order conditions of ROW methods can be derived using rooted trees. In RootedTrees.jl, this functionality is again implemented in residual_order_condition. Thus, a RosenbrockMethod is of order $p$ if the residual_order_condition vanishes for all rooted trees with order up to $p$.

For example, the method GRK4A of Kaps and Rentrop (1979) can be written as follows.

using RootedTrees

γ = [0.395 0 0 0;
     -0.767672395484 0.395 0 0;
     -0.851675323742 0.522967289188 0.395 0;
     0.288463109545 0.880214273381e-1 -.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

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

using Test

@testset "GRK4A, order 4" begin
  for o in 0:4
    for t in RootedTreeIterator(o)
      val = residual_order_condition(t, ros)
      @test abs(val) < 3000 * eps()
    end
  end
end

Test Summary:  | Pass  Total  Time
GRK4A, order 4 |    9      9  0.1s

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

@testset "GRK4A, not order 5" begin
  s = 0.0
  for t in RootedTreeIterator(5)
    s += abs(residual_order_condition(t, ros))
  end
  @test s > 0.06
end

Test Summary:      | Pass  Total  Time
GRK4A, not order 5 |    1      1  0.0s