Merge pull request #12 from ranocha/hr/compose

compose

Project.toml

@@ -13,7 +13,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 [compat]
 OrderedCollections = "1"
 Reexport = "1"
-RootedTrees = "2"
+RootedTrees = "2.1"
 Symbolics = "3"
 julia = "1"
 

src/BSeries.jl

@@ -11,7 +11,7 @@ using RootedTrees: RootedTree
 using Symbolics: Differential, expand_derivatives
 
 
-export bseries, substitute
+export bseries, substitute, compose
 
 export modified_equation, modifying_integrator
 
@@ -31,6 +31,14 @@ end
 
 Compute the coefficient correspoding to the tree `t` of the B-series that is
 formed by substituting the B-series `b` into the B-series `a`.
+
+# References
+
+Section 3.2 of
+- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
+  Algebraic Structures of B-series
+  Foundations of Computational Mathematics
+  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
 """
 function substitute(b, a, t::RootedTree)
   result = zero(first(values(a)) * first(values(b)))
@@ -43,6 +51,35 @@ function substitute(b, a, t::RootedTree)
 end
 
 
+"""
+    compose(b, a, t::RootedTree)
+
+Compute the coefficient correspoding to the tree `t` of the B-series that is
+formed by composing the B-series `a` with the B-series `b`.
+
+# References
+
+Section 3.1 of
+- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010)
+  Algebraic Structures of B-series
+  Foundations of Computational Mathematics
+  [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1)
+"""
+function compose(b, a, t::RootedTree)
+  result = zero(first(values(a)) * first(values(b)))
+
+  for (forest, subtree) in SplittingIterator(t)
+    if isempty(forest)
+      result += a[subtree]
+    else
+      result += reduce(*, b[tree] for tree in forest) * a[subtree]
+    end
+  end
+
+  return result
+end
+
+
 
 """
     bseries(A::AbstractMatrix, b::AbstractVector, c::AbstractVector, order)

test/runtests.jl

@@ -14,7 +14,8 @@ using Symbolics: Symbolics, @variables, Num, expand
   @test terms == [1//1, 1//2, 1//6, 1//3]
 end
 
-@testset "subsitution" begin
+
+@testset "substitute" begin
   @variables a1 a2 a31 a32
   a = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Num}(
         rootedtree(Int[])     => 1,
@@ -31,16 +32,68 @@ end
         rootedtree([1, 2, 2]) => b31,
         rootedtree([1, 2, 3]) => b32)
 
+  # See equation (6) of
+  # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
+  #   Numerical integrators based on modified differential equations
+  #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)
+
   t = rootedtree([1])
   @inferred substitute(b, a, t)
   @test isequal(substitute(b, a, t), a1 * b1)
 
+  t = rootedtree([1, 2])
+  @inferred substitute(b, a, t)
+  @test isequal(substitute(b, a, t), a1 * b2 + a2 * b1^2)
+
+  t = rootedtree([1, 2, 2])
+  @inferred substitute(b, a, t)
+  @test isequal(substitute(b, a, t), a1 * b31 + 2 * a2 * b1 * b2 + a31 * b1^3)
+
   t = rootedtree([1, 2, 3])
   @inferred substitute(b, a, t)
   @test isequal(substitute(b, a, t), a1 * b32 + 2 * a2 * b1 * b2 + a32 * b1^3)
 end
 
 
+@testset "compose" begin
+  @variables a0 a1 a2 a31 a32
+  a = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Num}(
+        rootedtree(Int[])     => a0,
+        rootedtree([1])       => a1,
+        rootedtree([1, 2])    => a2,
+        rootedtree([1, 2, 2]) => a31,
+        rootedtree([1, 2, 3]) => a32)
+
+  @variables b1 b2 b31 b32
+  b = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Num}(
+        rootedtree(Int[])     => 0,
+        rootedtree([1])       => b1,
+        rootedtree([1, 2])    => b2,
+        rootedtree([1, 2, 2]) => b31,
+        rootedtree([1, 2, 3]) => b32)
+
+  # See equation (5) of
+  # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
+  #   Numerical integrators based on modified differential equations
+  #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)
+  t = rootedtree([1])
+  @inferred compose(b, a, t)
+  @test isequal(compose(b, a, t), a0 * b1 + a1)
+
+  t = rootedtree([1, 2])
+  @inferred compose(b, a, t)
+  @test isequal(compose(b, a, t), a0 * b2 + a1 * b1 + a2)
+
+  t = rootedtree([1, 2, 2])
+  @inferred compose(b, a, t)
+  @test isequal(compose(b, a, t), a0 * b31 + a1 * b1^2 + 2 * a2 * b1 + a31)
+
+  t = rootedtree([1, 2, 3])
+  @inferred compose(b, a, t)
+  @test isequal(compose(b, a, t), a0 * b32 + a1 * b2 + a2 * b1 + a32)
+end
+
+
 @testset "modified_equation" begin
   t1  = rootedtree([1])
   t2  = rootedtree([1, 2])