diff --git a/.gitignore b/.gitignore index 74170cb2..268fc6c4 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,7 @@ **/.ipynb_checkpoints/ **/*.gif **/Manifest.toml + docs/build docs/src/contributing.md docs/src/license.md @@ -12,3 +13,7 @@ docs/src/benchmark_python_*.txt docs/dev/*.pdf docs/dev/*.png +run*/ +**/.CondaPkg/ + +**/*.code-workspace diff --git a/Project.toml b/Project.toml index 0899a461..8463da96 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.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" diff --git a/docs/src/benchmark_python_pybs.py b/docs/src/benchmark_python_pybs.py index 6859dbe4..c3862052 100644 --- a/docs/src/benchmark_python_pybs.py +++ b/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) diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md index b8072392..12fbe88c 100644 --- a/docs/src/benchmarks.md +++ b/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 ``` diff --git a/src/BSeries.jl b/src/BSeries.jl index d19385f7..010da84f 100644 --- a/src/BSeries.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 78107f5e..4577d992 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -21,10 +21,10 @@ using Aqua: Aqua @testset "BSeries" begin @testset "lazy representation of exact ODE solution" begin - exact = ExactSolution{Rational{Int}}() + exact = @inferred ExactSolution{Rational{Int}}() terms = collect(Iterators.take(exact, 5)) @test terms == [1 // 1, 1 // 1, 1 // 2, 1 // 6, 1 // 3] - @test exact == ExactSolution(exact) + @test exact == @inferred ExactSolution(exact) end @testset "non-conflicting exports" begin @@ -45,7 +45,7 @@ using Aqua: Aqua b = @SArray [0, 1 // 1] c = @SArray [0, 1 // 2] - series_integrator = bseries(A, b, c, 2) + series_integrator = @inferred bseries(A, b, c, 2) @test_nowarn latexify(series_integrator) @test_nowarn latexify(series_integrator, cdot = false) @test_nowarn latexify(series_integrator, dt = SymEngine.symbols("h")) @@ -2382,6 +2382,279 @@ using Aqua: Aqua end end + @testset "Symplecticity (preserving quadratic invariants)" begin + @testset "Classical RK4 method" begin + @testset "rational coefficients" begin + A = [0//1 0//1 0//1 0//1 + 1//2 0//1 0//1 0//1 + 0//1 1//2 0//1 0//1 + 0//1 0//1 1//1 0//1] + b = [1 // 6, 1 // 3, 1 // 3, 1 // 6] + rk = RungeKuttaMethod(A, b) + series = bseries(rk, 6) + @test @inferred(order_of_symplecticity(series)) == 4 + @test @inferred(is_symplectic(series)) == false + end + + @testset "floating point coefficients" begin + A = [0 0 0 0 + 0.5 0 0 0 + 0 0.5 0 0 + 0 0 1 0] + b = [1 / 6, 1 / 3, 1 / 3, 1 / 6] + rk = RungeKuttaMethod(A, b) + series = bseries(rk, 6) + @test @inferred(order_of_symplecticity(series)) == 4 + @test @inferred(is_symplectic(series)) == false + end + + @testset "rational coefficients, kwargs" begin + A = [0//1 0//1 0//1 0//1 + 1//2 0//1 0//1 0//1 + 0//1 1//2 0//1 0//1 + 0//1 0//1 1//1 0//1] + b = [1 // 6, 1 // 3, 1 // 3, 1 // 6] + rk = RungeKuttaMethod(A, b) + series = bseries(rk, 6) + # With big tolerances, the series is "symplectic" + @test @inferred(order_of_symplecticity(series; atol = 1)) == 6 + @test @inferred(is_symplectic(series; atol = 1)) == true + @test @inferred(order_of_symplecticity(series; rtol = 0.2)) == 5 + @test @inferred(is_symplectic(series; rtol = 0.2)) == false + end + end + + @testset "Average Vector Field (AVF)" begin + series = @inferred bseries(AverageVectorFieldMethod(), 6) + @test @inferred(order_of_symplecticity(series)) == + @inferred(order_of_accuracy(series)) == 2 + @test @inferred(is_symplectic(series)) == false + + @testset "$T" for T in [Float32, Float64, BigFloat] + series = @inferred bseries(AverageVectorFieldMethod(T), 6) + @test @inferred(order_of_symplecticity(series)) == + @inferred(order_of_accuracy(series)) == 2 + @test @inferred(is_symplectic(series)) == false + end + end + + @testset "Not a series of a consistent integrator" begin + series = @inferred bseries(AverageVectorFieldMethod(), 6) + series[rootedtree(Int[])] = 0 + @test @inferred(order_of_symplecticity(series)) == 0 + @test @inferred(is_symplectic(series)) == false + end + + @testset "Implicit midpoint method (symplectic)" begin + @testset "rational coefficients" begin + A = @SMatrix [1 // 2;;] + b = @SVector [1] + rk = RungeKuttaMethod(A, b) + series = @inferred bseries(rk, 9) + @test @inferred(order_of_accuracy(series)) == 2 + @test @inferred is_symplectic(series) + end + + @testset "floating point coefficients" begin + A = @SMatrix [0.5;;] + b = @SVector [1] + rk = RungeKuttaMethod(A, b) + series = @inferred bseries(rk, 9) + @test @inferred(order_of_accuracy(series)) == 2 + @test @inferred is_symplectic(series) + end + end + + @testset "Gauss method (s = 2)" begin + # Butcher (2016) + # Numerical methods for ordinary differential equations + # Section 342 + A = [1/4 (1 / 4-sqrt(3) / 6); + (1 / 4+sqrt(3) / 6) 1/4] + b = [1 / 2, 1 / 2] + rk = @inferred RungeKuttaMethod(A, b) + series = @inferred bseries(rk, 9) + @test @inferred(order_of_accuracy(series)) == 4 + atol = eps(valtype(series)) + rtol = sqrt(eps(valtype(series))) + @test @inferred(is_symplectic(series; atol = atol, rtol = rtol)) + end + + @testset "Gauss method (s = 3)" begin + # Butcher (2016) + # Numerical methods for ordinary differential equations + # Section 342 + A = [5/36 2 / 9-sqrt(15) / 15 5 / 36-sqrt(15) / 30; + 5 / 36+sqrt(15) / 24 2/9 5 / 36-sqrt(15) / 24; + 5 / 36+sqrt(15) / 30 2 / 9+sqrt(15) / 15 5/36] + b = [5 / 18, 4 / 9, 5 / 18] + rk = @inferred RungeKuttaMethod(A, b) + series = @inferred bseries(rk, 9) + @test @inferred(order_of_accuracy(series)) == 6 + atol = eps(valtype(series)) + rtol = sqrt(eps(valtype(series))) + @test @inferred(is_symplectic(series; atol = atol, rtol = rtol)) + end + + @testset "Pseudo-symplectic method PS(2, 4, 2)" begin + # Aubry, Chartier (1998) + # Pseudo-symplectic Runge-Kutta methods + k = 3 + A = [0 0 0; + (8 * k - 3)//(8k - 4) 0 0; + (16k^2 - 8 * k + 1)//(2 * k * (8k - 3)) (2k - 1)//(2 * k * (8k - 3)) 0] + b = [(3k - 1) // (8k - 3), (-2 * (2k - 1)^2) // (8k - 3), k] + c = [0, (8k - 3) // (8k - 4), 1] + rk = @inferred RungeKuttaMethod(A, b) + series = @inferred bseries(rk, 6) + @test @inferred(order_of_accuracy(series)) == 2 + @test @inferred(order_of_symplecticity(series)) == 4 + @test @inferred(is_symplectic(series)) == false + end + + @testset "Symplectic Euler method" begin + @testset "rational coefficients" begin + if VERSION >= v"1.10" + ex_euler = @inferred RungeKuttaMethod(@SMatrix([0 // 1]), @SVector [1]) + im_euler = @inferred RungeKuttaMethod(@SMatrix([1 // 1]), @SVector [1]) + else + # On Julia v1.6, the `@inferred` check fails in CI + ex_euler = RungeKuttaMethod(@SMatrix([0 // 1]), @SVector [1]) + im_euler = RungeKuttaMethod(@SMatrix([1 // 1]), @SVector [1]) + end + ark = @inferred AdditiveRungeKuttaMethod([ex_euler, im_euler]) + series = @inferred bseries(ark, 9) + # needs to be implemented for colored rooted trees + @test_throws ArgumentError @inferred is_symplectic(series) + end + + @testset "floating point coefficients" begin + if VERSION >= v"1.10" + ex_euler = @inferred RungeKuttaMethod(@SMatrix([0.0]), @SVector [1.0]) + im_euler = @inferred RungeKuttaMethod(@SMatrix([1.0]), @SVector [1.0]) + else + # On Julia v1.6, the `@inferred` check fails in CI + ex_euler = RungeKuttaMethod(@SMatrix([0.0]), @SVector [1.0]) + im_euler = RungeKuttaMethod(@SMatrix([1.0]), @SVector [1.0]) + end + ark = @inferred AdditiveRungeKuttaMethod([ex_euler, im_euler]) + series = @inferred bseries(ark, 9) + # needs to be implemented for colored rooted trees + @test_skip @inferred is_symplectic(series) + end + end + + @testset "Störmer-Verlet" begin + @testset "rational coefficients" begin + # Hairer, Lubich, Wanner (2002) + # Geometric numerical integration + # Table II.2.1 + As = [ + [0 0; 1//2 1//2], + [1//2 0; 1//2 0], + ] + bs = [ + [1 // 2, 1 // 2], + [1 // 2, 1 // 2], + ] + ark = AdditiveRungeKuttaMethod(As, bs) + series = @inferred bseries(ark, 9) + # needs to be implemented for colored rooted trees + @test_skip @inferred is_symplectic(series) + end + + @testset "floating point coefficients" begin + # Hairer, Lubich, Wanner (2002) + # Geometric numerical integration + # Table II.2.1 + As = [ + [0.0 0.0; 0.5 0.52], + [0.5 0.0; 0.5 0.0], + ] + bs = [ + [0.5, 0.5], + [0.5, 0.5], + ] + ark = AdditiveRungeKuttaMethod(As, bs) + series = @inferred bseries(ark, 9) + # needs to be implemented for colored rooted trees + @test_skip @inferred is_symplectic(series) + end + end + + @testset "Lobatto IIIA-IIIB pair (s = 3)" begin + # Hairer, Lubich, Wanner (2002) + # Geometric numerical integration + # Table II.2.2 + As = [ + [0 0 0; 5//24 1//3 -1//24; 1//6 2//3 1//6], + [1//6 -1//6 0; 1//6 1//3 0; 1//6 5//6 0], + ] + bs = [ + [1 // 6, 2 // 3, 1 // 6], + [1 // 6, 2 // 3, 1 // 6], + ] + ark = AdditiveRungeKuttaMethod(As, bs) + series = @inferred bseries(ark, 9) + # needs to be implemented for colored rooted trees + @test_skip @inferred is_symplectic(series) + end + + @testset "Symbolic coefficients" begin + @testset "SymEngine.jl" begin + # This method is second-order accurate. Thus, it is + # symplectic up to order two. + α = SymEngine.symbols("α") + A = [0 0; 1/(2 * α) 0] + b = [1 - α, α] + c = [0, 1 / (2 * α)] + series_integrator = @inferred(bseries(A, b, c, 3)) + @test @inferred(order_of_accuracy(series_integrator)) == 2 + @test @inferred(order_of_symplecticity(series_integrator)) == 2 + @test @inferred(is_symplectic(series_integrator)) == false + end + + @testset "SymPy.jl" begin + # This method is second-order accurate. Thus, it is + # energy-preserving up to order two. + α = SymPy.symbols("α", real = true) + A = [0 0; 1/(2 * α) 0] + b = [1 - α, α] + c = [0, 1 / (2 * α)] + series_integrator = @inferred(bseries(A, b, c, 3)) + @test @inferred(order_of_accuracy(series_integrator)) == 2 + @test @inferred(order_of_symplecticity(series_integrator)) == 2 + @test @inferred(is_symplectic(series_integrator)) == false + end + + @testset "SymPyPythonCall.jl" begin + # This method is second-order accurate. Thus, it is + # energy-preserving up to order two. + α = SymPyPythonCall.symbols("α", real = true) + A = [0 0; 1/(2 * α) 0] + b = [1 - α, α] + c = [0, 1 / (2 * α)] + series_integrator = @inferred(bseries(A, b, c, 3)) + @test @inferred(order_of_accuracy(series_integrator)) == 2 + @test @inferred(order_of_symplecticity(series_integrator)) == 2 + @test @inferred(is_symplectic(series_integrator)) == false + end + + @testset "Symbolics.jl" begin + # This method is second-order accurate. Thus, it is + # energy-preserving up to order two. + Symbolics.@variables α + A = [0 0; 1/(2 * α) 0] + b = [1 - α, α] + c = [0, 1 / (2 * α)] + series_integrator = @inferred(bseries(A, b, c, 3)) + @test @inferred(order_of_accuracy(series_integrator)) == 2 + @test @inferred(order_of_symplecticity(series_integrator)) == 2 + @test @inferred(is_symplectic(series_integrator)) == false + end + end + end + @testset "Aqua" begin if VERSION < v"1.7" Aqua.test_all(BSeries;