Implementing CSRK code in BSeries.jl

src/BSeries.jl

@@ -17,7 +17,7 @@ end
 
 using Latexify: Latexify, LaTeXString
 using Combinatorics: Combinatorics, permutations
-using LinearAlgebra: LinearAlgebra, rank
+using LinearAlgebra: LinearAlgebra, rank, dot
 using SparseArrays: SparseArrays, sparse
 
 @reexport using Polynomials: Polynomials, Polynomial
@@ -40,6 +40,8 @@ export renormalize!
 
 export is_energy_preserving, energy_preserving_order
 
+export elementary_differentials_csrk, CSRK
+
 # Types used for traits
 # These traits may decide between different algorithms based on the
 # corresponding complexity etc.
@@ -74,6 +76,20 @@ end
 
 TruncatedBSeries{T, V}() where {T, V} = TruncatedBSeries{T, V}(OrderedDict{T, V}())
 
+"""
+CSRK struct
+"""
+struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
+    matrix::MatT
+end
+
+function CSRK(matrix::AbstractMatrix)
+    T = promote_type(eltype(matrix))
+    _M = T.(matrix)
+    return ContinuousStageRungeKuttaMethod(_M)
+end
+
+
 # general interface methods of `AbstractDict` for `TruncatedBSeries`
 @inline Base.iterate(series::TruncatedBSeries) = iterate(series.coef)
 @inline Base.iterate(series::TruncatedBSeries, state) = iterate(series.coef, state)
@@ -612,6 +628,23 @@ function bseries(ros::RosenbrockMethod, order)
     return series
 end
 
+"""
+bseries CSRK
+"""
+function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
+    csrk = csrk.matrix
+    V = Rational{Int64}
+    series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()
+    series[rootedtree(Int[])] = one(Int64)
+    for o in 1:order
+        for t in RootedTreeIterator(o)
+            series[copy(t)] = elementary_differentials_csrk(csrk, t)
+        end
+    end
+
+    return series
+end
+
 # TODO: bseries(ros::RosenbrockMethod)
 # should create a lazy version, optionally a memoized one
 
@@ -2016,4 +2049,73 @@ function equivalent_trees(tree)
     return equivalent_trees_set
 end
 
+
+
+
+# This function generates a polynomial 
+#       A_{t,z} = [t,t^2/2,..., t^s/s]*M*[1, z, ..., z^(s-1)]^T
+# for a given square matrix M of dimension s and chars 't' and 'z'.  
+function PolynomialA(M,t,z)
+    s = size(M,1)
+    # we need variables to work with
+    variable1 = Sym(t)
+    variable2 = Sym(z)
+    # conjugate the variable 1, provided that this will be the variable
+    # of the left polynomial
+    variable1 = conjugate(variable1)
+    # generate the components of the polynomial with powers of t
+    poli_z = Array{SymPy.Sym}(undef, s)
+    for i in 1:s
+        poli_z[i] = variable2^(i-1)
+    end
+    # generate the components of the polynomial with powers of z
+    poli_t = Array{SymPy.Sym}(undef, s)
+    for i in 1:s
+        poli_t[i] = (1 // i)*(variable1^i)
+    end
+    # multiply matrix times vector
+    result = M * poli_z
+    # use dot product for the two vectors
+    result = dot(poli_t,result)
+    return result
+end
+
+
+"""
+    elementary_differentials_csrk(M,tree)
+
+This function calculates the CSRK elementary differential for a given 
+square matrix 'M' and a given RootedTree.
+
+"""
+function elementary_differentials_csrk(M,rootedtree)
+    # we'll work with the level sequence
+    tree = rootedtree.level_sequence
+    m = maximum(tree)
+    l = length(tree)
+    # create the variables called 'xi' for 1 <= i <= m
+    variables = []
+    for i in 1:m
+        var_name = "x$i"
+        var = Sym(var_name)
+        push!(variables, var)
+    end
+    inverse_counter = l-1
+    # stablish initial integrand, which is the rightmost leaf (last node of the level sequence)
+    if l > 1
+        integrand = integrate(PolynomialA(M,variables[tree[end]-1],variables[tree[end]]),(variables[tree[end]],0,1))
+    else
+        # if the RootedTree is [1] or [], the elementary differential will be 1.
+        return 1
+    end
+    while inverse_counter > 1
+        pseudo_integrand = PolynomialA(M,variables[tree[inverse_counter]-1],variables[tree[inverse_counter]])*integrand
+        integrand = integrate(pseudo_integrand,(variables[tree[inverse_counter]],0,1))
+        inverse_counter -= 1
+    end
+    # multiply for the Basis_Polynomial, i.e. the Polynomial B
+    # return the integral with respect to x1.
+    return integrate(PolynomialA(M,1,variables[1])*integrand,(variables[1],0,1))
+end
+
 end # module