diff --git a/Project.toml b/Project.toml index 1be88ed0..13f5ce05 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.59" +version = "0.1.60" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" @@ -17,12 +17,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [weakdeps] SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" +SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [extensions] -SymEngineExt = "SymEngine" -SymPyExt = "SymPy" -SymbolicsExt = "Symbolics" +BSeriesSymEngineExt = "SymEngine" +BSeriesSymPyExt = "SymPy" +BSeriesSymPyPythonCallExt = "SymPyPythonCall" +BSeriesSymbolicsExt = "Symbolics" [compat] Combinatorics = "1" @@ -35,11 +37,13 @@ Requires = "1" RootedTrees = "2.16" SparseArrays = "1" SymEngine = "0.8, 0.9.1, 0.10, 0.11" -SymPy = "1" +SymPy = "2" +SymPyPythonCall = "0.3" Symbolics = "4, 5" julia = "1.6" [extras] SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" +SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" diff --git a/docs/Project.toml b/docs/Project.toml index 83df925c..b29307f0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -7,7 +7,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" -SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" +SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" [compat] @@ -19,7 +19,7 @@ Plots = "1" PyCall = "1" StaticArrays = "1" SymEngine = "0.8, 0.9.1, 0.10, 0.11" -SymPy = "1.0.53" +SymPyPythonCall = "0.3" Symbolics = "3, 4, 5" [preferences.OrdinaryDiffEq] diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md index 1d2c9aa4..b8072392 100644 --- a/docs/src/benchmarks.md +++ b/docs/src/benchmarks.md @@ -15,7 +15,8 @@ Symbolic computations of [`modified_equation`](@ref)s and support - [SymEngine.jl](https://github.com/symengine/SymEngine.jl), -- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and +- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), +- [SymPyPythonCall.jl](https://github.com/jverzani/SymPyPythonCall.jl), and - [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) as symbolic backends. Here, we compare them in the context of the explicit @@ -70,12 +71,12 @@ end Next, we load the symbolic packages and run the benchmarks. ```@setup benchmark-nonlinear-oscillator -using SymPy # generates annoying output online when conda installs sympy +using SymPyPythonCall # generates annoying output online when conda installs sympy ``` ```@example benchmark-nonlinear-oscillator using SymEngine: SymEngine -using SymPy: SymPy +using SymPyPythonCall: SymPyPythonCall using Symbolics: Symbolics println("SymEngine") @@ -85,9 +86,9 @@ subs = SymEngine.subs benchmark(u, dt, subs, 8) println("SymPy") -dt = SymPy.symbols("dt") -u = SymPy.symbols("u1, u2") -subs = SymPy.subs +dt = SymPyPythonCall.symbols("dt") +u = SymPyPythonCall.symbols("u1, u2") +subs = SymPyPythonCall.subs benchmark(u, dt, subs, 8) println("Symbolics") @@ -104,7 +105,7 @@ using InteractiveUtils versioninfo() using Pkg -Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPy", "Symbolics"], +Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPyPythonCall", "Symbolics"], mode=PKGMODE_MANIFEST) nothing # hide ``` diff --git a/docs/src/tutorials/bseries_basics.md b/docs/src/tutorials/bseries_basics.md index 91a0aceb..7ab21995 100644 --- a/docs/src/tutorials/bseries_basics.md +++ b/docs/src/tutorials/bseries_basics.md @@ -9,7 +9,7 @@ methods, generically or when applied to a specific ordinary differential equatio # These must first be installed using: import Pkg; Pkg.add("package_name") using BSeries using Latexify # Only needed for some pretty-printing cells below using `latexify` -import SymPy; sp=SymPy; +import SymPyPythonCall; sp = SymPyPythonCall; ``` ## B-series for a generic ODE @@ -20,7 +20,7 @@ Here is a generic 2-stage, 2nd-order method: ```@example bseries-basics -α = sp.symbols("α", real=true) +α = sp.symbols("α", real = true) A = [0 0; 1/(2*α) 0]; b = [1-α, α]; c = [0, 1/(2*α)] coeffs2 = bseries(A, b, c, 3) latexify(coeffs2, cdot=false) @@ -155,9 +155,9 @@ u_2'(t) & = t. ```@example bseries-basics -λ = sp.symbols("λ", real=true) -y, t = sp.symbols("y t", real=true) -h = sp.symbols("h", real=true) +λ = sp.symbols("λ", real = true) +y, t = sp.symbols("y t", real = true) +h = sp.symbols("h", real = true) u = [y, t] ff = [λ*(u[1]-sin(t))+cos(t), 1] diff --git a/docs/src/tutorials/rk_order_conditions.md b/docs/src/tutorials/rk_order_conditions.md index f81b7660..9b1c1a9f 100644 --- a/docs/src/tutorials/rk_order_conditions.md +++ b/docs/src/tutorials/rk_order_conditions.md @@ -3,11 +3,11 @@ In this tutorial, we generate order conditions for a generic explicit Runge-Kutta method. First, we create symbolic coefficient arrays with the appropriate structure. -There are several symbolic packages you can use. Here, we will use SymPy.jl. +There are several symbolic packages you can use. Here, we will use SymPyPythonCall.jl. ```@example bseries-RK-order-conditions -using BSeries, SymPy, Latexify +using BSeries, SymPyPythonCall, Latexify s = 4 # Stages p = 4 # Desired order of accuracy diff --git a/docs/src/tutorials/symbolic_computations.md b/docs/src/tutorials/symbolic_computations.md index fedf8ead..ac1c8944 100644 --- a/docs/src/tutorials/symbolic_computations.md +++ b/docs/src/tutorials/symbolic_computations.md @@ -6,7 +6,8 @@ computations in [BSeries.jl](https://github.com/ranocha/BSeries.jl) support at least - [SymEngine.jl](https://github.com/symengine/SymEngine.jl), -- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), and +- [SymPy.jl](https://github.com/JuliaPy/SymPy.jl), +- [SymPyPythonCall.jl](https://github.com/jverzani/SymPyPythonCall.jl), and - [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) as symbolic backends. You can find some performance comparisons of them in the @@ -87,12 +88,12 @@ latexify(series / h, cdot=false) |> println ``` We can also use other packages for the symbolic computations, of course. -SymPy.jl often provides very clean expressions. +SymPyPythonCall.jl often provides very clean expressions. ```@example modified-equation-sympy -using BSeries, SymPy +using BSeries, SymPyPythonCall -α = SymPy.symbols("α", real=true) +α = SymPyPythonCall.symbols("α", real = true) A = [0 0; 1/(2α) 0] b = [1-α, α] c = [0, 1/(2α)] @@ -105,20 +106,21 @@ SymEngine.jl. ```@example modified-equation-sympy using Latexify -latexify(series, reduce_order_by=1, dt=SymPy.symbols("h"), cdot=false) |> println +latexify(series, reduce_order_by = 1, dt = SymPyPythonCall.symbols("h"), + cdot = false) |> println ``` We can also use the simplified versions. ```@example modified-equation-sympy using Latexify -latexify(series, reduce_order_by=1, cdot=false) |> println +latexify(series, reduce_order_by = 1, cdot = false) |> println ``` ```@example modified-equation-sympy using Latexify -h = SymPy.symbols("h", real=true) -latexify(series / h, cdot=false) |> println +h = SymPyPythonCall.symbols("h", real = true) +latexify(series / h, cdot = false) |> println ``` Alternatively, we can also use Symbolics.jl. @@ -137,7 +139,7 @@ series = modified_equation(A, b, c, 3) ```@example modified-equation-symbolics using Latexify Symbolics.@variables h -latexify(series / h, cdot=false) |> println +latexify(series / h, cdot = false) |> println ``` @@ -163,14 +165,14 @@ If even `Int128` is not enough, one can specify the type `BigInt`, which has adj You can also create purely symbolic B-series as starting point of exploratory research, e.g., ```@example ex:symbolic-series -using BSeries, SymPy +using BSeries, SymPyPythonCall series = bseries(5) do t, series - return symbols("a_$(butcher_representation(t))", real=true) + return symbols("a_$(butcher_representation(t))", real = true) end ``` -This B-series can be used as any other B-series, e.g., to compute a +This B-series can be used as any other B-series, e.g., to compute a modified equation: ```@example ex:symbolic-series diff --git a/ext/SymEngineExt.jl b/ext/BSeriesSymEngineExt.jl similarity index 93% rename from ext/SymEngineExt.jl rename to ext/BSeriesSymEngineExt.jl index 3e9d5316..4d177dc8 100644 --- a/ext/SymEngineExt.jl +++ b/ext/BSeriesSymEngineExt.jl @@ -1,4 +1,4 @@ -module SymEngineExt +module BSeriesSymEngineExt if isdefined(Base, :get_extension) using SymEngine: SymEngine diff --git a/ext/BSeriesSymPyExt.jl b/ext/BSeriesSymPyExt.jl new file mode 100644 index 00000000..9987dea0 --- /dev/null +++ b/ext/BSeriesSymPyExt.jl @@ -0,0 +1,20 @@ +module BSeriesSymPyExt + +if isdefined(Base, :get_extension) + using SymPy: SymPy +else + using ..SymPy: SymPy +end + +using BSeries: BSeries + +function BSeries.compute_derivative(expression::SymPy.Sym{SymPy.PyCall.PyObject}, + variable::SymPy.Sym{SymPy.PyCall.PyObject}) + SymPy.diff(expression, variable) +end + +function BSeries.latexify_default_dt(::Type{SymPy.Sym{SymPy.PyCall.PyObject}}) + SymPy.symbols("h", real = true) +end + +end # module diff --git a/ext/BSeriesSymPyPythonCallExt.jl b/ext/BSeriesSymPyPythonCallExt.jl new file mode 100644 index 00000000..bbb25495 --- /dev/null +++ b/ext/BSeriesSymPyPythonCallExt.jl @@ -0,0 +1,20 @@ +module BSeriesSymPyPythonCallExt + +if isdefined(Base, :get_extension) + using SymPyPythonCall: SymPyPythonCall +else + using ..SymPyPythonCall: SymPyPythonCall +end + +using BSeries: BSeries + +function BSeries.compute_derivative(expression::SymPyPythonCall.Sym{SymPyPythonCall.PythonCall.Core.Py}, + variable::SymPyPythonCall.Sym{SymPyPythonCall.PythonCall.Core.Py}) + SymPyPythonCall.diff(expression, variable) +end + +function BSeries.latexify_default_dt(::Type{SymPyPythonCall.Sym{SymPyPythonCall.PythonCall.Core.Py}}) + SymPyPythonCall.symbols("h", real = true) +end + +end # module diff --git a/ext/SymbolicsExt.jl b/ext/BSeriesSymbolicsExt.jl similarity index 93% rename from ext/SymbolicsExt.jl rename to ext/BSeriesSymbolicsExt.jl index 46c9b22d..18eb320d 100644 --- a/ext/SymbolicsExt.jl +++ b/ext/BSeriesSymbolicsExt.jl @@ -1,4 +1,4 @@ -module SymbolicsExt +module BSeriesSymbolicsExt if isdefined(Base, :get_extension) using Symbolics: Symbolics diff --git a/ext/SymPyExt.jl b/ext/SymPyExt.jl deleted file mode 100644 index 2bf0b216..00000000 --- a/ext/SymPyExt.jl +++ /dev/null @@ -1,17 +0,0 @@ -module SymPyExt - -if isdefined(Base, :get_extension) - using SymPy: SymPy -else - using ..SymPy: SymPy -end - -using BSeries: BSeries - -function BSeries.compute_derivative(expression::SymPy.Sym, variable::SymPy.Sym) - SymPy.diff(expression, variable) -end - -BSeries.latexify_default_dt(::Type{SymPy.Sym}) = SymPy.symbols("h", real = true) - -end # module diff --git a/src/BSeries.jl b/src/BSeries.jl index 8fa023ea..9b1fd24e 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1605,15 +1605,19 @@ function compute_derivative end @static if !isdefined(Base, :get_extension) function __init__() @require Symbolics="0c5d862f-8b57-4792-8d23-62f2024744c7" begin - include("../ext/SymbolicsExt.jl") + include("../ext/BSeriesSymbolicsExt.jl") end @require SymEngine="123dc426-2d89-5057-bbad-38513e3affd8" begin - include("../ext/SymEngineExt.jl") + include("../ext/BSeriesSymEngineExt.jl") end @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin - include("../ext/SymPyExt.jl") + include("../ext/BSeriesSymPyExt.jl") + end + + @require SymPyPythonCall="bc8888f7-b21e-4b7c-a06a-5d9c9496438c" begin + include("../ext/BSeriesSymPyPythonCallExt.jl") end end end diff --git a/test/Project.toml b/test/Project.toml index e0c1884e..9de3479b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,9 +1,11 @@ [deps] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8" SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" +SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -11,5 +13,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Aqua = "0.8" StaticArrays = "1" SymEngine = "0.8, 0.9.1, 0.10, 0.11" -SymPy = "1" +SymPy = "2" +SymPyPythonCall = "0.3" Symbolics = "4, 5" diff --git a/test/runtests.jl b/test/runtests.jl index c6987a2e..c492a27a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,9 +6,16 @@ using BSeries.Latexify: latexify using LinearAlgebra: I using StaticArrays: @SArray, @SMatrix, @SVector +using Symbolics: Symbolics using SymEngine: SymEngine +# see +# https://juliapy.github.io/PythonCall.jl/stable/pycall/#Tips +# https://github.com/JuliaPy/PyCall.jl/issues/1056 +using SymPyPythonCall: SymPyPythonCall +ENV["PYTHON"] = SymPyPythonCall.PythonCall.python_executable_path() +import Pkg +Pkg.build("SymPy") using SymPy: SymPy -using Symbolics: Symbolics using Aqua: Aqua @@ -76,7 +83,7 @@ using Aqua: Aqua A = [0 0; 1/(2 * α) 0] b = [1 - α, α] c = [0, 1 / (2 * α)] - series_integrator = bseries(A, b, c, 3) + series_integrator = @inferred bseries(A, b, c, 3) # Call this once to avoid failing tests with `@test_nowarn` # caused by deprecation warnings from PyCall.jl. See also @@ -97,6 +104,44 @@ using Aqua: Aqua coefficients = @inferred modified_equation(series_integrator) val1 = @test_nowarn latexify(coefficients, reduce_order_by = 1, cdot = false) + @show typeof(h) + @test typeof(BSeries.latexify_default_dt(valtype(coefficients / h))) == + typeof(h) + val2 = @test_nowarn latexify(coefficients / h, cdot = false) + @test val1 == val2 + end + end + + @testset "SymPyPythonCall.jl" begin + @testset "Symbolic coefficients" begin + α = SymPyPythonCall.symbols("α", real = true) + A = [0 0; 1/(2 * α) 0] + b = [1 - α, α] + c = [0, 1 / (2 * α)] + series_integrator = @inferred bseries(A, b, c, 3) + + # Call this once to avoid failing tests with `@test_nowarn` + # caused by deprecation warnings from PyCall.jl. See also + # https://github.com/JuliaPy/PyCall.jl/pull/1042 + deepcopy(α) + + @test_nowarn latexify(series_integrator) + end + + @testset "Divide by h" begin + A = [0 0; 1//2 0] + b = [0, 1] + c = [0, 1 // 2] + series_integrator = @inferred bseries(A, b, c, 1) + @test_nowarn latexify(series_integrator) + + h = SymPyPythonCall.symbols("h", real = true) + coefficients = @inferred modified_equation(series_integrator) + val1 = @test_nowarn latexify(coefficients, reduce_order_by = 1, + cdot = false) + @show typeof(h) + @test typeof(BSeries.latexify_default_dt(valtype(coefficients / h))) == + typeof(h) val2 = @test_nowarn latexify(coefficients / h, cdot = false) @test val1 == val2 end @@ -839,6 +884,75 @@ using Aqua: Aqua @test mapreduce(iszero ∘ SymPy.expand, &, series4 - series4_reference) end + @testset "SymPyPythonCall.jl" begin + # Lotka-Volterra model + dt = SymPyPythonCall.symbols("dt") + p, q = u = SymPyPythonCall.symbols("p, q") + f = [p * (2 - q), q * (p - 1)] + + # Explicit Euler method + A = @SMatrix [0 // 1;;] + b = @SArray [1 // 1] + c = @SArray [0 // 1] + + # tested with the Python package BSeries + series2 = modified_equation(f, u, dt, A, b, c, 2) + series2_reference = [ + -dt * (-p * q * (p - 1) + p * (2 - q)^2) / 2 + p * (2 - q), + -dt * (p * q * (2 - q) + q * (p - 1)^2) / 2 + q * (p - 1), + ] + @test mapreduce(isequal, &, series2, series2_reference) + + # tested with the Python package BSeries + series3 = modified_equation(f, u, dt, A, b, c, 3) + series3_reference = [ + -dt^2 * p * q * (2 - q) * (p - 1) / 6 + + dt^2 * + (-p^2 * q * (2 - q) - p * q * (2 - q) * (p - 1) - p * q * (p - 1)^2 + + p * (2 - q)^3) / 3 - dt * (-p * q * (p - 1) + p * (2 - q)^2) / 2 + + p * (2 - q), + dt^2 * p * q * (2 - q) * (p - 1) / 6 + + dt^2 * + (-p * q^2 * (p - 1) + p * q * (2 - q)^2 + p * q * (2 - q) * (p - 1) + + q * (p - 1)^3) / 3 - dt * (p * q * (2 - q) + q * (p - 1)^2) / 2 + + q * (p - 1), + ] + @test mapreduce(iszero ∘ SymPy.expand, &, series3 - series3_reference) + + # tested with the Python package BSeries + series4 = modified_equation(f, u, dt, A, b, c, 4) + series4_reference = [ + -dt^3 * + (-2 * p^2 * q * (2 - q) * (p - 1) + 2 * p * q * (2 - q) * (p - 1) * (q - 2)) / + 12 - + dt^3 * (-p^2 * q * (2 - q)^2 + p * q^2 * (p - 1)^2 + + p * q * (1 - p) * (2 - q) * (p - 1) + p * q * (2 - q) * (p - 1) * (q - 2)) / + 12 - + dt^3 * (p^2 * q^2 * (p - 1) - 2 * p^2 * q * (2 - q)^2 - + p^2 * q * (2 - q) * (p - 1) - p * q * (2 - q)^2 * (p - 1) - + p * q * (2 - q) * (p - 1)^2 - p * q * (p - 1)^3 + p * (2 - q)^4) / 4 - + dt^2 * p * q * (2 - q) * (p - 1) / 6 + + dt^2 * + (-p^2 * q * (2 - q) - p * q * (2 - q) * (p - 1) - p * q * (p - 1)^2 + + p * (2 - q)^3) / 3 - dt * (-p * q * (p - 1) + p * (2 - q)^2) / 2 + + p * (2 - q), + -dt^3 * + (-2 * p * q^2 * (2 - q) * (p - 1) + 2 * p * q * (2 - q) * (p - 1)^2) / 12 - + dt^3 * + (p^2 * q * (2 - q)^2 - p * q^2 * (p - 1)^2 + p * q * (2 - q)^2 * (p - 1) + + p * q * (2 - q) * (p - 1)^2) / 12 - + dt^3 * (-p^2 * q^2 * (2 - q) - p * q^2 * (2 - q) * (p - 1) - + 2 * p * q^2 * (p - 1)^2 + p * q * (2 - q)^3 + p * q * (2 - q)^2 * (p - 1) + + p * q * (2 - q) * (p - 1)^2 + q * (p - 1)^4) / 4 + + dt^2 * p * q * (2 - q) * (p - 1) / 6 + + dt^2 * + (-p * q^2 * (p - 1) + p * q * (2 - q)^2 + p * q * (2 - q) * (p - 1) + + q * (p - 1)^3) / 3 - dt * (p * q * (2 - q) + q * (p - 1)^2) / 2 + + q * (p - 1), + ] + @test mapreduce(iszero ∘ SymPy.expand, &, series4 - series4_reference) + end + @testset "Symbolics.jl" begin # Lotka-Volterra model Symbolics.@variables dt @@ -1103,6 +1217,56 @@ using Aqua: Aqua end end + @testset "SymPyPythonCall.jl" begin + dt = SymPyPythonCall.symbols("dt") + α, β, γ = SymPyPythonCall.symbols("alpha, beta, gamma") + u1, u2, u3 = u = SymPyPythonCall.symbols("u1 u2 u3") + f = [α * u[2] * u[3], β * u[3] * u[1], γ * u[1] * u[2]] + + # Implicit midpoint method + A = @SMatrix [1 // 2;;] + b = @SArray [1 // 1] + c = @SArray [1 // 2] + + # See eq. (12) of + # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007) + # Numerical integrators based on modified differential equations + # [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9) + + series = modifying_integrator(f, u, dt, A, b, c, 5) + # Dirty workaround used also for the other symbolic setups - just make + # it consistent here, although we could use another approach with SymPyPythonCall.jl. + # We differentiate the expression and set `dt` to zero to get the corresponding + # coefficient divided by factorial(how often we needed to differentiate). + s3_reference = -(α * β * u3^2 + α * γ * u2^2 + β * γ * u1^2) / 12 + for i in eachindex(f) + s3 = SymPyPythonCall.subs(1 // 2 * + SymPyPythonCall.diff(SymPyPythonCall.diff((series[i] - + f[i]) / + f[i], + dt), + dt), + dt => 0) + @test iszero(SymPyPythonCall.expand(s3 - s3_reference)) + end + + s5_reference = 6 // 5 * s3_reference^2 + + 1 // 60 * α * β * γ * + (β * u1^2 * u3^2 + γ * u2^2 * u1^2 + α * u3^2 * u2^2) + for i in eachindex(f) + s5 = SymPyPythonCall.subs(1 // 24 * + SymPyPythonCall.diff(SymPyPythonCall.diff(SymPyPythonCall.diff(SymPyPythonCall.diff((series[i] - + f[i]) / + f[i], + dt), + dt), + dt), + dt), + dt => 0) + @test iszero(SymPyPythonCall.expand(s5 - s5_reference)) + end + end + @testset "Symbolics.jl" begin Symbolics.@variables dt Symbolics.@variables α β γ @@ -1646,7 +1810,7 @@ using Aqua: Aqua end # TODO: This is currently not implemented - @test_broken is_energy_preserving(series) + @test_skip is_energy_preserving(series) end @testset "SymPy.jl" begin @@ -1665,7 +1829,26 @@ using Aqua: Aqua @test @inferred(order_of_accuracy(series)) == 4 # TODO: This is currently not implemented - @test_broken is_energy_preserving(series) + @test_skip is_energy_preserving(series) + end + + @testset "SymPyPythonCall.jl" begin + # Examples in Section 5.3.1 + α = SymPyPythonCall.symbols("α", real = true) + α1 = 1 / (36 * α - 7) + M = [α1+4 -6 * α1-6 6*α1 + -6 * α1-6 36 * α1+12 -36*α1 + 6*α1 -36*α1 36*α1] + csrk = @inferred ContinuousStageRungeKuttaMethod(M) + + # TODO: This is no type stable at the moment + # series = @inferred bseries(csrk, 5) + series = bseries(csrk, 5) + + @test @inferred(order_of_accuracy(series)) == 4 + + # TODO: This is currently not implemented + @test_skip is_energy_preserving(series) end @testset "Symbolics.jl" begin @@ -1686,7 +1869,7 @@ using Aqua: Aqua # @test @inferred(order_of_accuracy(series)) == 4 # # TODO: This is currently not implemented - @test_broken is_energy_preserving(series) + @test_skip is_energy_preserving(series) end end @@ -2156,7 +2339,20 @@ using Aqua: Aqua series_integrator = @inferred(bseries(A, b, c, 2)) @test @inferred(order_of_accuracy(series_integrator)) == 2 # TODO: This test is currently broken and throws an error - @test_broken is_energy_preserving(series_integrator) + @test_skip is_energy_preserving(series_integrator) + 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, 2)) + @test @inferred(order_of_accuracy(series_integrator)) == 2 + # TODO: This test is currently broken and throws an error + @test_skip is_energy_preserving(series_integrator) end @testset "Symbolics.jl" begin