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

• SymEngine.jl,
• SymPy.jl,
• SymPyPythonCall.jl, and
• Symbolics.jl

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, StaticArrays

function 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")
end

benchmark (generic function with 1 method)

Next, we load the symbolic packages and run the benchmarks.

using SymEngine: SymEngine
using SymPyPythonCall: SymPyPythonCall
using Symbolics: Symbolics

println("SymEngine")
dt   = SymEngine.symbols("dt")
u    = SymEngine.symbols("u1, u2")
subs = SymEngine.subs
benchmark(u, dt, subs, 8)

println("SymPy")
dt   = SymPyPythonCall.symbols("dt")
u    = SymPyPythonCall.symbols("u1, u2")
subs = SymPyPythonCall.subs
benchmark(u, dt, subs, 8)

println("Symbolics")
Symbolics.@variables dt
u = Symbolics.@variables u1 u2
subs = Symbolics.substitute
benchmark(u, dt, subs, 8)

SymEngine

 Computing the series coefficients:
  0.005025 seconds (470 allocations: 81.156 KiB)
  0.005043 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  1.072181 seconds (1.63 M allocations: 75.848 MiB, 98.55% compilation time)
  0.012635 seconds (50.95 k allocations: 1.306 MiB)

 Substituting the initial condition:
  0.090434 seconds (196.03 k allocations: 9.242 MiB, 94.75% compilation time)
  0.004070 seconds (10.24 k allocations: 239.422 KiB)


SymPy

 Computing the series coefficients:
  0.005024 seconds (470 allocations: 81.156 KiB)
  0.004980 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  3.352420 seconds (3.90 M allocations: 209.453 MiB, 1.58% gc time, 64.26% compilation time)
  0.812454 seconds (61.49 k allocations: 1.671 MiB, 3.98% compilation time)

 Substituting the initial condition:
  1.373046 seconds (613.54 k allocations: 30.951 MiB, 3.61% gc time, 53.93% compilation time)
  0.649574 seconds (142 allocations: 3.641 KiB)


Symbolics

 Computing the series coefficients:
  0.005037 seconds (470 allocations: 81.156 KiB)
  0.004974 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  4.015688 seconds (9.03 M allocations: 383.602 MiB, 1.31% gc time, 90.98% compilation time: 2% of which was recompilation)
  0.329197 seconds (3.34 M allocations: 98.931 MiB, 8.41% gc time, 0.26% compilation time)

 Substituting the initial condition:
  4.663091 seconds (20.95 M allocations: 710.624 MiB, 1.47% gc time, 71.23% compilation time)
  1.305344 seconds (17.02 M allocations: 526.948 MiB, 3.57% gc time)

These results were obtained using the following versions.

using InteractiveUtils
versioninfo()

using Pkg
Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPyPythonCall", "Symbolics"],
           mode=PKGMODE_MANIFEST)

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  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.71 `~/work/BSeries.jl/BSeries.jl`
  [47965b36] RootedTrees v2.24.1
  [123dc426] SymEngine v0.13.0
  [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

• BSeries
• pybs
• orderconditions

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
 15.052567958831787 seconds
19063/26880
 15.334388017654419 seconds

Modifying integrator
5460293/241920
 13.897258520126343 seconds
5460293/241920
 13.864242315292358 seconds

Python 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.735898733139038 seconds
19063/26880
 5.686330556869507 seconds

Energy preservation
9
 5.811593532562256 seconds
9
 5.644437074661255 seconds

Symplecticity (conservation of quadratic invariants)
9
 0.03485894203186035 seconds
9
 0.019673585891723633 seconds

Python 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.1016900539398193 seconds
19063/26880
 1.1030910015106201 seconds

Modifying integrator
5460293/241920
 0.9527595043182373 seconds
5460293/241920
 0.9525721073150635 seconds

This Julia package BSeries.jl

Finally, we perform the same task using BSeries.jl in Julia.

using BSeries, StaticArrays

A = @SArray [0 0; 1//2 0]
b = @SArray [0, 1//1]
c = @SArray [0, 1//2]
up_to_order = 9


println("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)))
end


println("\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)))
end


println("\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))
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

Modified equation
19063//26880
  0.146119 seconds (90.60 k allocations: 4.916 MiB, 63.05% compilation time)
19063//26880
  0.053850 seconds (1.10 k allocations: 451.148 KiB)

Modifying integrator
5460293//241920
  0.030499 seconds (1.12 k allocations: 424.953 KiB, 16.05% compilation time)
5460293//241920
  0.025544 seconds (1.08 k allocations: 422.359 KiB)

Energy preservation
true
  2.347073 seconds (4.28 M allocations: 212.087 MiB, 1.48% gc time, 96.96% compilation time)
true
  0.057068 seconds (64.27 k allocations: 4.866 MiB)

Symplecticity (conservation of quadratic invariants)
true
  0.432599 seconds (442.70 k allocations: 21.798 MiB, 99.76% compilation time)
true
  0.000728 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