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. SymPy.jl often provides very clean expressions.
using BSeries, SymPy
α = SymPy.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}}, SymPy.Sym} 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=SymPy.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 = SymPy.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//2If 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, SymPy
series = bseries(5) do t, series
return symbols("a_$(butcher_representation(t))", real=true)
endTruncatedBSeries{RootedTree{Int64, Vector{Int64}}, SymPy.Sym} 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}}, SymPy.Sym} 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…