Benchmarks
Here, we collect some simple benchmarks of BSeries.jl. Take them with a grain of salt since they run on virtual machines in the cloud to generate the documentation automatically. You can of course also copy the code and run the benchmarks locally yourself.
Comparing different symbolic packages
Symbolic computations of modified_equations and modifying_integrators in BSeries.jl support
as symbolic backends. Here, we compare them in the context of the explicit midpoint method and the nonlinear oscillator ODE
\[u'(t) = \frac{1}{\| u(t) \|^2} \begin{pmatrix} -u_2(t) \\ u_1(t) \end{pmatrix}.\]
This particular combination of explicit Runge-Kutta method and ODE is special since the explicit midpoint method is unconditionally energy-conserving for this problem[RanochaKetcheson2020].
First, we set up some code to perform the benchmarks. Here, we use a very naive approach, run the code twice (to see the effect of compilation) and use @time to print the runtime. More sophisticated approaches should of course use something like @benchmark from BenchmarkTools.jl. However, this simple and cheap version suffices to compare the orders of magnitude.
using BSeries, StaticArraysfunction benchmark(u, dt, subs, order) # explicit midpoint method A = @SArray [0 0; 1//2 0] b = @SArray [0, 1//1] c = @SArray [0, 1//2] # nonlinear oscillator f = [-u[2], u[1]] / (u[1]^2 + u[2]^2) println("\n Computing the series coefficients:") @time coefficients = modifying_integrator(A, b, c, order) @time coefficients = modifying_integrator(A, b, c, order) println("\n Computing the series including elementary differentials:") @time series = modifying_integrator(f, u, dt, A, b, c, order) @time series = modifying_integrator(f, u, dt, A, b, c, order) substitution_variables = Dict(u[1] => 1//1, u[2] => 0//1) println("\n Substituting the initial condition:") @time subs.(series, (substitution_variables, )) @time subs.(series, (substitution_variables, )) println("\n")endbenchmark (generic function with 1 method)Next, we load the symbolic packages and run the benchmarks.
using SymEngine: SymEngineusing SymPyPythonCall: SymPyPythonCallusing Symbolics: Symbolicsprintln("SymEngine")dt = SymEngine.symbols("dt")u = SymEngine.symbols("u1, u2")subs = SymEngine.subsbenchmark(u, dt, subs, 8)println("SymPy")dt = SymPyPythonCall.symbols("dt")u = SymPyPythonCall.symbols("u1, u2")subs = SymPyPythonCall.subsbenchmark(u, dt, subs, 8)println("Symbolics")Symbolics.@variables dtu = Symbolics.@variables u1 u2subs = Symbolics.substitutebenchmark(u, dt, subs, 8)SymEngine
Computing the series coefficients:
0.005102 seconds (470 allocations: 81.156 KiB)
0.005085 seconds (470 allocations: 81.156 KiB)
Computing the series including elementary differentials:
1.151186 seconds (1.62 M allocations: 75.311 MiB, 4.12% gc time, 98.42% compilation time)
0.013488 seconds (50.95 k allocations: 1.307 MiB)
Substituting the initial condition:
0.093762 seconds (196.02 k allocations: 9.246 MiB, 94.25% compilation time)
0.004389 seconds (10.24 k allocations: 239.188 KiB)
SymPy
Computing the series coefficients:
0.005141 seconds (470 allocations: 81.156 KiB)
0.005032 seconds (470 allocations: 81.156 KiB)
Computing the series including elementary differentials:
3.124516 seconds (3.89 M allocations: 208.755 MiB, 67.13% compilation time)
0.674193 seconds (53.68 k allocations: 1.247 MiB)
Substituting the initial condition:
1.224201 seconds (619.53 k allocations: 31.229 MiB, 55.37% compilation time)
0.550444 seconds (162 allocations: 4.016 KiB)
Symbolics
Computing the series coefficients:
0.005114 seconds (470 allocations: 81.156 KiB)
0.005086 seconds (470 allocations: 81.156 KiB)
Computing the series including elementary differentials:
4.289571 seconds (8.98 M allocations: 380.199 MiB, 2.18% gc time, 90.26% compilation time: 2% of which was recompilation)
0.292043 seconds (3.36 M allocations: 99.201 MiB, 0.32% compilation time)
Substituting the initial condition:
5.056002 seconds (20.97 M allocations: 711.310 MiB, 2.39% gc time, 72.09% compilation time)
1.399053 seconds (17.02 M allocations: 526.948 MiB, 5.62% gc time)These results were obtained using the following versions.
using InteractiveUtilsversioninfo()using PkgPkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPyPythonCall", "Symbolics"], mode=PKGMODE_MANIFEST)Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Build Info:
Official https://julialang.org release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 9V74 80-Core Processor
WORD_SIZE: 64
LLVM: libLLVM-18.1.7 (ORCJIT, znver4)
GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.9.25/x64/lib
JULIA_PYTHONCALL_EXE = /home/runner/work/BSeries.jl/BSeries.jl/docs/.CondaPkg/.pixi/envs/default/bin/python
Status `~/work/BSeries.jl/BSeries.jl/docs/Manifest.toml`
[ebb8d67c] BSeries v0.1.72 `~/work/BSeries.jl/BSeries.jl`
[47965b36] RootedTrees v2.25.4
[123dc426] SymEngine v0.13.2
[bc8888f7] SymPyPythonCall v0.5.1
⌅ [0c5d862f] Symbolics v6.58.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`Comparison with other packages
There are also other open source packages for B-series. Currently, we are aware of the Python packages
If you know about similar open source packages out there, please inform us, e.g., by creating an issue on GitHub.
The packages listed above and BSeries.jl all use different approaches and have different features. Thus, comparisons must be restricted to their common subset of features. Here, we present some simple performance comparisons. Again, we just use (the equivalent of) @time twice to get an idea of the performance after compilation, allowing us to compare orders of magnitude.
Python package BSeries
First, we start with the Python package BSeries and the following benchmark script.
import sys
from importlib.metadata import version
print("Python version", sys.version)
print("Package version", version('pybs'))
import time
import BSeries.bs as bs
import nodepy.runge_kutta_method as rk
midpoint_method = rk.loadRKM("Mid22")
up_to_order = 9
print("\nModified equation")
start_time = time.time()
series = bs.modified_equation(None, None,
midpoint_method.A, midpoint_method.b,
up_to_order, True)
result = sum(series.values())
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
series = bs.modified_equation(None, None,
midpoint_method.A, midpoint_method.b,
up_to_order, True)
result = sum(series.values())
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
print("\nModifying integrator")
start_time = time.time()
series = bs.modifying_integrator(None, None,
midpoint_method.A, midpoint_method.b,
up_to_order, True)
result = sum(series.values())
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
series = bs.modifying_integrator(None, None,
midpoint_method.A, midpoint_method.b,
up_to_order, True)
result = sum(series.values())
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")The results are as follows.
Python version 3.9.25 (main, Nov 3 2025, 15:16:36)
[GCC 13.3.0]
Package version 0.3
Modified equation
19063/26880
17.076056241989136 seconds
19063/26880
17.922030210494995 seconds
Modifying integrator
5460293/241920
15.959421634674072 seconds
5460293/241920
16.22330093383789 secondsPython package pybs
Next, we look at the Python package pybs and the following benchmark script. Note that this package does not provide functionality for modifying integrators.
import sys
from importlib.metadata import version
print("Python version", sys.version)
print("Package version", version('pybs'))
import time
import pybs
from pybs.rungekutta import methods as rk_methods
midpoint_method = rk_methods.RKmidpoint
up_to_order = 9
number_of_terms = pybs.unordered_tree.number_of_trees_up_to_order(up_to_order+1)
from itertools import islice
def first_values(f, n):
return (f(tree) for tree in islice(pybs.unordered_tree.tree_generator(), 0, n))
print("\nModified equation")
start_time = time.time()
midpoint_series = midpoint_method.phi()
series = pybs.series.modified_equation(midpoint_series)
result = sum(first_values(series, number_of_terms))
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
midpoint_series = midpoint_method.phi()
series = pybs.series.modified_equation(midpoint_series)
result = sum(first_values(series, number_of_terms))
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
print("\nEnergy preservation")
start_time = time.time()
a = pybs.series.AVF
b = pybs.series.modified_equation(a)
result = pybs.series.energy_preserving_upto_order(b, up_to_order)
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
a = pybs.series.AVF
b = pybs.series.modified_equation(a)
result = pybs.series.energy_preserving_upto_order(b, up_to_order)
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
print("\nSymplecticity (conservation of quadratic invariants)")
start_time = time.time()
a = rk_methods.RKimplicitMidpoint.phi()
result = pybs.series.symplectic_up_to_order(a, up_to_order)
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
a = rk_methods.RKimplicitMidpoint.phi()
result = pybs.series.symplectic_up_to_order(a, up_to_order)
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")The results are as follows.
Python version 3.9.25 (main, Nov 3 2025, 15:16:36)
[GCC 13.3.0]
Package version 0.3
Modified equation
19063/26880
5.36523699760437 seconds
19063/26880
5.335699558258057 seconds
Energy preservation
9
5.39194130897522 seconds
9
5.292031288146973 seconds
Symplecticity (conservation of quadratic invariants)
9
0.03091144561767578 seconds
9
0.018796205520629883 secondsPython package orderconditions
Next, we look at the Python package orderconditions of Valentin Dallerit and the following benchmark script.
import sys
print("Python version", sys.version)
import time
from orderConditions import BSeries
import nodepy.runge_kutta_method as rk
midpoint_method = rk.loadRKM("Mid22")
up_to_order = 9
print("\nModified equation")
start_time = time.time()
BSeries.set_order(up_to_order)
Y1 = BSeries.y() + midpoint_method.A[1,0] * BSeries.hf()
rk2 = BSeries.y() + midpoint_method.b[0] * BSeries.hf() + midpoint_method.b[1] * BSeries.compo_hf(Y1)
series = BSeries.modified_equation(rk2)
result = series.sum()
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
BSeries.set_order(up_to_order)
Y1 = BSeries.y() + midpoint_method.A[1,0] * BSeries.hf()
rk2 = BSeries.y() + midpoint_method.b[0] * BSeries.hf() + midpoint_method.b[1] * BSeries.compo_hf(Y1)
series = BSeries.modified_equation(rk2)
result = series.sum()
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
print("\nModifying integrator")
start_time = time.time()
BSeries.set_order(up_to_order)
Y1 = BSeries.y() + midpoint_method.A[1,0] * BSeries.hf()
rk2 = BSeries.y() + midpoint_method.b[0] * BSeries.hf() + midpoint_method.b[1] * BSeries.compo_hf(Y1)
series = BSeries.modifying_integrator(rk2)
result = series.sum()
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")
start_time = time.time()
BSeries.set_order(up_to_order)
Y1 = BSeries.y() + midpoint_method.A[1,0] * BSeries.hf()
rk2 = BSeries.y() + midpoint_method.b[0] * BSeries.hf() + midpoint_method.b[1] * BSeries.compo_hf(Y1)
series = BSeries.modifying_integrator(rk2)
result = series.sum()
end_time = time.time()
print(result)
print("", end_time - start_time, "seconds")The results are as follows.
Python version 3.9.25 (main, Nov 3 2025, 15:16:36)
[GCC 13.3.0]
Modified equation
19063/26880
1.0805635452270508 seconds
19063/26880
1.0779643058776855 seconds
Modifying integrator
5460293/241920
0.9396345615386963 seconds
5460293/241920
0.9289748668670654 secondsThis Julia package BSeries.jl
Finally, we perform the same task using BSeries.jl in Julia.
using BSeries, StaticArraysA = @SArray [0 0; 1//2 0]b = @SArray [0, 1//1]c = @SArray [0, 1//2]up_to_order = 9println("Modified equation")@time begin series = modified_equation(A, b, c, up_to_order) println(sum(values(series)))end@time begin series = modified_equation(A, b, c, up_to_order) println(sum(values(series)))endprintln("\nModifying integrator")@time begin series = modifying_integrator(A, b, c, up_to_order) println(sum(values(series)))end@time begin series = modifying_integrator(A, b, c, up_to_order) println(sum(values(series)))endprintln("\nEnergy preservation")@time begin series = bseries(AverageVectorFieldMethod(), up_to_order) println(is_energy_preserving(series))end@time begin series = bseries(AverageVectorFieldMethod(), up_to_order) println(is_energy_preserving(series))endprintln("\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))endModified equation
19063//26880
0.159033 seconds (90.20 k allocations: 4.895 MiB, 63.84% compilation time)
19063//26880
0.056355 seconds (1.10 k allocations: 450.867 KiB)
Modifying integrator
5460293//241920
0.031082 seconds (1.12 k allocations: 424.672 KiB, 17.93% compilation time)
5460293//241920
0.025630 seconds (1.08 k allocations: 422.359 KiB)
Energy preservation
true
2.880952 seconds (4.38 M allocations: 216.478 MiB, 15.95% gc time, 97.88% compilation time)
true
0.060262 seconds (64.27 k allocations: 4.867 MiB)
Symplecticity (conservation of quadratic invariants)
true
0.433055 seconds (442.73 k allocations: 21.805 MiB, 99.73% compilation time)
true
0.000791 seconds (1.85 k allocations: 238.867 KiB)References
Hendrik Ranocha and David Ketcheson (2020) Energy Stability of Explicit Runge-Kutta Methods for Nonautonomous or Nonlinear Problems. SIAM Journal on Numerical Analysis DOI: 10.1137/19M1290346