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.005020 seconds (470 allocations: 81.156 KiB)
  0.005025 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  1.078929 seconds (1.63 M allocations: 75.788 MiB, 3.55% gc time, 98.59% compilation time)
  0.012520 seconds (50.95 k allocations: 1.306 MiB)

 Substituting the initial condition:
  0.088976 seconds (196.03 k allocations: 9.244 MiB, 94.79% compilation time)
  0.003901 seconds (10.24 k allocations: 238.891 KiB)


SymPy

 Computing the series coefficients:
  0.005068 seconds (470 allocations: 81.156 KiB)
  0.005016 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  3.273391 seconds (3.90 M allocations: 209.576 MiB, 1.33% gc time, 63.70% compilation time)
  0.806500 seconds (61.49 k allocations: 1.671 MiB, 3.96% compilation time)

 Substituting the initial condition:
  1.302128 seconds (613.55 k allocations: 30.954 MiB, 51.44% compilation time)
  0.629089 seconds (142 allocations: 3.641 KiB)


Symbolics

 Computing the series coefficients:
  0.005066 seconds (470 allocations: 81.156 KiB)
  0.005008 seconds (470 allocations: 81.156 KiB)

 Computing the series including elementary differentials:
  3.959271 seconds (9.08 M allocations: 385.428 MiB, 0.69% gc time, 91.42% compilation time: 2% of which was recompilation)
  0.330285 seconds (3.53 M allocations: 104.589 MiB, 7.56% gc time)

 Substituting the initial condition:
  4.590427 seconds (20.95 M allocations: 710.623 MiB, 1.91% gc time, 70.72% compilation time)
  1.286335 seconds (17.02 M allocations: 526.949 MiB, 3.60% 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.1
Commit ba1e628ee49 (2025-10-17 13:02 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.70 `~/work/BSeries.jl/BSeries.jl`
  [47965b36] RootedTrees v2.24.0
  [123dc426] SymEngine v0.13.0
  [bc8888f7] SymPyPythonCall v0.5.1
⌅ [0c5d862f] Symbolics v6.57.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
 14.673510074615479 seconds
19063/26880
 15.11969518661499 seconds

Modifying integrator
5460293/241920
 12.970076322555542 seconds
5460293/241920
 13.045676708221436 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.642908573150635 seconds
19063/26880
 5.582064151763916 seconds

Energy preservation
9
 5.661517858505249 seconds
9
 5.554283857345581 seconds

Symplecticity (conservation of quadratic invariants)
9
 0.03452801704406738 seconds
9
 0.019342899322509766 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.0750401020050049 seconds
19063/26880
 1.0734894275665283 seconds

Modifying integrator
5460293/241920
 0.9333069324493408 seconds
5460293/241920
 0.9260282516479492 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.145002 seconds (90.60 k allocations: 4.928 MiB, 62.39% compilation time)
19063//26880
  0.054307 seconds (1.10 k allocations: 450.867 KiB)

Modifying integrator
5460293//241920
  0.030659 seconds (1.12 k allocations: 424.672 KiB, 15.95% compilation time)
5460293//241920
  0.025980 seconds (1.08 k allocations: 422.641 KiB)

Energy preservation
true
  2.300955 seconds (4.29 M allocations: 212.229 MiB, 2.94% gc time, 97.47% compilation time)
true
  0.057163 seconds (64.27 k allocations: 4.866 MiB)

Symplecticity (conservation of quadratic invariants)
true
  0.407041 seconds (442.70 k allocations: 21.797 MiB, 99.75% compilation time)
true
  0.000737 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