rosenbrock methods (#89)

* RosenbrockMethod

* RosenbrockMethod tests

* add tutorial

* format

* test show and add new tutorial to docs/make.jl

* set version to 2.13.0

* mention ROW in RosenbrockMethod docstring

* fix typos in new tutorial

Author: ranocha

Project.toml

@@ -1,7 +1,7 @@
 name = "RootedTrees"
 uuid = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
 authors = ["Hendrik Ranocha <mail@ranocha.de> and contributors"]
-version = "2.12.2"
+version = "2.13.0"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"

docs/make.jl

@@ -53,6 +53,7 @@ makedocs(modules = [RootedTrees],
              "Tutorials" => [
                  "tutorials/RK_order_conditions.md",
                  "tutorials/ARK_order_conditions.md",
+                 "tutorials/Rosenbrock_order_conditions.md",
              ],
              # "Benchmarks" => "benchmarks.md",
              "API reference" => "api_reference.md",

docs/src/tutorials/ARK_order_conditions.md

@@ -37,7 +37,7 @@ ARK methods are represented as
 The order conditions of RK methods can be derived using colored rooted trees.
 In [RootedTrees.jl](https://github.com/SciML/RootedTrees.jl), this
 functionality is implemented in [`residual_order_condition`](@ref).
-Thus, a [`RungeKuttaMethod`](@ref) is of order ``p`` if the
+Thus, an [`AdditiveRungeKuttaMethod`](@ref) is of order ``p`` if the
 [`residual_order_condition`](@ref) vanishes for all colored rooted trees
 with [`order`](@ref) up to ``p`` and ``N`` colors. The most important case
 is ``N = 2``, i.e., [`BicoloredRootedTree`](@ref)s as special case of

docs/src/tutorials/Rosenbrock_order_conditions.md

@@ -0,0 +1,93 @@
+# Rosenbrock methods
+
+Consider an ordinary differential equation (ODE) of the form
+```math
+u'(t) = f(t, u(t)).
+```
+
+A Rosenbrock (or Rosenbrock-Wanner, ROW) method with ``s`` stages is given
+by its coefficients
+```math
+\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
+```math
+\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
+```math
+\begin{aligned}
+  k^i &= \Delta t f\bigbl(u^n + \sum_j a_{i,j} k^j \bibgr) + \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](https://github.com/SciML/RootedTrees.jl),
+ROW methods are represented as
+[`RosenbrockMethod`](@ref)s.
+
+
+## Order conditions
+
+The order conditions of ROW methods can be derived using rooted trees.
+In [RootedTrees.jl](https://github.com/SciML/RootedTrees.jl), this
+functionality is again implemented in [`residual_order_condition`](@ref).
+Thus, a [`RosenbrockMethod`](@ref) is of order ``p`` if the
+[`residual_order_condition`](@ref) vanishes for all rooted trees
+with [`order`](@ref) up to ``p``.
+
+For example, the method GRK4A of
+[Kaps and Rentrop (1979)](https://doi.org/10.1007/BF01396495) can be
+written as follows.
+
+```@example GRK4A
+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)
+```
+
+To verify that this method is at least fourth-order accurate, we can
+check the [`residual_order_condition`](@ref)s up to this order.
+
+```@example GRK4A
+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
+nothing # hide
+```
+
+To verify that this method does not satisfy the order conditions
+for an order of accuracy of five, we can use the following code.
+
+```@example GRK4A
+
+@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
+nothing # hide
+```

src/RootedTrees.jl

@@ -27,7 +27,7 @@ export partition_forest, PartitionForestIterator,
 
 export all_splittings, SplittingIterator
 
-export RungeKuttaMethod, AdditiveRungeKuttaMethod
+export RungeKuttaMethod, AdditiveRungeKuttaMethod, RosenbrockMethod
 
 abstract type AbstractRootedTree end
 

src/time_integration_methods.jl

@@ -287,3 +287,131 @@ function residual_order_condition(t::ColoredRootedTree, ark::AdditiveRungeKuttaM
 
     (ew - one(T) / γ(t)) / σ(t)
 end
+
+"""
+    RosenbrockMethod(γ, A, b, c=vec(sum(A, dims=2)))
+
+Represent a Rosenbrock (or Rosenbrock-Wanner, ROW) method with
+coefficients `γ`, `A`, `b`, and `c`.
+If `c` is not provided, the usual "row sum" requirement of consistency with
+autonomous problems is applied.
+
+# Reference
+
+- Ernst Hairer, Gerhard Wanner.
+  Solving ordinary differential equations II: Stiff and differential-algebraic problems.
+  Springer, 2010.
+  Section IV.7
+"""
+struct RosenbrockMethod{T, MatT <: AbstractMatrix{T}, VecT <: AbstractVector{T}} <:
+       AbstractTimeIntegrationMethod
+    γ::MatT
+    A::MatT
+    b::VecT
+    c::VecT
+end
+
+function RosenbrockMethod(γ::AbstractMatrix, A::AbstractMatrix,
+                          b::AbstractVector,
+                          c::AbstractVector = vec(sum(A, dims = 2)))
+    T = promote_type(eltype(γ), eltype(A), eltype(b), eltype(c))
+    _γ = T.(γ)
+    _A = T.(A)
+    _b = T.(b)
+    _c = T.(c)
+    return RosenbrockMethod(_γ, _A, _b, _c)
+end
+
+Base.eltype(ros::RosenbrockMethod{T}) where {T} = T
+
+function Base.show(io::IO, ros::RosenbrockMethod)
+    print(io, "RosenbrockMethod{", eltype(ros), "}")
+    if get(io, :compact, false)
+        print(io, "(")
+        show(io, ros.γ)
+        print(io, ", ")
+        show(io, ros.A)
+        print(io, ", ")
+        show(io, ros.b)
+        print(io, ", ")
+        show(io, ros.c)
+        print(io, ")")
+    else
+        print(io, " with\nγ: ")
+        show(io, MIME"text/plain"(), ros.γ)
+        print(io, "\nA: ")
+        show(io, MIME"text/plain"(), ros.A)
+        print(io, "\nb: ")
+        show(io, MIME"text/plain"(), ros.b)
+        print(io, "\nc: ")
+        show(io, MIME"text/plain"(), ros.c)
+        print(io, "\n")
+    end
+end
+
+"""
+    elementary_weight(t::RootedTree, ros::RosenbrockMethod)
+
+Compute the elementary weight Φ(`t`) of the [`RosenbrockMethod`](@ref) `ros`
+for a rooted tree `t``.
+"""
+function elementary_weight(t::RootedTree, ros::RosenbrockMethod)
+    dot(ros.b, derivative_weight(t, ros))
+end
+
+"""
+    derivative_weight(t::RootedTree, ros::RosenbrockMethod)
+
+Compute the derivative weight (ΦᵢD)(`t`) of the [`RosenbrockMethod`](@ref) `ros`
+for the rooted tree `t`.
+"""
+function derivative_weight(t::RootedTree, ros::RosenbrockMethod)
+    γ = ros.γ
+    A = ros.A
+    c = ros.c
+
+    # Initialize `result` with the identity element of pointwise multiplication `.*`
+    result = zero(c) .+ one(eltype(c))
+
+    # Count the number of subtrees to decide which matrix to use for multiplications
+    num_subtrees = 0
+    for subtree in SubtreeIterator(t)
+        num_subtrees += 1
+    end
+    if num_subtrees == 1
+        matrix = A + γ
+    else
+        matrix = A
+    end
+
+    # Iterate over all subtrees and update the `result` using recursion
+    for subtree in SubtreeIterator(t)
+        tmp = matrix * derivative_weight(subtree, ros)
+        result = result .* tmp
+    end
+
+    return result
+end
+
+"""
+    residual_order_condition(t::RootedTree, ros::RosenbrockMethod)
+
+The residual of the order condition
+  `(Φ(t) - 1/γ(t)) / σ(t)`
+with [`elementary_weight`](@ref) `Φ(t)`, [`density`](@ref) `γ(t)`, and
+[`symmetry`](@ref) `σ(t)` of the [`RosenbrockMethod`](@ref) `ros`
+for the rooted tree `t`.
+
+# Reference
+
+- Ernst Hairer, Gerhard Wanner.
+  Solving ordinary differential equations II: Stiff and differential-algebraic problems.
+  Springer, 2010.
+  Section IV.7
+"""
+function residual_order_condition(t::RootedTree, ros::RosenbrockMethod)
+    ew = elementary_weight(t, ros)
+    T = typeof(ew)
+
+    (ew - one(T) / γ(t)) / σ(t)
+end

test/runtests.jl

@@ -1275,6 +1275,62 @@ using Aqua: Aqua
                 @test abs(residual_order_condition(t, ark)) > 10 * eps()
             end
         end
+
+        @testset "RosenbrockMethod, original Rosenbrock" begin
+            γ = [1-sqrt(2) / 2 0;
+                 0 1-sqrt(2) / 2]
+            A = [0 0;
+                 (sqrt(2) - 1)/2 0]
+            b = [0, 1]
+            ros = @inferred RosenbrockMethod(γ, A, b)
+
+            # second-order accurate
+            @test abs(@inferred(residual_order_condition(rootedtree(Int[]), ros))) <
+                  10 * eps()
+            @test abs(@inferred(residual_order_condition(rootedtree(Int[1]), ros))) <
+                  10 * eps()
+            @test abs(@inferred(residual_order_condition(rootedtree(Int[1, 2]), ros))) <
+                  10 * eps()
+            @test abs(@inferred(residual_order_condition(rootedtree(Int[1, 2, 3]), ros))) >
+                  0.04
+            @test abs(@inferred(residual_order_condition(rootedtree(Int[1, 2, 2]), ros))) >
+                  0.14
+        end
+
+        @testset "RosenbrockMethod, GRK4A (Kaps and Rentrop, 1979)" begin
+            # Kaps, Rentrop (1979)
+            # Generalized Runge-Kutta methods of order four with stepsize control
+            # for stiff ordinary differential equations
+            # https://doi.org/10.1007/BF01396495
+            γ = [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 = @inferred RosenbrockMethod(γ, A, b)
+
+            @test_nowarn show(ros)
+            @test_nowarn show(IOContext(stdout, :compact => true), ros)
+
+            # fourth-order accurate
+            for o in 0:4
+                for t in RootedTreeIterator(o)
+                    val = @inferred residual_order_condition(t, ros)
+                    @test abs(val) < 3000 * eps()
+                end
+            end
+
+            # not fifth-order accurate
+            s = 0.0
+            for t in RootedTreeIterator(5)
+                s += abs(residual_order_condition(t, ros))
+            end
+            @test s > 0.06
+        end
     end # @testset "Order conditions"
 
     @testset "plots" begin