diff --git a/Project.toml b/Project.toml index 092f265c..4c56c5fa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "BSeries" uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32" authors = ["Hendrik Ranocha and contributors"] -version = "0.1.28" +version = "0.1.29" [deps] Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" diff --git a/docs/src/tutorials/bseries_basics.md b/docs/src/tutorials/bseries_basics.md index 2191610a..3a240b75 100644 --- a/docs/src/tutorials/bseries_basics.md +++ b/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) diff --git a/docs/src/tutorials/bseries_creation.md b/docs/src/tutorials/bseries_creation.md index 38017a71..c1e90f96 100644 --- a/docs/src/tutorials/bseries_creation.md +++ b/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) ``` diff --git a/src/BSeries.jl b/src/BSeries.jl index c8e31be2..3de03da6 100644 --- a/src/BSeries.jl +++ b/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) diff --git a/test/runtests.jl b/test/runtests.jl index a43f40c7..474c8cde 100644 --- a/test/runtests.jl +++ b/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)