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_equation
s and modifying_integrator
s 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, 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 SymPy: SymPy
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 = SymPy.symbols("dt")
u = SymPy.symbols("u1, u2")
subs = SymPy.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.008967 seconds (242 allocations: 95.047 KiB)
0.009692 seconds (242 allocations: 95.047 KiB)
Computing the series including elementary differentials:
1.515342 seconds (1.05 M allocations: 68.152 MiB, 3.31% gc time, 97.84% compilation time)
0.028136 seconds (47.59 k allocations: 1.075 MiB)
Substituting the initial condition:
0.113300 seconds (124.11 k allocations: 7.858 MiB, 92.53% compilation time)
0.006343 seconds (10.24 k allocations: 119.344 KiB)
SymPy
Computing the series coefficients:
0.009056 seconds (242 allocations: 95.047 KiB)
0.009706 seconds (242 allocations: 95.047 KiB)
Computing the series including elementary differentials:
4.678071 seconds (751.99 k allocations: 48.134 MiB, 0.47% gc time, 25.47% compilation time)
2.631552 seconds (80.78 k allocations: 3.011 MiB)
Substituting the initial condition:
1.273507 seconds (319.43 k allocations: 20.959 MiB, 2.10% gc time, 23.41% compilation time)
0.963464 seconds (126 allocations: 4.141 KiB)
Symbolics
Computing the series coefficients:
0.009019 seconds (242 allocations: 95.047 KiB)
0.009103 seconds (242 allocations: 95.047 KiB)
Computing the series including elementary differentials:
13.401722 seconds (17.74 M allocations: 1.019 GiB, 4.29% gc time, 84.77% compilation time)
1.640169 seconds (6.24 M allocations: 299.055 MiB, 10.00% gc time)
Substituting the initial condition:
10.769903 seconds (22.33 M allocations: 840.703 MiB, 3.59% gc time, 43.08% compilation time)
5.990176 seconds (19.77 M allocations: 677.387 MiB, 3.45% gc time)
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["BSeries", "RootedTrees", "SymEngine", "SymPy", "Symbolics"],
mode=PKGMODE_MANIFEST)
Julia Version 1.9.3
Commit bed2cd540a1 (2023-08-24 14:43 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 2 × Intel(R) Xeon(R) CPU E5-2673 v4 @ 2.30GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-14.0.6 (ORCJIT, broadwell)
Threads: 1 on 2 virtual cores
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
LD_LIBRARY_PATH = /opt/hostedtoolcache/Python/3.9.18/x64/lib
Status `~/work/BSeries.jl/BSeries.jl/docs/Manifest.toml`
[ebb8d67c] BSeries v0.1.58 `~/work/BSeries.jl/BSeries.jl`
[47965b36] RootedTrees v2.19.2
[123dc426] SymEngine v0.11.0
[24249f21] SymPy v1.1.14
[0c5d862f] Symbolics v5.8.0
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.18 (main, Aug 28 2023, 08:38:32)
[GCC 11.4.0]
Package version 0.3
Modified equation
19063/26880
21.6888689994812 seconds
19063/26880
22.329633235931396 seconds
Modifying integrator
5460293/241920
20.179504871368408 seconds
5460293/241920
20.840371131896973 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")
The results are as follows.
Python version 3.9.18 (main, Aug 28 2023, 08:38:32)
[GCC 11.4.0]
Package version 0.3
Modified equation
19063/26880
9.608121395111084 seconds
19063/26880
9.22347617149353 seconds
Energy preservation
9
9.146897077560425 seconds
9
8.655707836151123 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.18 (main, Aug 28 2023, 08:38:32)
[GCC 11.4.0]
Modified equation
19063/26880
1.5211694240570068 seconds
19063/26880
1.4823462963104248 seconds
Modifying integrator
5460293/241920
1.262559175491333 seconds
5460293/241920
1.2672080993652344 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
Modified equation
19063//26880
0.301151 seconds (75.95 k allocations: 5.662 MiB, 66.02% compilation time)
19063//26880
0.103507 seconds (598 allocations: 367.750 KiB)
Modifying integrator
5460293//241920
0.053055 seconds (628 allocations: 334.273 KiB, 13.01% compilation time)
5460293//241920
0.046612 seconds (583 allocations: 332.633 KiB)
Energy preservation
true
3.274216 seconds (3.40 M allocations: 220.021 MiB, 5.99% gc time, 96.53% compilation time: 4% of which was recompilation)
true
0.106893 seconds (38.74 k allocations: 5.070 MiB)
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