From 8ee52041d71b6eea1d30773bd759ca3b56fec22c Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Tue, 5 Sep 2023 17:54:16 +0300 Subject: [PATCH 01/13] Absolute tolerance in new PR --- src/BSeries.jl | 40 +++++++++++++++++++++++----------------- test/runtests.jl | 12 ++++++++++++ 2 files changed, 35 insertions(+), 17 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 0d1d8c9e..3a9d3a08 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1748,16 +1748,18 @@ function renormalize!(series) end """ - energy_preserving_order(rk::RungeKuttaMethod, max_order) + energy_preserving_order(rk::RungeKuttaMethod, max_order; tol::Real=0) This function checks up to which order a Runge-Kutta method `rk` is 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. +Keyword Arguments: +- `tol::Real=1e-14`: The absolute tolerance for energy preservation. See also [`is_energy_preserving`](@ref) """ -function energy_preserving_order(rk::RungeKuttaMethod, max_order) +function energy_preserving_order(rk::RungeKuttaMethod, max_order; tol::Real=0) p = 0 not_energy_preserving = false while not_energy_preserving == false @@ -1766,7 +1768,7 @@ function energy_preserving_order(rk::RungeKuttaMethod, max_order) return max_order end # Check energy preservation up to the given order - if is_energy_preserving(rk, p + 1) == false + if is_energy_preserving(rk, p + 1, tol=tol) == false not_energy_preserving = true end p = p + 1 @@ -1775,22 +1777,26 @@ function energy_preserving_order(rk::RungeKuttaMethod, max_order) end """ - is_energy_preserving(rk::RungeKuttaMethod, order)::Bool + is_energy_preserving(rk::RungeKuttaMethod, order; tol::Real=0)::Bool This function checks whether the Runge-Kutta method `rk` is energy-preserving for Hamiltonian systems up to a given `order`. +Keyword Arguments: +- `tol::Real=1e-14`: The absolute tolerance for energy preservation. """ -function is_energy_preserving(rk::RungeKuttaMethod, order) +function is_energy_preserving(rk::RungeKuttaMethod, order; tol::Real=0) series = bseries(rk, order) - return is_energy_preserving(series) + return is_energy_preserving(series, tol=tol) end """ - is_energy_preserving(series_integrator)::Bool + is_energy_preserving(series_integrator; tol::Real=0)::Bool This function checks whether the B-series `series_integrator` of a time integration method is energy-preserving for Hamiltonian systems - up to the [`order`](@ref) of `series_integrator`. +Keyword Arguments: +- `tol::Real=1e-14`: The absolute tolerance for energy preservation. # References @@ -1800,7 +1806,7 @@ This code is based on the Theorem 2 of Foundations of Computational Mathematics 10 (2010): 673-693. [DOI: 10.1007/s10208-010-9073-1](https://link.springer.com/article/10.1007/s10208-010-9073-1) """ -function is_energy_preserving(series_integrator) +function is_energy_preserving(series_integrator; tol::Real=0) # This method requires the B-series of a map as input let t = first(keys(series_integrator)) @assert isempty(t) @@ -1835,7 +1841,7 @@ function is_energy_preserving(series_integrator) # for order less than 4 because the energy-preserving basis has only one element. # The following function checks the energy-preserving condition up to order 4 # together. - if _is_energy_preserving_low_order(trees, coefficients) == false + if _is_energy_preserving_low_order(trees, coefficients, tol) == false return false end @@ -1852,7 +1858,7 @@ function is_energy_preserving(series_integrator) push!(same_order_trees, trees[j]) end end - if _is_energy_preserving(same_order_trees, same_order_coeffs) == false + if _is_energy_preserving(same_order_trees, same_order_coeffs, tol) == false return false end end @@ -1864,7 +1870,7 @@ end # and the set of 'coefficients' corresponding to a B-series is # energy_preserving up to the 'uppermost_order', which we choose # to be 4 for computational optimization -function _is_energy_preserving_low_order(trees, coefficients) +function _is_energy_preserving_low_order(trees, coefficients, tol) # Set the limit up to which this function will work uppermost_order = 4 same_order_trees = empty(trees) @@ -1877,7 +1883,7 @@ function _is_energy_preserving_low_order(trees, coefficients) push!(same_order_trees, t) push!(same_order_coeffs, coefficients[index]) end - return _is_energy_preserving(same_order_trees, same_order_coeffs) + return _is_energy_preserving(same_order_trees, same_order_coeffs, tol) end # This function is the application of Theorem 2 of the paper @@ -1886,7 +1892,7 @@ end # of coefficients. # Checks whether `trees` and `ocefficients` satisfify the energy_preserving # condition. -function _is_energy_preserving(trees, coefficients) +function _is_energy_preserving(trees, coefficients, tol) # TODO: `Float32` would also be nice to have. However, the default tolerance # of the rank computation is based on `Float64`. Thus, it will usually # not work with coefficients given only in 32 bit precision. @@ -1897,13 +1903,13 @@ function _is_energy_preserving(trees, coefficients) Rational{Int8}, Rational{Int16}, Rational{Int32}, Rational{Int64}, Rational{Int128}}) # These types support efficient computations in sparse matrices - _is_energy_preserving_sparse(trees, coefficients) + _is_energy_preserving_sparse(trees, coefficients, tol) else - _is_energy_preserving_dense(trees, coefficients) + _is_energy_preserving_dense(trees, coefficients, tol) end end -function _is_energy_preserving_dense(trees, coefficients) +function _is_energy_preserving_dense(trees, coefficients, tol) # For every tree, obtain the adjoint and check if it exists length_coeff = length(trees) # Provided that a tree can have several adjoints, for every tree `t` @@ -1969,7 +1975,7 @@ end # This method is specialized for coefficients that can be used efficiently in # sparse arrays. -function _is_energy_preserving_sparse(trees, coefficients) +function _is_energy_preserving_sparse(trees, coefficients, tol) # For every tree, obtain the adjoint and check if it exists. # Provided that a tree can have several adjoints, for every tree `t` # we need to check if there is a linear combination of the coefficients diff --git a/test/runtests.jl b/test/runtests.jl index acc0c234..a7077754 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2115,6 +2115,18 @@ using Aqua: Aqua series = bseries(AverageVectorFieldMethod(BigFloat), 7) # TODO: This test is currently broken and throws an error @test_broken is_energy_preserving(series) + + # Generated from MatLab + A = [0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 + -2.0314781009861151E-02 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 + 1.0721250651661693E+01 -1.0078547105299583E+01 0.0000000000000000E+00 0.0000000000000000E+00 + -4.1299845612430317E+01 4.0239538957660457E+01 2.0603066547698594E+00 0.0000000000000000E+00] + b = [2.0623108969503225E+00 + -1.7306596124709785E+00 + 5.6957376564633877E-01 + 9.8774949874317328E-02 ] + rk = RungeKuttaMethod(A,b) + @test energy_preserving_order(rk, 10) == 4 end @testset "Symbolic coefficients" begin From 990df436d90d3ce6a5cba0379291c67381e6f011 Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Tue, 5 Sep 2023 18:02:44 +0300 Subject: [PATCH 02/13] Tolerances for sparse and dense matrixes --- src/BSeries.jl | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 3a9d3a08..488be55d 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1927,8 +1927,19 @@ function _is_energy_preserving_dense(trees, coefficients, tol) # If so, we can exit early since the coefficients of bushy trees # must be zero in the modified equation of energy-preserving methods. if length(t) > 1 && is_bushy(t) - if !iszero(coefficients[t_index]) - return false + if (typeof(coefficients[t_index]) == Float64 || typeof(coefficients[t_index]) == Float32 || typeof(coefficients[t_index]) == Float16) && tol == 0 + tol = 1e-14 + if !isapprox(coefficients[t_index], 0, atol = tol) + return false + end + elseif tol != 0 + if !isapprox(coefficients[t_index], 0, atol = tol) + return false + end + else + if !iszero(coefficients[t_index]) + return false + end end else # If the tree is not a bush, then generate all the equivalent trees @@ -1997,8 +2008,19 @@ function _is_energy_preserving_sparse(trees, coefficients, tol) # If so, we can exit early since the coefficients of bushy trees # must be zero in the modified equation of energy-preserving methods. if length(t) > 1 && is_bushy(t) - if !iszero(coefficients[t_index]) - return false + if (typeof(coefficients[t_index]) == Float64 || typeof(coefficients[t_index]) == Float32 || typeof(coefficients[t_index]) == Float16) && tol == 0 + tol = 1e-14 + if !isapprox(coefficients[t_index], 0, atol = tol) + return false + end + elseif tol != 0 + if !isapprox(coefficients[t_index], 0, atol = tol) + return false + end + else + if !iszero(coefficients[t_index]) + return false + end end else # If the tree is not a bush, then generate all the equivalent trees From fcc646d2956e05b1da44f4789436a02076ac2a74 Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Wed, 6 Sep 2023 20:53:18 +0300 Subject: [PATCH 03/13] New function 'energy_preserving_default_tolerance' --- src/BSeries.jl | 56 +++++++++++++++++++++--------------------------- test/runtests.jl | 5 ++++- 2 files changed, 29 insertions(+), 32 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 488be55d..f805e984 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1755,7 +1755,8 @@ 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. Keyword Arguments: -- `tol::Real=1e-14`: The absolute tolerance for energy preservation. +- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +`rk` are Floating Points, the default tolerance is set to tol = 1e-14. See also [`is_energy_preserving`](@ref) """ @@ -1782,7 +1783,8 @@ end This function checks whether the Runge-Kutta method `rk` is energy-preserving for Hamiltonian systems up to a given `order`. Keyword Arguments: -- `tol::Real=1e-14`: The absolute tolerance for energy preservation. +- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +`rk` are Floating Points, the default tolerance is set to tol = 1e-14. """ function is_energy_preserving(rk::RungeKuttaMethod, order; tol::Real=0) series = bseries(rk, order) @@ -1796,7 +1798,8 @@ This function checks whether the B-series `series_integrator` of a time integration method is energy-preserving for Hamiltonian systems - up to the [`order`](@ref) of `series_integrator`. Keyword Arguments: -- `tol::Real=1e-14`: The absolute tolerance for energy preservation. +- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +`rk` are Floating Points, the default tolerance is set to tol = 1e-14. # References @@ -1813,6 +1816,9 @@ function is_energy_preserving(series_integrator; tol::Real=0) @assert isone(series_integrator[t]) end + # tolerance + V = eltype(coefficients) + tol = energy_preserving_default_tolerance(V, tol) # Theorem 2 of Celledoni et al. (2010) requires working with the modified # equation. The basic idea is to check whether the modified equation lies # in a certain subspace spanned by energy-preserving pairs of trees. @@ -1922,24 +1928,12 @@ function _is_energy_preserving_dense(trees, coefficients, tol) for (t_index, t) in enumerate(trees) # We do not need to consider the empty tree isempty(t) && continue - # Check whether the level sequence corresponds to a bushy tree. # If so, we can exit early since the coefficients of bushy trees # must be zero in the modified equation of energy-preserving methods. if length(t) > 1 && is_bushy(t) - if (typeof(coefficients[t_index]) == Float64 || typeof(coefficients[t_index]) == Float32 || typeof(coefficients[t_index]) == Float16) && tol == 0 - tol = 1e-14 - if !isapprox(coefficients[t_index], 0, atol = tol) - return false - end - elseif tol != 0 - if !isapprox(coefficients[t_index], 0, atol = tol) - return false - end - else - if !iszero(coefficients[t_index]) - return false - end + if !isapprox(coefficients[t_index], 0, atol = tol) + return false end else # If the tree is not a bush, then generate all the equivalent trees @@ -2003,24 +1997,12 @@ function _is_energy_preserving_sparse(trees, coefficients, tol) for (t_index, t) in enumerate(trees) # We do not need to consider the empty tree isempty(t) && continue - # Check whether the level sequence corresponds to a bushy tree. # If so, we can exit early since the coefficients of bushy trees # must be zero in the modified equation of energy-preserving methods. if length(t) > 1 && is_bushy(t) - if (typeof(coefficients[t_index]) == Float64 || typeof(coefficients[t_index]) == Float32 || typeof(coefficients[t_index]) == Float16) && tol == 0 - tol = 1e-14 - if !isapprox(coefficients[t_index], 0, atol = tol) - return false - end - elseif tol != 0 - if !isapprox(coefficients[t_index], 0, atol = tol) - return false - end - else - if !iszero(coefficients[t_index]) - return false - end + if !isapprox(coefficients[t_index], 0, atol = tol) + return false end else # If the tree is not a bush, then generate all the equivalent trees @@ -2246,4 +2228,16 @@ function equivalent_trees(tree) return equivalent_trees_set end +# This function calculates the default tolernace for 'is_energy_preserving()' +function energy_preserving_default_tolerance(V, tol) + if tol == 0 + if V <: AbstractFloat + tol = 100 * eps(T) + else + tol = zero(T) + end + end + return tol +end + end # module diff --git a/test/runtests.jl b/test/runtests.jl index a7077754..19063230 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2116,7 +2116,10 @@ using Aqua: Aqua # TODO: This test is currently broken and throws an error @test_broken is_energy_preserving(series) - # Generated from MatLab + # Generated using RK-Opt + # References: + # Ketcheson et al., (2020). RK-Opt: A package for the design of numerical ODE solvers. + # Journal of Open Source Software, 5(54), 2514, https://doi.org/10.21105/joss.02514 A = [0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 -2.0314781009861151E-02 0.0000000000000000E+00 0.0000000000000000E+00 0.0000000000000000E+00 1.0721250651661693E+01 -1.0078547105299583E+01 0.0000000000000000E+00 0.0000000000000000E+00 From c120e1ef8a8a7fe121412a9f25a671f64fa8b36d Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Wed, 6 Sep 2023 21:06:56 +0300 Subject: [PATCH 04/13] Correction --- src/BSeries.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index f805e984..487b0731 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1817,7 +1817,7 @@ function is_energy_preserving(series_integrator; tol::Real=0) end # tolerance - V = eltype(coefficients) + V = eltype(series_integrator[rootedtree([0])]) tol = energy_preserving_default_tolerance(V, tol) # Theorem 2 of Celledoni et al. (2010) requires working with the modified # equation. The basic idea is to check whether the modified equation lies @@ -2232,9 +2232,9 @@ end function energy_preserving_default_tolerance(V, tol) if tol == 0 if V <: AbstractFloat - tol = 100 * eps(T) + tol = 100 * eps(V) else - tol = zero(T) + tol = zero(V) end end return tol From a5541b3b1c9296b90bf7d5c363cc2c6ee556e1ba Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Thu, 7 Sep 2023 11:29:02 +0300 Subject: [PATCH 05/13] Tol = nothing --- src/BSeries.jl | 14 +++++++------- test/runtests.jl | 2 ++ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 487b0731..b8b22cdd 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1755,12 +1755,12 @@ 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. Keyword Arguments: -- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of `rk` are Floating Points, the default tolerance is set to tol = 1e-14. See also [`is_energy_preserving`](@ref) """ -function energy_preserving_order(rk::RungeKuttaMethod, max_order; tol::Real=0) +function energy_preserving_order(rk::RungeKuttaMethod, max_order; tol=nothing) p = 0 not_energy_preserving = false while not_energy_preserving == false @@ -1783,10 +1783,10 @@ end This function checks whether the Runge-Kutta method `rk` is energy-preserving for Hamiltonian systems up to a given `order`. Keyword Arguments: -- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of `rk` are Floating Points, the default tolerance is set to tol = 1e-14. """ -function is_energy_preserving(rk::RungeKuttaMethod, order; tol::Real=0) +function is_energy_preserving(rk::RungeKuttaMethod, order; tol=nothing) series = bseries(rk, order) return is_energy_preserving(series, tol=tol) end @@ -1798,7 +1798,7 @@ This function checks whether the B-series `series_integrator` of a time integration method is energy-preserving for Hamiltonian systems - up to the [`order`](@ref) of `series_integrator`. Keyword Arguments: -- `tol::Real=1e-0`: The absolute tolerance for energy preservation. If the elements of +- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of `rk` are Floating Points, the default tolerance is set to tol = 1e-14. # References @@ -1809,7 +1809,7 @@ This code is based on the Theorem 2 of Foundations of Computational Mathematics 10 (2010): 673-693. [DOI: 10.1007/s10208-010-9073-1](https://link.springer.com/article/10.1007/s10208-010-9073-1) """ -function is_energy_preserving(series_integrator; tol::Real=0) +function is_energy_preserving(series_integrator; tol=nothing) # This method requires the B-series of a map as input let t = first(keys(series_integrator)) @assert isempty(t) @@ -2230,7 +2230,7 @@ end # This function calculates the default tolernace for 'is_energy_preserving()' function energy_preserving_default_tolerance(V, tol) - if tol == 0 + if tol === nothing if V <: AbstractFloat tol = 100 * eps(V) else diff --git a/test/runtests.jl b/test/runtests.jl index 19063230..03c73de4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2129,6 +2129,8 @@ using Aqua: Aqua 5.6957376564633877E-01 9.8774949874317328E-02 ] rk = RungeKuttaMethod(A,b) + series = bseries(rk) + @test order_of_accuracy(series) == 4 @test energy_preserving_order(rk, 10) == 4 end From 1e9dc5e47cb0b1b2106ae5a6db03c48d36402278 Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Sun, 10 Sep 2023 12:01:10 +0300 Subject: [PATCH 06/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index b8b22cdd..0c88a5d5 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -2228,7 +2228,7 @@ function equivalent_trees(tree) return equivalent_trees_set end -# This function calculates the default tolernace for 'is_energy_preserving()' +# This function calculates the default tolerance for 'is_energy_preserving' function energy_preserving_default_tolerance(V, tol) if tol === nothing if V <: AbstractFloat From 619f261295bebf3299fe78ec129e782ecb733f9b Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Sun, 10 Sep 2023 12:01:21 +0300 Subject: [PATCH 07/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index 0c88a5d5..d55addb6 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -2230,6 +2230,7 @@ end # This function calculates the default tolerance for 'is_energy_preserving' function energy_preserving_default_tolerance(V, tol) + # We use `nothing` as default value if tol === nothing if V <: AbstractFloat tol = 100 * eps(V) From 3221375e990ffce4c60f86a43596efa2c399e692 Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Sun, 10 Sep 2023 12:01:53 +0300 Subject: [PATCH 08/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index d55addb6..9a19ce56 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -2233,6 +2233,7 @@ function energy_preserving_default_tolerance(V, tol) # We use `nothing` as default value if tol === nothing if V <: AbstractFloat + # For floating point numbers, we cannot expect exact calculations in general. tol = 100 * eps(V) else tol = zero(V) From 8969aa760458fedcd73b556a4435d659c5b59113 Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Sun, 10 Sep 2023 12:02:09 +0300 Subject: [PATCH 09/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index 9a19ce56..9c12ce80 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -2236,6 +2236,7 @@ function energy_preserving_default_tolerance(V, tol) # For floating point numbers, we cannot expect exact calculations in general. tol = 100 * eps(V) else + # If something like rational numbers are used, we do not expect numerical errors. tol = zero(V) end end From f0eae7061926897647884f4cd4af28b57fad37d9 Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Sun, 10 Sep 2023 12:02:22 +0300 Subject: [PATCH 10/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 9c12ce80..c8e2a077 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1817,7 +1817,7 @@ function is_energy_preserving(series_integrator; tol=nothing) end # tolerance - V = eltype(series_integrator[rootedtree([0])]) + V = valtype(series_integrator) tol = energy_preserving_default_tolerance(V, tol) # Theorem 2 of Celledoni et al. (2010) requires working with the modified # equation. The basic idea is to check whether the modified equation lies From 47791f698d603df1855bfcdbb06d2ac92d6dcbbe Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Sun, 10 Sep 2023 12:21:12 +0300 Subject: [PATCH 11/13] Update docstrings --- src/BSeries.jl | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index b8b22cdd..550f5ca6 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1748,15 +1748,16 @@ function renormalize!(series) end """ - energy_preserving_order(rk::RungeKuttaMethod, max_order; tol::Real=0) + energy_preserving_order(rk::RungeKuttaMethod, max_order; tol=nothing) This function checks up to which order a Runge-Kutta method `rk` is 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. Keyword Arguments: -- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of -`rk` are Floating Points, the default tolerance is set to tol = 1e-14. +- `tol = nothing`: The absolute tolerance for energy preservation. Default value +is `tol = 100 * eps(V)` if `valtype(rk) == AbstractFloat`, and `tol = zero(V)` +for rational numbers. See also [`is_energy_preserving`](@ref) """ @@ -1778,13 +1779,14 @@ function energy_preserving_order(rk::RungeKuttaMethod, max_order; tol=nothing) end """ - is_energy_preserving(rk::RungeKuttaMethod, order; tol::Real=0)::Bool + is_energy_preserving(rk::RungeKuttaMethod, order; tol=nothing)::Bool This function checks whether the Runge-Kutta method `rk` is energy-preserving for Hamiltonian systems up to a given `order`. Keyword Arguments: -- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of -`rk` are Floating Points, the default tolerance is set to tol = 1e-14. +- `tol = nothing`: The absolute tolerance for energy preservation. Default value +is `tol = 100 * eps(V)` if `valtype(rk) == AbstractFloat`, and `tol = zero(V)` +for rational numbers. """ function is_energy_preserving(rk::RungeKuttaMethod, order; tol=nothing) series = bseries(rk, order) @@ -1792,15 +1794,16 @@ function is_energy_preserving(rk::RungeKuttaMethod, order; tol=nothing) end """ - is_energy_preserving(series_integrator; tol::Real=0)::Bool + is_energy_preserving(series_integrator; tol=nothing)::Bool This function checks whether the B-series `series_integrator` of a time integration method is energy-preserving for Hamiltonian systems - up to the [`order`](@ref) of `series_integrator`. Keyword Arguments: -- `tol = nothing`: The absolute tolerance for energy preservation. If the elements of -`rk` are Floating Points, the default tolerance is set to tol = 1e-14. - +- `tol = nothing`: The absolute tolerance for energy preservation. Default value +is `tol = 100 * eps(V)` if `valtype(series_integrator) == AbstractFloat`, and `tol = zero(V)` +for rational numbers. + # References This code is based on the Theorem 2 of From b61fe070675ea7f41abdadce1b97270074e6bfb3 Mon Sep 17 00:00:00 2001 From: Sondar74 <79122393+Sondar74@users.noreply.github.com> Date: Tue, 12 Sep 2023 11:32:11 +0300 Subject: [PATCH 12/13] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 43783816..23c3ff14 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1756,7 +1756,7 @@ 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. Keyword Arguments: - `tol = nothing`: The absolute tolerance for energy preservation. Default value -is `tol = 100 * eps(V)` if `valtype(rk) == AbstractFloat`, and `tol = zero(V)` +is `tol = 100 * eps(V)` if `valtype(rk) <: AbstractFloat`, and `tol = zero(V)` for rational numbers. See also [`is_energy_preserving`](@ref) From 3c8bf2ffc39445789ee3bfa94d3fc8847a302445 Mon Sep 17 00:00:00 2001 From: Sondar74 Date: Tue, 12 Sep 2023 11:33:55 +0300 Subject: [PATCH 13/13] changed '==' for '<:' --- src/BSeries.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 43783816..1431e1ab 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1756,7 +1756,7 @@ 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. Keyword Arguments: - `tol = nothing`: The absolute tolerance for energy preservation. Default value -is `tol = 100 * eps(V)` if `valtype(rk) == AbstractFloat`, and `tol = zero(V)` +is `tol = 100 * eps(V)` if `valtype(rk) <: AbstractFloat`, and `tol = zero(V)` for rational numbers. See also [`is_energy_preserving`](@ref) @@ -1785,7 +1785,7 @@ This function checks whether the Runge-Kutta method `rk` is energy-preserving for Hamiltonian systems up to a given `order`. Keyword Arguments: - `tol = nothing`: The absolute tolerance for energy preservation. Default value -is `tol = 100 * eps(V)` if `valtype(rk) == AbstractFloat`, and `tol = zero(V)` +is `tol = 100 * eps(V)` if `valtype(rk) <: AbstractFloat`, and `tol = zero(V)` for rational numbers. """ function is_energy_preserving(rk::RungeKuttaMethod, order; tol=nothing) @@ -1801,7 +1801,7 @@ integration method is energy-preserving for Hamiltonian systems - up to the [`order`](@ref) of `series_integrator`. Keyword Arguments: - `tol = nothing`: The absolute tolerance for energy preservation. Default value -is `tol = 100 * eps(V)` if `valtype(series_integrator) == AbstractFloat`, and `tol = zero(V)` +is `tol = 100 * eps(V)` if `valtype(series_integrator) <: AbstractFloat`, and `tol = zero(V)` for rational numbers. # References