Rosenbrock methods (#73)

* support Rosenbrock methods

* set version to 0.1.31

* mention Rosenbrock methods in tutorials

* format

* fix typo

* linking from tutorials to stuff defined in RootedTrees does not work

Project.toml

@@ -1,7 +1,7 @@
 name = "BSeries"
 uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
 authors = ["Hendrik Ranocha <[email protected]> and contributors"]
-version = "0.1.30"
+version = "0.1.31"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
@@ -18,5 +18,5 @@ OrderedCollections = "1"
 Polynomials = "2.0.23, 3"
 Reexport = "1"
 Requires = "1"
-RootedTrees = "2.12"
+RootedTrees = "2.13"
 julia = "1.6"

docs/src/tutorials/bseries_creation.md

@@ -91,6 +91,56 @@ order_of_accuracy(series)
   Springer, 2002.
 
 
+## B-series for Rosenbrock methods
+
+[BSeries.jl](https://github.com/ranocha/BSeries.jl) and
+[RootedTrees.jl](https://github.com/SciML/RootedTrees.jl) also support
+Rosenbrock (Rosenbrock-Wanner, ROW) methods via the wrapper
+`RosebrockMethod`. For example, a classical ROW method of
+Kaps and Rentrop (1979) can be parameterized as follows.
+
+```@example ex:ROW
+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)
+```
+
+We can create the B-series as usual, truncated to order 5.
+
+```@example ex:ROW
+series = bseries(ros, 5)
+```
+
+Again, we can check the order of accuracy by comparing the coefficients to
+the exact solution.
+
+```@example ex:ROW
+series - ExactSolution(series)
+```
+
+```@example ex:ROW
+order_of_accuracy(series)
+```
+
+
+### References
+
+- 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](https://doi.org/10.1007/BF01396495)
+
+
 ## [B-series for the average vector field method](@id tutorial-bseries-creation-AVF)
 
 Consider the autonomous ODE

src/BSeries.jl

@@ -501,6 +501,44 @@ end
 # TODO: bseries(ark::AdditiveRungeKuttaMethod)
 # should create a lazy version, optionally a memoized one
 
+"""
+    bseries(ros::RosenbrockMethod, order)
+
+Compute the B-series of the Rosenbrock method `ros` up to a prescribed
+integer `order`.
+
+!!! note "Normalization by elementary differentials"
+    The coefficients of the B-series returned by this method need to be
+    multiplied by a power of the time step divided by the `symmetry` of the
+    rooted tree and multiplied by the corresponding elementary differential
+    of the input vector field ``f``.
+    See also [`evaluate`](@ref).
+"""
+function bseries(ros::RosenbrockMethod, order)
+    V_tmp = eltype(ros)
+    if V_tmp <: Integer
+        # If people use integer coefficients, they will likely want to have results
+        # as exact as possible. However, general terms are not integers. Thus, we
+        # use rationals instead.
+        V = Rational{V_tmp}
+    else
+        V = V_tmp
+    end
+    series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()
+
+    series[rootedtree(Int[])] = one(V)
+    for o in 1:order
+        for t in RootedTreeIterator(o)
+            series[copy(t)] = elementary_weight(t, ros)
+        end
+    end
+
+    return series
+end
+
+# TODO: bseries(ros::RosenbrockMethod)
+# should create a lazy version, optionally a memoized one
+
 # TODO: Documentation, Base.show, export etc.
 """
     MultirateInfinitesimalSplitMethod(A, D, G, c)

test/runtests.jl

@@ -1432,6 +1432,32 @@ using Aqua: Aqua
         end
     end # @testset "additive Runge-Kutta methods interface"
 
+    @testset "Rosenbrock methods interface" 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)
+
+        # fourth-order accurate
+        series_integrator = @inferred bseries(ros, 5)
+        @test @inferred(order_of_accuracy(series_integrator)) == 4
+
+        # not fifth-order accurate
+        series_exact = @inferred ExactSolution(series_integrator)
+        @test mapreduce(isapprox, &, values(series_integrator), values(series_exact)) ==
+              false
+    end # @testset "Rosenbrock methods interface"
+
     @testset "multirate infinitesimal split methods interface" begin
         # NOTE: This is still considered experimental and might change at any time!