order_of_symplecticity and is_symplectic (#216)

* update .gitignore

* order_of_symplecticity and is_symplectic

* first tests

* format

* kwargs for is_symplectic

* test AVF

* symbolic tests

* more tests

* test Gauss methods

* format

* pseudo-symplectic method

* ref.

* more inferred tests

* WIP: buffers for butcher_product! ???

* no buffers required

* skip failing test

* bump version of RootedTrees

* add TODO notes for colored trees

* test consistency

* benchmark symplectic check with pybs

* fix tests on Julia v1.6

* fix highlighing of md

* fix test

* fix

.gitignore

@@ -5,10 +5,15 @@
 **/.ipynb_checkpoints/
 **/*.gif
 **/Manifest.toml
+
 docs/build
 docs/src/contributing.md
 docs/src/license.md
 docs/src/benchmark_python_*.txt
 docs/dev/*.pdf
 docs/dev/*.png
 
+run*/
+**/.CondaPkg/
+
+**/*.code-workspace

Project.toml

@@ -1,7 +1,7 @@
 name = "BSeries"
 uuid = "ebb8d67c-85b4-416c-b05f-5f409e808f32"
 authors = ["Hendrik Ranocha <[email protected]> and contributors"]
-version = "0.1.61"
+version = "0.1.62"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -34,7 +34,7 @@ OrderedCollections = "1"
 Polynomials = "2.0.23, 3 - 3.2.8, 3.2.10"
 Reexport = "1"
 Requires = "1"
-RootedTrees = "2.16"
+RootedTrees = "2.23"
 SparseArrays = "1"
 SymEngine = "0.8, 0.9.1, 0.10, 0.11"
 SymPy = "2"

docs/src/benchmark_python_pybs.py

@@ -63,3 +63,23 @@ def first_values(f, n):
 with open(io_file, 'a') as io:
   print(result, file=io)
   print("", end_time - start_time, "seconds", file=io)
+
+
+with open(io_file, 'a') as io:
+  print("\nSymplecticity (conservation of quadratic invariants)", file=io)
+
+start_time = time.time()
+a = rk_methods.RKimplicitMidpoint.phi()
+result = pybs.series.symplectic_up_to_order(a, up_to_order)
+end_time = time.time()
+with open(io_file, 'a') as io:
+  print(result, file=io)
+  print("", end_time - start_time, "seconds", file=io)
+
+start_time = time.time()
+a = rk_methods.RKimplicitMidpoint.phi()
+result = pybs.series.symplectic_up_to_order(a, up_to_order)
+end_time = time.time()
+with open(io_file, 'a') as io:
+  print(result, file=io)
+  print("", end_time - start_time, "seconds", file=io)

docs/src/benchmarks.md

@@ -118,6 +118,7 @@ of the Python packages
 
 - [`BSeries`](https://github.com/ketch/BSeries)
 - [`pybs`](https://github.com/henriksu/pybs)
+- [`orderconditions`](https://gitlab.com/v_dallerit/orderconditions)
 
 If you know about similar open source packages out there, please inform us, e.g.,
 by [creating an issue](https://github.com/ranocha/BSeries.jl/issues/new/choose)
@@ -304,6 +305,26 @@ end
   series = bseries(AverageVectorFieldMethod(), up_to_order)
   println(is_energy_preserving(series))
 end
+
+
+println("\nSymplecticity (conservation of quadratic invariants)")
+@time begin
+  # implicit midpoint method = first Gauss method
+  A = @SArray [1//2;;]
+  b = @SArray [1//1]
+  rk = RungeKuttaMethod(A, b)
+  series = bseries(rk, up_to_order)
+  println(is_symplectic(series))
+end
+
+@time begin
+  # implicit midpoint method = first Gauss method
+  A = @SArray [1//2;;]
+  b = @SArray [1//1]
+  rk = RungeKuttaMethod(A, b)
+  series = bseries(rk, up_to_order)
+  println(is_symplectic(series))
+end
 ```
 
 

src/BSeries.jl

@@ -44,6 +44,8 @@ export renormalize!
 
 export is_energy_preserving, energy_preserving_order
 
+export order_of_symplecticity, is_symplectic
+
 # Types used for traits
 # These traits may decide between different algorithms based on the
 # corresponding complexity etc.
@@ -330,11 +332,13 @@ end
 
 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 performed using
-`isapprox` and the keyword arguments `kwargs...` are forwarded.
+`isequal`. If keyword arguments such as absolute/relative tolerances
+`atol`/`rtol` are given or floating point numbers are used, the comparison
+is performed using `isapprox` and the keyword arguments `kwargs...` are
+forwarded.
 
-See also [`order`](@ref), [`ExactSolution`](@ref).
+See also [`order`](@ref), [`ExactSolution`](@ref),
+[`order_of_symplecticity`](@ref).
 """
 function order_of_accuracy(series::TruncatedBSeries; kwargs...)
     if isempty(kwargs) && !(valtype(series) <: AbstractFloat)
@@ -1770,7 +1774,7 @@ energy-preserving for Hamiltonian problems.
 It requires a `max_order` so that it does not run forever if the order up to
 which the method is energy-preserving is too big or infinite.
 
-See also [`is_energy_preserving`](@ref)
+See also [`is_energy_preserving`](@ref).
 """
 function energy_preserving_order(rk::RungeKuttaMethod, max_order)
     p = 0
@@ -2233,4 +2237,127 @@ function equivalent_trees(tree)
     return equivalent_trees_set
 end
 
+#####################################################################
+# Symplectic = preserving quadratic invariants
+
+"""
+    satisfied_for_trees_of_order(condition, series, order,
+                                 iterator = RootedTreeIterator)
+
+Checks whether a given `condition` is satisfied for all pairs of trees
+`t1` and `t2` with given `order == order(t1) + order(t2)` for a given
+`series`. Which trees are considered is determined by the `iterator`.
+
+The `condition` is called as `condition(series, t1, t2)` and should return
+`true` if the condition is satisfied and `false` otherwise.
+"""
+function satisfied_for_trees_up_to_order(condition, series, order,
+                                         iterator = RootedTreeIterator)
+    for o1 in 1:(order - 1)
+        o2 = order - o1
+        for t1 in iterator(o1)
+            for t2 in iterator(o2)
+                if !condition(series, t1, t2)
+                    return false
+                end
+            end
+        end
+    end
+    return true
+end
+
+"""
+    order_of_symplecticity(series_integrator; kwargs...)
+
+Determine the order of symplecticity of the B-series `series_integrator`,
+i.e., the order up to which quadratic invariants are conserved.
+By default, the comparison of the coefficients entering the conditions 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 performed using `isapprox` and the keyword arguments
+`kwargs...` are forwarded.
+
+See also [`is_symplectic`](@ref),
+[`order`](@ref), [`order_of_accuracy`](@ref).
+"""
+function order_of_symplecticity(series::TruncatedBSeries; kwargs...)
+    # TODO: Implement this for colored trees
+    #       Theorem IV.7.2 of Hairer, Lubich, Wanner (2006) states the
+    #       following conditions.
+    #       1) If the coefficients a(t) satisfy
+    #            a(u ∘ v) + a(v ∘ u) = a(u) * a(v) for u ∈ TP_p, v ∈ TP_q,
+    #          and
+    #            a(τ) is independent of the colour of the root of τ,
+    #          the method exactly conserves quadratic invariants Q(p, q)
+    #          and it is symplectic for general Hamiltonian systems.
+    #        2) If the coefficients a(τ) satisfy only
+    #             a(u ∘ v) + a(v ∘ u) = a(u) * a(v) for u ∈ TP_p, v ∈ TP_q,
+    #           the method exactly conserves quadratic invariants Q(p,q)
+    #           for problems of the form p' = f1(q), q' = f2(p), and it
+    #           is symplectic for separable Hamiltonian systems where
+    #           H(p,q) = T(p) + U(q).
+    #        How do we want to handle the two cases? The general interface
+    #        should be that we check the first condition by default. Then,
+    #        we need to different `iterator`s for the two cases.
+    #        It should also be possible to specify that we want to check
+    #        only the second case, i.e., separable Hamiltonian systems.
+    if !(keytype(series) <: RootedTree)
+        throw(ArgumentError("This method is only implemented for single-colored rooted trees"))
+    end
+
+    if isempty(kwargs) && !(valtype(series) <: AbstractFloat)
+        compare = isequal
+    else
+        compare = (a, b) -> isapprox(a, b; kwargs...)
+    end
+
+    t12 = copy(first(keys(series)))
+    t21 = copy(first(keys(series)))
+
+    condition = let compare = compare, t12 = t12, t21 = t21
+        function (series, t1, t2)
+            butcher_product!(t12, t1, t2)
+            butcher_product!(t21, t2, t1)
+            compare(series[t12] + series[t21], series[t1] * series[t2])
+        end
+    end
+    iterator = _iterator_type(series)
+
+    if order(series) < 1
+        return 0
+    end
+
+    t = first(iterator(0))
+    if !compare(series[t], one(valtype(series)))
+        return 0
+    end
+
+    for o in 1:order(series)
+        if !satisfied_for_trees_up_to_order(condition, series, o, iterator)
+            return o - 1
+        end
+    end
+
+    return order(series)
+end
+
+"""
+    is_symplectic(series_integrator; kwargs...)::Bool
+
+This function checks whether the B-series `series_integrator` of a time
+integration method is symplectic (conserves quadratic invariants) - up to the
+[`order`](@ref) of `series_integrator`.
+
+By default, the comparison of the coefficients entering the conditions 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 performed using `isapprox` and the keyword arguments
+`kwargs...` are forwarded.
+
+See also [`order_of_symplecticity`](@ref).
+"""
+function is_symplectic(series::TruncatedBSeries; kwargs...)
+    order_of_symplecticity(series; kwargs...) == order(series)
+end
+
 end # module