Symbolic computations

This tutorial describes some possibilities for symbolic computations based on BSeries.jl. Currently, symbolic computations in BSeries.jl support at least

as symbolic backends. You can find some performance comparisons of them in the benchmarks.

Generating LaTeX code

BSeries.jl allow fully symbolic computations of modified_equations and modifying_integrators. More information about these features are available in the respective tutorials. Here, we will also demonstrate how LaTeX code can be generated via Latexify.jl.

At the core, BSeries.jl is based on the representation of rooted trees via RootedTrees.jl. Thus, you need to use the setup described in the docstring of RootedTrees.latexify.

LaTeX code of a general modified equation

An explicit two-stage, second-order Runge-Kutta method is given by the Butcher tableau

\[\begin{array}{c|cc} 0 & 0 & \\ \frac{1}{2\alpha} & \frac{1}{2\alpha} & 0 \\ \hline & 1 - \alpha & \alpha\\ \end{array}\]

with real parameter $\alpha$. We can compute the B-series of its modified equation symbolically as follows.

using BSeries, SymEngine

α = SymEngine.symbols("α")
A = [0 0; 1/(2α) 0]
b = [1-α, α]
c = [0, 1/(2α)]

series = modified_equation(A, b, c, 3)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymEngine.Basic} with 5 entries:
  RootedTree{Int64}: Int64[]   => 0
  RootedTree{Int64}: [1]       => 1
  RootedTree{Int64}: [1, 2]    => 0
  RootedTree{Int64}: [1, 2, 3] => -1/6
  RootedTree{Int64}: [1, 2, 2] => -1/3 + (1/4)*α^(-1)

To generate the appropriate LaTeX code, we need to remember that the B-series coefficients shown above represent the modified vector field multiplied by the time step size. Thus, we need to pass reduce_order_by = 1 to latexify. Other keyword arguments implemented for B-series are

  • f = "f": The symbol used for the indices of the elementary differentials.
  • dt = "h": The symbol used for the time step size. Alternatively, symbolic variables can also be used.
  • reduce_order_by = 0: Reduce the power of the time step size accordingly.
using Latexify
latexify(series, reduce_order_by=1, dt=SymEngine.symbols("h"), cdot=false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{-1}{6} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{2} h^{2} \left( \frac{-1}{3} + \frac{1}{4} \alpha^{-1} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

The argument dt=SymEngine.symbols("h") is not required here.

using Latexify
latexify(series, reduce_order_by=1, cdot=false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{-1}{6} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{2} h^{2} \left( \frac{-1}{3} + \frac{1}{4} \alpha^{-1} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

We can obtain basically the same result (up to other notations of the coefficients) by dividing the B-series series by a symbolic variable for the time step size as follows.

using Latexify
h = SymEngine.symbols("h")
latexify(series / h, cdot=false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{-1}{6} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{2} h^{2} \left( \frac{-1}{3} + \frac{1}{4} \alpha^{-1} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

We can also use other packages for the symbolic computations, of course. SymPyPythonCall.jl often provides very clean expressions.

using BSeries, SymPyPythonCall

α = SymPyPythonCall.symbols("α", real = true)
A = [0 0; 1/(2α) 0]
b = [1-α, α]
c = [0, 1/(2α)]

series = modified_equation(A, b, c, 3)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPyCore.Sym{PythonCall.Core.Py}} with 5 entries:
  RootedTree{Int64}: Int64[]   => 0
  RootedTree{Int64}: [1]       => 1
  RootedTree{Int64}: [1, 2]    => 0
  RootedTree{Int64}: [1, 2, 3] => -1/6
  RootedTree{Int64}: [1, 2, 2] => -1/3 + 1/(4*α)

We can also generate LaTeX code as follows, using the same approach as for SymEngine.jl.

using Latexify
latexify(series, reduce_order_by = 1, dt = SymPyPythonCall.symbols("h"),
         cdot = false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{ - h^{2}}{6} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + h^{2} \left( \frac{-1}{6} + \frac{1}{8 \alpha} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

We can also use the simplified versions.

using Latexify
latexify(series, reduce_order_by = 1, cdot = false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{ - h^{2}}{6} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + h^{2} \left( \frac{-1}{6} + \frac{1}{8 \alpha} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$
using Latexify
h = SymPyPythonCall.symbols("h", real = true)
latexify(series / h, cdot = false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} + \frac{ - h^{2}}{6} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{h^{2} \left( \frac{-1}{3} + \frac{1}{4 \alpha} \right)}{2} F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

Alternatively, we can also use Symbolics.jl.

using BSeries, Symbolics

Symbolics.@variables α
A = [0 0; 1/(2α) 0]
b = [1-α, α]
c = [0, 1/(2α)]

series = modified_equation(A, b, c, 3)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Symbolics.Num} with 5 entries:
  RootedTree{Int64}: Int64[]   => 0
  RootedTree{Int64}: [1]       => 1
  RootedTree{Int64}: [1, 2]    => 0//1
  RootedTree{Int64}: [1, 2, 3] => -1//6
  RootedTree{Int64}: [1, 2, 2] => -(1//3) + 1 / (4α)
using Latexify
Symbolics.@variables h
latexify(series / h, cdot = false) |> println
$F_{f}\mathopen{}\left( \rootedtree[.] \right)\mathclose{} - \frac{1}{6} h^{2} F_{f}\mathopen{}\left( \rootedtree[.[.[.]]] \right)\mathclose{} + \frac{1}{2} h^{2} \left( \frac{-1}{3} + \frac{1}{4 \alpha} \right) F_{f}\mathopen{}\left( \rootedtree[.[.][.]] \right)\mathclose{}$

Working with rational coefficients

High-order Runge-Kutta methods are often given in terms of rational coefficients that have large denominators. It's often best to use the rational coefficients in calculations, rather than converting to floating-point, since otherwise terms that ought to cancel may not cancel exactly. Rational coefficients can be entered by using // for the division operation.

By default, Julia uses 64-bit integers for the numerator and denominator of a rational on 64-bit systems. In practical calculations with high-order RK methods, the denominators may become too large to be represented with 64 bits, leading to overflow. This can be remedied by specifying higher precision when entering the coefficients:

b = Rational{Int128}[(1081252805//134140608),(2639189439//74522560),(33646441//4191894),(-7873511875//210792384),(-504040617//14904512),(2110843561//115277085),(13//7),(1//2)]
8-element Vector{Rational{Int128}}:
  1081252805//134140608
  2639189439//74522560
    33646441//4191894
 -7873511875//210792384
  -504040617//14904512
  2110843561//115277085
          13//7
           1//2

If even Int128 is not enough, one can specify the type BigInt, which has adjustable precision.

Setting up symbolic B-series

You can also create purely symbolic B-series as starting point of exploratory research, e.g.,

using BSeries, SymPyPythonCall

series = bseries(5) do t, series
    return symbols("a_$(butcher_representation(t))", real = true)
end
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPyCore.Sym{PythonCall.Core.Py}} with 18 entries:
  RootedTree{Int64}: Int64[]         => a_∅
  RootedTree{Int64}: [1]             => a_τ
  RootedTree{Int64}: [1, 2]          => a_[τ]
  RootedTree{Int64}: [1, 2, 3]       => a_[[τ]]
  RootedTree{Int64}: [1, 2, 2]       => a_[τ²]
  RootedTree{Int64}: [1, 2, 3, 4]    => a_[[[τ]]]
  RootedTree{Int64}: [1, 2, 3, 3]    => a_[[τ²]]
  RootedTree{Int64}: [1, 2, 3, 2]    => a_[[τ]τ]
  RootedTree{Int64}: [1, 2, 2, 2]    => a_[τ³]
  RootedTree{Int64}: [1, 2, 3, 4, 5] => a_[[[[τ]]]]
  RootedTree{Int64}: [1, 2, 3, 4, 4] => a_[[[τ²]]]
  RootedTree{Int64}: [1, 2, 3, 4, 3] => a_[[[τ]τ]]
  RootedTree{Int64}: [1, 2, 3, 4, 2] => a_[[[τ]]τ]
  RootedTree{Int64}: [1, 2, 3, 3, 3] => a_[[τ³]]
  RootedTree{Int64}: [1, 2, 3, 3, 2] => a_[[τ²]τ]
  RootedTree{Int64}: [1, 2, 3, 2, 3] => a_[[τ][τ]]
  RootedTree{Int64}: [1, 2, 3, 2, 2] => a_[[τ]τ²]
  RootedTree{Int64}: [1, 2, 2, 2, 2] => a_[τ⁴]

This B-series can be used as any other B-series, e.g., to compute a modified equation:

modified_equation(series)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPyCore.Sym{PythonCall.Core.Py}} with 18 entries:
  RootedTree{Int64}: Int64[]         => 0
  RootedTree{Int64}: [1]             => a_τ
  RootedTree{Int64}: [1, 2]          => a_[τ] - a_τ^2/2
  RootedTree{Int64}: [1, 2, 3]       => a_[[τ]] - a_τ^3/6 - a_τ*(a_[τ]/2 - a_τ^…
  RootedTree{Int64}: [1, 2, 2]       => a_[τ²] - a_τ^3/3 - a_τ*(a_[τ] - a_τ^2/2)
  RootedTree{Int64}: [1, 2, 3, 4]    => a_[[[τ]]] - a_τ^4/24 - a_τ^2*(a_[τ]/6 -…
  RootedTree{Int64}: [1, 2, 3, 3]    => a_[[τ²]] - a_τ^4/12 - 2*a_τ^2*(a_[τ] - …
  RootedTree{Int64}: [1, 2, 3, 2]    => a_[[τ]τ] - a_τ^4/8 - 5*a_τ^2*(a_[τ] - a…
  RootedTree{Int64}: [1, 2, 2, 2]    => a_[τ³] - a_τ^4/4 - a_τ^2*(a_[τ] - a_τ^2…
  RootedTree{Int64}: [1, 2, 3, 4, 5] => a_[[[[τ]]]] - a_τ^5/120 - a_τ^3*(a_[τ]/…
  RootedTree{Int64}: [1, 2, 3, 4, 4] => a_[[[τ²]]] - a_τ^5/60 - a_τ^3*(a_[τ] - …
  RootedTree{Int64}: [1, 2, 3, 4, 3] => a_[[[τ]τ]] - a_τ^5/40 - a_τ^3*(a_[τ] - …
  RootedTree{Int64}: [1, 2, 3, 4, 2] => a_[[[τ]]τ] - a_τ^5/30 - 5*a_τ^3*(a_[τ] …
  RootedTree{Int64}: [1, 2, 3, 3, 3] => a_[[τ³]] - a_τ^5/20 - a_τ^3*(a_[τ] - a_…
  RootedTree{Int64}: [1, 2, 3, 3, 2] => a_[[τ²]τ] - a_τ^5/15 - 7*a_τ^3*(a_[τ] -…
  RootedTree{Int64}: [1, 2, 3, 2, 3] => a_[[τ][τ]] - a_τ^5/20 - a_τ^3*(a_[τ]/8 …
  RootedTree{Int64}: [1, 2, 3, 2, 2] => a_[[τ]τ²] - a_τ^5/10 - 3*a_τ^3*(a_[τ] -…
  RootedTree{Int64}: [1, 2, 2, 2, 2] => a_[τ⁴] - a_τ^5/5 - a_τ^3*(a_[τ] - a_τ^2…