fix evaluate for B-series of numerical integrators (#214)

* fix bug from issue 213

* add test

Project.toml

@@ -1,7 +1,7 @@
 name = "BSeries"
 uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
 authors = ["Hendrik Ranocha <[email protected]> and contributors"]
-version = "0.1.60"
+version = "0.1.61"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/BSeries.jl

@@ -1537,7 +1537,16 @@ function elementary_differentials(f::AbstractVector, u, order)
     differentials = OrderedDict{RootedTree{Int, Vector{Int}}, typeof(f)}()
 
     t = rootedtree(Int[])
-    differentials[t] = one.(f)
+    if typeof(u) == typeof(f)
+        differentials[t] = u
+    else
+        # `u` may be a tuple of symbolic variables
+        # Since we determine the type of the values of `differentials` based on
+        # the type of `f`, we need to create a vector of the same type as `f`.
+        val = zero(f)
+        val .= u
+        differentials[t] = val
+    end
 
     t = rootedtree([1])
     differentials[t] = f

test/runtests.jl

@@ -640,7 +640,19 @@ using Aqua: Aqua
         end
     end # @testset "modified_equation"
 
-    @testset "elementary differentials" begin
+    @testset "elementary differentials and evaluate" begin
+        @testset "Rooted trees" begin
+            @testset "Issue #213" begin
+                y, h = SymPyPythonCall.symbols("y, h")
+                u = [y]
+                f = [-y]
+                series = @inferred bseries(AverageVectorFieldMethod(), 3)
+                expansion = first(evaluate(f, u, h, series))
+                reference = y - h * y + h^2 * y / 2 - h^3 * y / 4
+                @test expansion == reference
+            end
+        end
+
         @testset "Bicolored trees" begin
             @testset "Lotka-Volterra" begin
                 # Verified with Mathematica