Code generation

This tutorial shows how to generate C code to compute expressions found using BSeries.jl. Although BSeries.jl is compatible with three symbolic backends, it's currently easiest to perform code generation using Symbolics.jl.

First, we generate the B-series that we want work with. Here we take a generic 2nd-order RK method and generate terms only up to 3rd order, in order to work with the leading truncation error.

using BSeries, Latexify, Symbolics

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

rk22 = bseries(A, b, c, 3)
exact = ExactSolution(rk22)
truncation_error = rk22 - exact

TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Symbolics.Num} with 5 entries:
  RootedTree{Int64}: Int64[]   => 0
  RootedTree{Int64}: [1]       => 0//1
  RootedTree{Int64}: [1, 2]    => 0//1
  RootedTree{Int64}: [1, 2, 3] => -1//6
  RootedTree{Int64}: [1, 2, 2] => -(1//3) + 1 / (4α)

Next we set up the ODE of interest, and evaluate the B-series with that right-hand side.

@variables h # time step size
@variables p q # variables of the ODE
f = [p * (2 - q), q * (p - 1)]

du = evaluate(f, [p, q], h, truncation_error)

\[ \begin{equation} \left[ \begin{array}{c} - \frac{1}{6} h^{3} \left( - p \left( \left( -1 + p \right)^{2} q + p q \left( 2 - q \right) \right) + \left( - p \left( -1 + p \right) q + \left( 2 - q \right)^{2} p \right) \left( 2 - q \right) \right) - h^{3} p \left( -1 + p \right) q \left( 2 - q \right) \left( \frac{-1}{3} + \frac{1}{4 \alpha} \right) \\ - \frac{1}{6} h^{3} \left( \left( -1 + p \right) \left( \left( -1 + p \right)^{2} q + p q \left( 2 - q \right) \right) + \left( - p \left( -1 + p \right) q + \left( 2 - q \right)^{2} p \right) q \right) + h^{3} p \left( -1 + p \right) q \left( 2 - q \right) \left( \frac{-1}{3} + \frac{1}{4 \alpha} \right) \\ \end{array} \right] \end{equation} \]

Finally, we generate a C function that evaluates the expressions above.

println(build_function(du, α, p, q, h, target = Symbolics.CTarget()))

#include <math.h>
void diffeqf(double* du, const double RHS1, const double RHS2, const double RHS3, const double RHS4) {
  du[0] = -0.16666666666666666 * (RHS4 * RHS4 * RHS4) * (-1 * RHS2 * (pow(-1 + RHS2, 2) * RHS3 + RHS2 * RHS3 * (2 + -1 * RHS3 * 1)) + (-1 * RHS2 * (-1 + RHS2) * RHS3 + RHS2 * pow(2 + -1 * RHS3 * 1, 2)) * (2 + -1 * RHS3 * 1)) + -1.0 * (RHS4 * RHS4 * RHS4) * RHS2 * (-1 + RHS2) * RHS3 * (2 + -1 * RHS3 * 1) * (-0.3333333333333333 + 1 / (4 * RHS1 * 1));
  du[1] = -0.16666666666666666 * (RHS4 * RHS4 * RHS4) * ((-1 + RHS2) * (pow(-1 + RHS2, 2) * RHS3 + RHS2 * RHS3 * (2 + -1 * RHS3 * 1)) + (-1 * RHS2 * (-1 + RHS2) * RHS3 + RHS2 * pow(2 + -1 * RHS3 * 1, 2)) * RHS3) + (RHS4 * RHS4 * RHS4) * RHS2 * (-1 + RHS2) * RHS3 * (2 + -1 * RHS3 * 1) * (-0.3333333333333333 + 1 / (4 * RHS1 * 1));
}

The Symbolics.jl function build_function can also generate code in Julia, MATLAB, and Stan; see the documentation for details and other options.