implement order_of_accuracy (#71)

* remove debug code

* implement order_of_accuracy

* set version to v0.1.29

* use order_of_accuracy in tutorials

* format

* iterate over rooted tree iterators in order_of_accuracy and add more tests

Project.toml

@@ -1,7 +1,7 @@
 name = "BSeries"
 uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
 authors = ["Hendrik Ranocha <[email protected]> and contributors"]
-version = "0.1.28"
+version = "0.1.29"
 
 [deps]
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"

docs/src/tutorials/bseries_basics.md

@@ -82,9 +82,20 @@ latexify(coeffs4-coeffs_ex, cdot=false)
 This confirms that the method is of 4th order, since all terms involving
 smaller powers of $h$ vanish exactly.  We don't see the $h^6$ and higher
 order terms since we only generated the truncated B-series up to 5th order.
+We can also obtain the order of accuracy awithout comparing the coefficients
+to the exact solution manually:
 
-For the 2nd-order method, we get:
+```@example bseries-basics
+order_of_accuracy(coeffs4)
+```
+
+For the 2nd-order method, we get
+
+```@example bseries-basics
+order_of_accuracy(coeffs2)
+```
 
+with the following leading error terms:
 
 ```@example bseries-basics
 latexify(coeffs2-coeffs_ex, cdot=false)

docs/src/tutorials/bseries_creation.md

@@ -38,6 +38,10 @@ We can check that the classical Runge-Kutta method is indeed fourth-order accura
 series - ExactSolution(series)
 ```
 
+```@example ex:RK4
+order_of_accuracy(series)
+```
+
 
 ## B-series for additive Runge-Kutta methods
 
@@ -75,6 +79,10 @@ the exact solution.
 series - ExactSolution(series)
 ```
 
+```@example ex:SV
+order_of_accuracy(series)
+```
+
 
 ### References
 
@@ -136,6 +144,10 @@ the B-series of the exact solution, truncated at the same order.
 series - ExactSolution(series)
 ```
 
+```@example ex:AVF
+order_of_accuracy(series)
+```
+
 ### References
 
 - Robert I. McLachlan, G. Reinout W. Quispel, and Nicolas Robidoux.
@@ -182,7 +194,7 @@ series_a = bseries(rk_a, 6)
 Note that this method is fourth-order accurate:
 
 ```@example ex:compose
-series_a - ExactSolution(series_a)
+order_of_accuracy(series_a)
 ```
 
 Next, we set up the starting procedure (method "b" in Butcher's paper):
@@ -198,6 +210,10 @@ rk_b = RungeKuttaMethod(A, b)
 series_b = bseries(rk_b, 6)
 ```
 
+```@example ex:compose
+order_of_accuracy(series_b)
+```
+
 Note that this method is only third-order accurate - as is the finishing
 procedure given by
 
@@ -212,6 +228,10 @@ rk_c = RungeKuttaMethod(A, b)
 series_c = bseries(rk_c, 6)
 ```
 
+```@example ex:compose
+order_of_accuracy(series_c)
+```
+
 Finally, we can compose the three methods to obtain
 
 ```@example ex:compose
@@ -221,6 +241,10 @@ series = compose(series_b, series_a, series_c, normalize_stepsize = true)
 Note that this composition has to be read from left to right. Finally, we check
 that the resulting `series` is indeed fifth-order accurate:
 
+```@example ex:compose
+order_of_accuracy(series)
+```
+
 ```@example ex:compose
 series - ExactSolution(series)
 ```

src/BSeries.jl

@@ -19,6 +19,8 @@ using Latexify: Latexify, LaTeXString
 
 export TruncatedBSeries, ExactSolution
 
+export order_of_accuracy
+
 export bseries, substitute, compose, evaluate
 
 export modified_equation, modifying_integrator
@@ -120,6 +122,8 @@ end
 
 The maximal `order` of a rooted tree with non-vanishing coefficient in the
 truncated B-series `series`.
+
+See also [`order_of_accuracy`](@ref).
 """
 RootedTrees.order(series::TruncatedBSeries) = order(series.coef.keys[end])
 
@@ -187,6 +191,16 @@ for op in (:*, :\)
     end
 end
 
+# Return a function returning an iterator over all rooted trees used as
+# keys for the B-series when given an order of the trees.
+function _iterator_type(::TruncatedBSeries{<:RootedTree})
+    return RootedTreeIterator
+end
+
+function _iterator_type(::TruncatedBSeries{<:BicoloredRootedTree})
+    return BicoloredRootedTreeIterator
+end
+
 """
     ExactSolution{V}()
 
@@ -274,7 +288,6 @@ for op in (:+, :-)
         V = promote_type(valtype(series1), valtype(series2))
         series_result = TruncatedBSeries{T, V}()
         for (key, val1, val2) in zip(series_keys, values(series1), values(series2))
-            @info "op" key val1 val2
             series_result[key] = ($op)(val1, val2)
         end
 
@@ -296,6 +309,42 @@ for op in (:+, :-)
     end
 end
 
+# investigate properties of B-series
+"""
+    order_of_accuracy(series; kwargs...)
+
+Determine the order of accuracy of the B-series `series`. By default, the
+comparison with the coefficients of the exact solution is performed using
+`isequal`. If keyword arguments such as absolute/relative tolerances `atol`/`rtol`
+are given or floating point numbers are used, the comparison is perfomed using
+`isapprox` and the keyword arguments `kwargs...` are forwarded.
+
+See also [`order`](@ref), [`ExactSolution`](@ref).
+"""
+function order_of_accuracy(series::TruncatedBSeries; kwargs...)
+    if isempty(kwargs) && !(valtype(series) <: AbstractFloat)
+        compare = isequal
+    else
+        compare = (a, b) -> isapprox(a, b; kwargs...)
+    end
+
+    exact = ExactSolution{valtype(series)}()
+    iterator = _iterator_type(series)
+
+    for o in 0:order(series)
+        # Iterate over all rooted trees used as keys in `series`
+        # of a given order `o`.
+        for t in iterator(o)
+            if compare(series[t], exact[t]) == false
+                return order(t) - 1
+            end
+        end
+    end
+
+    return order(series)
+end
+
+# construct B-series
 """
     bseries(f::Function, order, iterator_type=RootedTreeIterator)
 

test/runtests.jl

@@ -128,6 +128,39 @@ using Aqua: Aqua
         @test mapreduce(iszero, &, values(diff))
     end # @testset "vector space interface"
 
+    @testset "order_of_accuracy" begin
+        @testset "RK4, rational coefficients" begin
+            # classical fourth-order Runge-Kutta method
+            A = [0 0 0 0;
+                 1//2 0 0 0;
+                 0 1//2 0 0;
+                 0 0 1 0]
+            b = [1 // 6, 1 // 3, 1 // 3, 1 // 6]
+            rk = RungeKuttaMethod(A, b)
+            series = bseries(rk, 5)
+            @test @inferred(order_of_accuracy(series)) == 4
+
+            series = bseries(rk, 3)
+            @test @inferred(order_of_accuracy(series)) == 3
+        end
+
+        @testset "RK4, floating point coefficients" begin
+            # classical fourth-order Runge-Kutta method
+            A = [0 0 0 0;
+                 1/2 0 0 0;
+                 0 1/2 0 0;
+                 0 0 1 0]
+            b = [1 / 6, 1 / 3, 1 / 3, 1 / 6]
+            rk = RungeKuttaMethod(A, b)
+            series = bseries(rk, 5)
+            @test @inferred(order_of_accuracy(series)) == 4
+            @test @inferred(order_of_accuracy(series; atol = 10 * eps())) == 4
+
+            series = bseries(rk, 3)
+            @test @inferred(order_of_accuracy(series)) == 3
+        end
+    end
+
     @testset "substitute" begin
         Symbolics.@variables a1 a2 a31 a32
         a = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Symbolics.Num}(rootedtree(Int[]) => 1,
@@ -1104,6 +1137,7 @@ using Aqua: Aqua
 
         diff = series - ExactSolution(series)
         @test mapreduce(abs ∘ last, +, diff) == 1 // 12
+        @test @inferred(order_of_accuracy(series)) == 2
     end
 
     @testset "Runge-Kutta methods interface" begin
@@ -1119,11 +1153,13 @@ using Aqua: Aqua
         series_integrator = @inferred bseries(rk, 4)
         series_exact = @inferred ExactSolution(series_integrator)
         @test mapreduce(==, &, values(series_integrator), values(series_exact))
+        @test @inferred(order_of_accuracy(series_integrator)) == 4
 
         # not fifth-order accurate
         series_integrator = @inferred bseries(rk, 5)
         series_exact = @inferred ExactSolution(series_integrator)
         @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false
+        @test @inferred(order_of_accuracy(series_integrator)) == 4
     end # @testset "Runge-Kutta methods interface"
 
     @testset "additive Runge-Kutta methods interface" begin
@@ -1136,6 +1172,7 @@ using Aqua: Aqua
             series_integrator = @inferred bseries(ark, 1)
             series_exact = @inferred ExactSolution(series_integrator)
             @test mapreduce(==, &, values(series_integrator), values(series_exact))
+            @test @inferred(order_of_accuracy(series_integrator)) == 1
 
             # not second-order accurate
             series_integrator = @inferred bseries(ark, 2)
@@ -1198,6 +1235,7 @@ using Aqua: Aqua
             series_integrator = @inferred bseries(ark, 2)
             series_exact = @inferred ExactSolution(series_integrator)
             @test mapreduce(==, &, values(series_integrator), values(series_exact))
+            @test @inferred(order_of_accuracy(series_integrator)) == 2
 
             # not third-order accurate
             series_integrator = @inferred bseries(ark, 3)
@@ -1306,6 +1344,7 @@ using Aqua: Aqua
             series_integrator = @inferred bseries(ark, 3)
             series_exact = @inferred ExactSolution(series_integrator)
             @test mapreduce(==, &, values(series_integrator), values(series_exact))
+            @test @inferred(order_of_accuracy(series_integrator)) == 3
 
             # not fourth-order accurate
             series_integrator = @inferred bseries(ark, 4)