From 44485370c4c1c43262cc69d353dc0a86b0cd500c Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Wed, 12 Nov 2025 02:11:59 +0000 Subject: [PATCH 01/25] Add support for Two-Derivative and Addative Two-Derivative RK Methods Implement TwoDerivativeRungeKuttaMethod and TwoDerivativeAdditiveRungeKuttaMethod structures with associated functions for B-series computations. --- src/BSeries.jl | 599 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 599 insertions(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index 35dd90d5..5ef6ffe8 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1093,6 +1093,333 @@ function RootedTrees.derivative_weight(t::RootedTree, return result end + +##################################################################### +# defining the new methods, not sure where best place for this code is. + +""" + TwoDerivativeRungeKuttaMethod{T, RKm} <: AbstractTimeIntegrationMethod + +Two-derivative Runge–Kutta (TDRK) method composed of multiple +`RungeKuttaMethod`s defining base and derivative components. +""" +struct TwoDerivativeRungeKuttaMethod{T, RKm <: AbstractVector{<:RungeKuttaMethod{T}}} <: RootedTrees.AbstractTimeIntegrationMethod + rkm::RKm +end + +""" + TwoDerivativeRungeKuttaMethod(rkm) + +Construct from a vector of `RungeKuttaMethod`s. +Promotes all coefficient types to a common numeric type. +""" +function TwoDerivativeRungeKuttaMethod(rkm) + T = mapreduce(eltype, promote_type, rkm) # promote coefficient type + As = map(rk -> T.(rk.A), rkm) + bs = map(rk -> T.(rk.b), rkm) + cs = map(rk -> T.(rk.c), rkm) + return TwoDerivativeRungeKuttaMethod(As, bs, cs) +end + +""" + TwoDerivativeRungeKuttaMethod(As, bs[, cs]) + +Construct directly from coefficient arrays. +If c not explicitly given, computed as `sum(A, dims = 2)` per usual. +""" +function TwoDerivativeRungeKuttaMethod(As, bs, cs = map(A -> vec(sum(A, dims = 2)), As)) + rkm = map(RungeKuttaMethod, As, bs, cs) + return TwoDerivativeRungeKuttaMethod(rkm) +end + +""" + eltype(tdrk::TwoDerivativeRungeKuttaMethod) + +Return the element type of the coefficients. +""" +Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T}) where {T} = T + + +""" + TwoDerivativeAdditiveRungeKuttaMethod{T, RKm} <: AbstractTimeIntegrationMethod + +Additive two-derivative Runge–Kutta (TDARK) method composed of multiple +`RungeKuttaMethod`s, typically representing `(f, f′, g)`. +""" +struct TwoDerivativeAdditiveRungeKuttaMethod{T, RKm <: AbstractVector{<:RungeKuttaMethod{T}}} <: RootedTrees.AbstractTimeIntegrationMethod + rkm::RKm +end + +""" + TwoDerivativeAdditiveRungeKuttaMethod(rkm) + +Construct from a vector of `RungeKuttaMethod`s (e.g., `(f, f′, g)`). +Promotes coefficient types to a common numeric type. +""" +function TwoDerivativeAdditiveRungeKuttaMethod(rkm) + T = mapreduce(eltype, promote_type, rkm) # promote element type + As = map(rk -> T.(rk.A), rkm) + bs = map(rk -> T.(rk.b), rkm) + cs = map(rk -> T.(rk.c), rkm) + return TwoDerivativeAdditiveRungeKuttaMethod(As, bs, cs) +end + +""" + TwoDerivativeAdditiveRungeKuttaMethod(As, bs[, cs]) + +Construct directly from coefficient arrays. +""" +function TwoDerivativeAdditiveRungeKuttaMethod(As, bs, cs = map(A -> vec(sum(A, dims = 2)), As)) + rkm = map(RungeKuttaMethod, As, bs, cs) + return TwoDerivativeAdditiveRungeKuttaMethod(rkm) +end + +""" + eltype(tdark::TwoDerivativeAdditiveRungeKuttaMethod) + +Return the coefficient element type. +""" +Base.eltype(tdark::TwoDerivativeAdditiveRungeKuttaMethod{T}) where {T} = T + +##################################################################### + +""" + elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + +Compute the elementary weight associated with the colored rooted tree `t` +for a two-derivative Runge–Kutta method `tdrk`. + +This function first generates all possible collapsed forms of `t` using +[`collapse_tree`](@ref), then sums their contributions weighted by their +multiplicities. Each tree’s contribution is obtained from [`tree_weight`](@ref). + + +This follows from the Butcher type order conditions exhibited in, +Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. +- Numer Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 + + +# Arguments +- `t`: A `ColoredRootedTree` representing the current term. +- `tdrk`: The `TwoDerivativeRungeKuttaMethod` whose coefficients define the + weights. + +# Returns +A scalar weight equal to the sum over all collapsed trees. +""" +function elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + trees, multiplicities = collapse_tree(t) + + sum = zero(eltype(multiplicities)) # accumulator for total weight + for (i, tree) in enumerate(trees) # sum over tree and all collapsed forms + sum += multiplicities[i] * tree_weight(tree, tdrk) + end + return sum +end + + +""" + tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + +Compute the weight contribution of a single colored rooted tree `t` +for the two-derivative Runge–Kutta method `tdrk`. + +The function recursively evaluates the contribution of each subtree +using the appropriate `(A, b, c)` coefficients based on node color: + +- Color `1` → derivative part, F. (uses `a1`, `b1`, `c1`) +- Color `2` → second derivative part, F'. (uses `a2`, `b2`, `c2`) + +# Arguments +- `t`: A `ColoredRootedTree`. +- `tdrk`: A `TwoDerivativeRungeKuttaMethod` containing two embedded RK tableau. + +# Returns +A scalar weight equal to the B-series coefficient associated with the tree. +""" +function tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + b1 = (tdrk.rkm)[1].b # regular b coefficients + b2 = (tdrk.rkm)[2].b # derivative b coefficients + a1 = (tdrk.rkm)[1].A # regular A coefficients + a2 = (tdrk.rkm)[2].A # derivative A coefficients + c1 = (tdrk.rkm)[1].c # regular c coefficients + c2 = (tdrk.rkm)[2].c # derivative c coefficients + + level_sequence = t.level_sequence + color_sequence = t.color_sequence + + # Choose root b vector based on root color + b_root = color_sequence[1] == 1 ? b1 : b2 + + # Recursive evaluation of the product of subtree contributions + function helper(t::ColoredRootedTree) + # Compute contribution of a subtree + function subtree_contribution(t::ColoredRootedTree) + A, c = t.color_sequence[1] == 1 ? (a1, c1) : (a2, c2) + + # Leaf node contributes only its c-value + if length(t.level_sequence) == 1 + return c + else + # Internal node: multiply A by recursive subtree product + return A * helper(t) + end + end + + # Initialize product vector (same length as b1) + product = ones(eltype(b1), length(b1)) + + # Get subtrees of the current node + subtrees_array = subtrees(t) + + # If no subtrees, return ones vector + if isempty(subtrees_array) + return product + else + # Multiply elementwise by each subtree’s contribution + for n in subtrees_array + product = product .* subtree_contribution(n) + end + end + + return product + end + + # Final weight = b_root ⋅ product of subtree contributions + product = helper(t) + result = LinearAlgebra.dot(b_root, product) + return result +end + +""" + elementary_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) -> Number + +Compute the elementary weight associated with the colored rooted tree `t` +for a two-derivative additive Runge–Kutta method `tdark`. + +This function first generates all possible collapsed forms of `t` using +[`collapse_tree`](@ref), then sums their contributions weighted by their +multiplicities. Each tree’s contribution is obtained from +[`tree_weight`](@ref). + +# Arguments +- `t`: A `ColoredRootedTree` representing the current term. +- `tdark`: The `TwoDerivativeAdditiveRungeKuttaMethod` whose coefficients define the weights. + +# Returns +A scalar weight equal to the sum over all collapsed trees. +""" +function elementary_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) + trees, multiplicities = collapse_tree(t) + + sum = 0 + # Sum over the tree and all its collapsed forms, weighted by multiplicity + for (i, tree) in enumerate(trees) + sum += multiplicities[i] * tree_weight(tree, tdark) + end + return sum +end + + +""" + tree_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) -> Number + +Compute the contribution (elementary weight) of a single colored rooted tree `t` +for the two-derivative additive Runge–Kutta method `tdark`. + +Each node color determines which coefficient tableau `(A, b, c)` is used: + +| Color | Meaning | Coefficients used | +|:------|:--------------|:-----------------| +| `0` | Additive `g` | `(a3, b3, c3)` | +| `1` | Regular `f` | `(a1, b1, c1)` | +| `2` | Derivative `f′` | `(a2, b2, c2)` | + +The tree is evaluated recursively by multiplying the contributions of its +subtrees and combining them according to the Butcher recursion. +With the proper A,b or c values according to the level sequence and colors of each according to the color sequence. + +# Arguments +- `t`: A `ColoredRootedTree` representing a single B-series term. +- `tdark`: A `TwoDerivativeAdditiveRungeKuttaMethod` defining all three tableau. + +# Returns +A scalar value representing the coefficient of the tree in the B-series expansion. +""" +function tree_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) + # Key: 0 = additive (g), 1 = regular (f), 2 = derivative (f′) + b1 = (tdark.rkm)[1].b # regular f + b2 = (tdark.rkm)[2].b # derivative f′ + b3 = (tdark.rkm)[3].b # additive g + + a1 = (tdark.rkm)[1].A # regular f + a2 = (tdark.rkm)[2].A # derivative f′ + a3 = (tdark.rkm)[3].A # additive g + + c1 = (tdark.rkm)[1].c # regular f + c2 = (tdark.rkm)[2].c # derivative f′ + c3 = (tdark.rkm)[3].c # additive g + + level_sequence = t.level_sequence + color_sequence = t.color_sequence + + # Choose the root b-vector based on root color + if color_sequence[1] == 0 + b_root = b3 + elseif color_sequence[1] == 1 + b_root = b1 + else + b_root = b2 + end + + # Recursive computation of subtree products + function helper(t::ColoredRootedTree) + # Compute contribution of a subtree recursively + function subtree_contribution(t::ColoredRootedTree) + # Determine coefficient set based on subtree color + if t.color_sequence[1] == 0 + A, c = a3, c3 + elseif t.color_sequence[1] == 1 + A, c = a1, c1 + else + A, c = a2, c2 + end + + # A leaf contributes only its c-vector + if length(t.level_sequence) == 1 + return c + else + # Internal node: multiply A by the recursive subtree product + return A * helper(t) + end + end + + # Initialize with a ones-vector of appropriate length + product = ones(eltype(b1), length(b1)) + + # Get all subtrees of the current node + subtrees_array = subtrees(t) + + # No subtrees → return ones-vector + if isempty(subtrees_array) + return product + else + # Multiply elementwise by each subtree’s contribution + for n in subtrees_array + product = product .* subtree_contribution(n) + end + end + return product + end + + # Compute final result: b_root ⋅ (product of subtree contributions) + product = helper(t) + result = LinearAlgebra.dot(b_root, product) + return result +end + + + # Equation (2.6) of Miyatake & Butcher (2016) function _matrix_a(csrk::ContinuousStageRungeKuttaMethod) M = csrk.matrix @@ -1253,6 +1580,105 @@ end # TODO: bseries(mis::MultirateInfinitesimalSplitMethod) # should create a lazy version, optionally a memoized one + + + +""" + bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) -> TruncatedBSeries + +Construct the truncated B-series of a two-derivative Runge–Kutta method `tdrk` +up to the specified `order`. + +This implementation uses **bicolored rooted trees** to represent the method, +where color 1 corresponds to the the ability for all trees to collapse, which is essential. + +Returns a `TruncatedBSeries{ColoredRootedTree, V}` where `V` is inferred from +the element type of `tdrk`. +""" +function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) + # sanity check: we expect two Runge–Kutta components + if length(tdrk.rkm) != 2 + throw(ArgumentError("Only TwoDerivativeRungeKuttaMethod with a dual splitting is supported. Got $(length(tdrk.rkm)) components.")) + end + + # determine coefficient type + V_tmp = eltype(tdrk) + V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp + + # construct an empty B-series with colored rooted trees + series = TruncatedBSeries{ColoredRootedTree{Int, Vector{Int}, Vector{Int}}, V}() + + # empty (order-0) tree + series[rootedtree(Int[], Int[])] = one(V) + + # iterate over orders + for o in 1:order + for t in BicoloredRootedTreeIterator(o) + # make a colored copy (for now all nodes color=1) since these arnt addative + color_sequence = fill(1, length(t.level_sequence)) + colored_t = rootedtree(t.level_sequence, color_sequence) + + # assign coefficient via elementary weight + series[copy(colored_t)] = elementary_weight(t, tdrk) + end + end + + return series +end + + + +""" + bseries(tdark::TwoDerivativeAdditiveRungeKuttaMethod, order, Flinear, Glinear) + -> TruncatedBSeries + +Construct the truncated B-series for a two-derivative additive Runge–Kutta +method (`TwoDerivativeAdditiveRungeKuttaMethod`) up to the given `order`. +Again this is defined as an addative (think IMEX) method where just one component is paired with its second derivative, could be extended to both having second derivatives in future. + +This version accounts for possible linearity in the functions `F` and `G`, where 'F' is the component with second derivative. + +- `Flinear = true` skips all trees with nonlinear dependence on `F`, + i.e. trees whose root is color `1` with two or more children. +- `Glinear = true` skips all trees with nonlinear dependence on `G`, + i.e. trees whose root is color `0` with two or more children. + +Returns a `TruncatedBSeries{ColoredRootedTree, V}` where `V` is inferred +from the coefficient type of `tdark`. +""" +function bseries(tdark::TwoDerivativeAdditiveRungeKuttaMethod, order, Flinear, Glinear) + if length(tdark.rkm) != 3 + throw(ArgumentError("Expected a TwoDerivativeAdditiveRungeKuttaMethod with three .")) + end + + V_tmp = eltype(tdark) + V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp + + # Initialize empty B-series over colored rooted trees + series = TruncatedBSeries{ColoredRootedTree{Int, Vector{Int}, Vector{Int}}, V}() + + # Add the empty (order-0) tree + series[rootedtree(Int[], Int[])] = one(V) + + # Iterate over all bicolored trees up to the given order + for o in 1:order + for t in BicoloredRootedTreeIterator(o) + # Convert iterator output to a general colored rooted tree + t = rootedtree(t.level_sequence, Int.(t.color_sequence)) + + # Skip trees violating F/G linearity assumptions + if should_skip_tree(t, Flinear, Glinear) + continue + end + + # Assign the elementary weight corresponding to this tree + series[copy(t)] = elementary_weight(t, tdark) + end + end + + return series +end + """ substitute(b, a, t::AbstractRootedTree) @@ -2858,4 +3284,177 @@ function is_symplectic(series::TruncatedBSeries; kwargs...) order_of_symplecticity(series; kwargs...) == order(series) end +##################################################################### +# Multi-Derivative Features + +""" + collapse_tree_at_index(t::ColoredRootedTree, index::Int) -> ColoredRootedTree + +Collapse the node at position `index` in the colored rooted tree `t`. + +This operation is only valid for nodes of color `1`. The parent node of the +collapsed node must also have color `1`; otherwise, the tree is left unchanged. + +The collapse: +- Changes the color of the collapsed node to `-1`. +- Changes the color of its parent node to `2`. +- Decrements the level of all subsequent nodes at deeper levels until the + next node at the same level. +- Removes the collapsed node from the tree. + +Returns the new `ColoredRootedTree` with the collapsed structure. +Throws an error if the node at `index` is not of color `1`. +""" + +function collapse_tree_at_index(t::ColoredRootedTree, index::Int) + if t.color_sequence[index] != 1 + error("Node at index $index is not of type 1 and cannot be collapsed.") + end + level_of_node = t.level_sequence[index] + + #make a copy of the tree to work with + new_tree = deepcopy(t) + + #find the parent by looking for the last node with level one less than the current node + parent = findlast(x -> x == level_of_node - 1, new_tree.level_sequence[1:index-1]) + + #if the parent node cant collapse then we return, that is if the color of the parent is not 1 + if new_tree.color_sequence[parent] != 1 + return + end + + + #main idea + + #change the color of the node to -1 + new_tree.color_sequence[index] = -1 + + #change the color of the parent to 2 + new_tree.color_sequence[parent] = 2 + + #decrement the level of each node after the index untill the next node of the same level + for i in index+1:length(new_tree.level_sequence) + if new_tree.level_sequence[i] > level_of_node + new_tree.level_sequence[i] -= 1 + elseif new_tree.level_sequence[i] == level_of_node + break + end + end + + #take the part of the level and color sequence before and after index and concat them into a new tree + level_sequence = vcat(new_tree.level_sequence[1:index-1], new_tree.level_sequence[index+1:end]) + color_sequence = vcat(new_tree.color_sequence[1:index-1], new_tree.color_sequence[index+1:end]) + new_tree = rootedtree(level_sequence, color_sequence) + + return new_tree +end + + +""" + collapse_tree(t::ColoredRootedTree) -> (trees, multiplicities) + +Generate all possible trees obtained by collapsing any combination of +collapsible nodes in `t`. That is we obtain the set of collapsed variants + +A node is collapsible if its color is `1`. +Each combination of collapses yields a new tree structure; identical trees +are merged and their multiplicities accumulated. + +Returns a tuple `(trees, multiplicities)` where: +- `trees` is a vector of unique collapsed trees, +- `multiplicities` gives how many collapse combinations produce each tree, which is important for order conditions. +""" + +function collapse_tree(t::ColoredRootedTree) + trees = [] + multiplicities = [] + + # Get list of all collapsible node indices (excluding the root at index 1) + collapsable_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices + n = length(collapsable_nodes) + + # Go through all 2^n combinations of collapses + for i in 0:(2^n - 1) + new_tree = deepcopy(t) + num_collapses = 0 + skip_combination = false + + for j in 1:n + if ((i >> (j - 1)) & 1) == 1 + # we lose an index every time we collapse a node so we need to adjust the index if its after ours + adjusted_index = collapsable_nodes[j] - num_collapses + + #check if the adjusted index is still valid + if new_tree.color_sequence[adjusted_index] != 1 + skip_combination = true + break + end + + + new_tree = collapse_tree_at_index(new_tree, adjusted_index) + + if new_tree === nothing + skip_combination = true + break + end + + num_collapses += 1 + end + end + + if skip_combination || new_tree === nothing + continue + end + + # Check if tree is already in the list update tree list or multiplicty respectively + index = findfirst(t -> t == new_tree, trees) + if isnothing(index) + push!(trees, new_tree) + push!(multiplicities, 1) + else + multiplicities[index] += 1 + end + end + + return trees, multiplicities +end + + +""" + should_skip_tree(t::ColoredRootedTree, Flinear::Bool, Glinear::Bool) -> Bool + +Determine whether the colored rooted tree `t` should be skipped when constructing +the B-series, based on the linearity assumptions for `F` and `G`. This allows for development of more methods when one of the functions is linear. + +If `F` (color 1) is linear, any tree whose root node is of color `1` and has +two or more child subtrees (corresponding to higher derivatives of `F`) +is skipped. + +If `G` (color 0) is linear, any tree whose root node is of color `0` and has +two or more child subtrees (corresponding to higher derivatives of `G`) +is skipped. + +Returns `true` if the tree should be ignored, and `false` otherwise. +""" +function should_skip_tree(t::ColoredRootedTree, Flinear::Bool, Glinear::Bool) + subtrees_array = subtrees(t) + + if Flinear + # If F is linear, skip trees where the F-root has more than one child + if t.color_sequence[1] == 1 && length(subtrees_array) > 1 + return true + end + end + + if Glinear + # If G is linear, skip trees where the G-root has more than one child + if t.color_sequence[1] == 0 && length(subtrees_array) > 1 + return true + end + end + + return false +end + + end # module From 0f00454d2ebfc11b49d5f174d3247a1aa99c07c3 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Wed, 12 Nov 2025 15:34:05 +0000 Subject: [PATCH 02/25] Fix typos in BSeries.jl --- src/BSeries.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 5ef6ffe8..d141483c 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1196,7 +1196,7 @@ multiplicities. Each tree’s contribution is obtained from [`tree_weight`](@ref This follows from the Butcher type order conditions exhibited in, Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. -- Numer Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 +- Numer. Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 # Arguments @@ -3332,7 +3332,7 @@ function collapse_tree_at_index(t::ColoredRootedTree, index::Int) #change the color of the parent to 2 new_tree.color_sequence[parent] = 2 - #decrement the level of each node after the index untill the next node of the same level + #decrement the level of each node after the index until the next node of the same level for i in index+1:length(new_tree.level_sequence) if new_tree.level_sequence[i] > level_of_node new_tree.level_sequence[i] -= 1 @@ -3370,8 +3370,8 @@ function collapse_tree(t::ColoredRootedTree) multiplicities = [] # Get list of all collapsible node indices (excluding the root at index 1) - collapsable_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices - n = length(collapsable_nodes) + collapsible_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices + n = length(collapsible_nodes) # Go through all 2^n combinations of collapses for i in 0:(2^n - 1) @@ -3382,7 +3382,7 @@ function collapse_tree(t::ColoredRootedTree) for j in 1:n if ((i >> (j - 1)) & 1) == 1 # we lose an index every time we collapse a node so we need to adjust the index if its after ours - adjusted_index = collapsable_nodes[j] - num_collapses + adjusted_index = collapsible_nodes[j] - num_collapses #check if the adjusted index is still valid if new_tree.color_sequence[adjusted_index] != 1 @@ -3406,7 +3406,7 @@ function collapse_tree(t::ColoredRootedTree) continue end - # Check if tree is already in the list update tree list or multiplicty respectively + # Check if tree is already in the list update tree list or multiplicity respectively index = findfirst(t -> t == new_tree, trees) if isnothing(index) push!(trees, new_tree) From c0d4d78f4a95e7fc4a792adb1a9f30d4ff394081 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Wed, 12 Nov 2025 15:55:36 +0000 Subject: [PATCH 03/25] Clean up BSeries.jl by removing blank lines Removed unnecessary blank lines and comments in BSeries.jl. --- src/BSeries.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index d141483c..8421a380 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1093,7 +1093,6 @@ function RootedTrees.derivative_weight(t::RootedTree, return result end - ##################################################################### # defining the new methods, not sure where best place for this code is. @@ -3456,5 +3455,4 @@ function should_skip_tree(t::ColoredRootedTree, Flinear::Bool, Glinear::Bool) return false end - end # module From e5886fcb0152e79851b7a0f2ec46c134ea78a8a8 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 24 Nov 2025 00:12:30 +0000 Subject: [PATCH 04/25] Implement tests for TwoDerivativeRungeKuttaMethod Added tests for the TwoDerivativeRungeKuttaMethod including constructor validation and accuracy checks. --- test/runtests.jl | 55 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 05acafb8..6c03428a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using BSeries using BSeries.Latexify: latexify -using LinearAlgebra: I +using LinearAlgebra: I, dot using StaticArrays: @SArray, @SMatrix, @SVector using Symbolics: Symbolics @@ -3245,4 +3245,57 @@ using Aqua: Aqua @test substituted == s end end + + @testset "TwoDerivativeRungeKuttaMethod" begin + + # Example two-derivative RK method (order 3) from "On explicit two-derivative Runge-Kutta methods" from R.P.K. Chan and A.Y.J. Tsai (2010) + A1 = [0 0; + 1//2 0] + + b1 = [1, 0] + + A2 = [0 0; + 1//8 0] + + b2 = [1//6, 1//3] + + #constructor + + tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2) + + @test tdrk isa TwoDerivativeRungeKuttaMethod + @test size(tdrk.A1) == (2, 2) + @test size(tdrk.A2) == (2, 2) + @test length(tdrk.b1) == 2 + @test length(tdrk.b2) == 2 + + # row-sum default for c1 + @test tdrk.c1 == vec(sum(A1, dims = 2)) + + # bseries + tdrk_series = @inferred bseries(tdrk, 5) + + #should be 4th order + @test @inferred(order_of_accuracy(tdrk_series)) == 4 + + #test a collapsing tree as well + tree = ColoredRootedTree([1,2,3,4], [1,1,1,1]) + + # Compute elementary weight + w = collapse_elementary_weight(tree, tdrk) + + # Manual expected weight + c = vec(sum(A1, dims=2)) + cp = vec(sum(A2, dims=2)) + b = b1 + bp = b2 + A = A1 + Ap = A2 + + w_expected = dot(b, A * (A * c)) + dot(bp, A * c) + dot(b, Ap * c) + dot(b, A * cp) + dot(bp, cp) + + @test w == w_expected + + end + end # @testset "BSeries" From 9ea2fcbf475b929e390ff74cc74f34feecaa76ae Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 24 Nov 2025 00:15:27 +0000 Subject: [PATCH 05/25] Condensed all multi-derivative functionality into a single section Placed all new functionality in one section. This adds support for multi-derivative RK methods and defines these methods explicitly from two tableau --- src/BSeries.jl | 847 +++++++++++++++++-------------------------------- 1 file changed, 290 insertions(+), 557 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 8421a380..b1fa90f7 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -46,6 +46,8 @@ export is_energy_preserving, energy_preserving_order export order_of_symplecticity, is_symplectic +export TwoDerivativeRungeKuttaMethod, collapse_elementary_weight + # Types used for traits # These traits may decide between different algorithms based on the # corresponding complexity etc. @@ -1093,332 +1095,6 @@ function RootedTrees.derivative_weight(t::RootedTree, return result end -##################################################################### -# defining the new methods, not sure where best place for this code is. - -""" - TwoDerivativeRungeKuttaMethod{T, RKm} <: AbstractTimeIntegrationMethod - -Two-derivative Runge–Kutta (TDRK) method composed of multiple -`RungeKuttaMethod`s defining base and derivative components. -""" -struct TwoDerivativeRungeKuttaMethod{T, RKm <: AbstractVector{<:RungeKuttaMethod{T}}} <: RootedTrees.AbstractTimeIntegrationMethod - rkm::RKm -end - -""" - TwoDerivativeRungeKuttaMethod(rkm) - -Construct from a vector of `RungeKuttaMethod`s. -Promotes all coefficient types to a common numeric type. -""" -function TwoDerivativeRungeKuttaMethod(rkm) - T = mapreduce(eltype, promote_type, rkm) # promote coefficient type - As = map(rk -> T.(rk.A), rkm) - bs = map(rk -> T.(rk.b), rkm) - cs = map(rk -> T.(rk.c), rkm) - return TwoDerivativeRungeKuttaMethod(As, bs, cs) -end - -""" - TwoDerivativeRungeKuttaMethod(As, bs[, cs]) - -Construct directly from coefficient arrays. -If c not explicitly given, computed as `sum(A, dims = 2)` per usual. -""" -function TwoDerivativeRungeKuttaMethod(As, bs, cs = map(A -> vec(sum(A, dims = 2)), As)) - rkm = map(RungeKuttaMethod, As, bs, cs) - return TwoDerivativeRungeKuttaMethod(rkm) -end - -""" - eltype(tdrk::TwoDerivativeRungeKuttaMethod) - -Return the element type of the coefficients. -""" -Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T}) where {T} = T - - -""" - TwoDerivativeAdditiveRungeKuttaMethod{T, RKm} <: AbstractTimeIntegrationMethod - -Additive two-derivative Runge–Kutta (TDARK) method composed of multiple -`RungeKuttaMethod`s, typically representing `(f, f′, g)`. -""" -struct TwoDerivativeAdditiveRungeKuttaMethod{T, RKm <: AbstractVector{<:RungeKuttaMethod{T}}} <: RootedTrees.AbstractTimeIntegrationMethod - rkm::RKm -end - -""" - TwoDerivativeAdditiveRungeKuttaMethod(rkm) - -Construct from a vector of `RungeKuttaMethod`s (e.g., `(f, f′, g)`). -Promotes coefficient types to a common numeric type. -""" -function TwoDerivativeAdditiveRungeKuttaMethod(rkm) - T = mapreduce(eltype, promote_type, rkm) # promote element type - As = map(rk -> T.(rk.A), rkm) - bs = map(rk -> T.(rk.b), rkm) - cs = map(rk -> T.(rk.c), rkm) - return TwoDerivativeAdditiveRungeKuttaMethod(As, bs, cs) -end - -""" - TwoDerivativeAdditiveRungeKuttaMethod(As, bs[, cs]) - -Construct directly from coefficient arrays. -""" -function TwoDerivativeAdditiveRungeKuttaMethod(As, bs, cs = map(A -> vec(sum(A, dims = 2)), As)) - rkm = map(RungeKuttaMethod, As, bs, cs) - return TwoDerivativeAdditiveRungeKuttaMethod(rkm) -end - -""" - eltype(tdark::TwoDerivativeAdditiveRungeKuttaMethod) - -Return the coefficient element type. -""" -Base.eltype(tdark::TwoDerivativeAdditiveRungeKuttaMethod{T}) where {T} = T - -##################################################################### - -""" - elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number - -Compute the elementary weight associated with the colored rooted tree `t` -for a two-derivative Runge–Kutta method `tdrk`. - -This function first generates all possible collapsed forms of `t` using -[`collapse_tree`](@ref), then sums their contributions weighted by their -multiplicities. Each tree’s contribution is obtained from [`tree_weight`](@ref). - - -This follows from the Butcher type order conditions exhibited in, -Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. -- Numer. Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 - - -# Arguments -- `t`: A `ColoredRootedTree` representing the current term. -- `tdrk`: The `TwoDerivativeRungeKuttaMethod` whose coefficients define the - weights. - -# Returns -A scalar weight equal to the sum over all collapsed trees. -""" -function elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) - trees, multiplicities = collapse_tree(t) - - sum = zero(eltype(multiplicities)) # accumulator for total weight - for (i, tree) in enumerate(trees) # sum over tree and all collapsed forms - sum += multiplicities[i] * tree_weight(tree, tdrk) - end - return sum -end - - -""" - tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number - -Compute the weight contribution of a single colored rooted tree `t` -for the two-derivative Runge–Kutta method `tdrk`. - -The function recursively evaluates the contribution of each subtree -using the appropriate `(A, b, c)` coefficients based on node color: - -- Color `1` → derivative part, F. (uses `a1`, `b1`, `c1`) -- Color `2` → second derivative part, F'. (uses `a2`, `b2`, `c2`) - -# Arguments -- `t`: A `ColoredRootedTree`. -- `tdrk`: A `TwoDerivativeRungeKuttaMethod` containing two embedded RK tableau. - -# Returns -A scalar weight equal to the B-series coefficient associated with the tree. -""" -function tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) - b1 = (tdrk.rkm)[1].b # regular b coefficients - b2 = (tdrk.rkm)[2].b # derivative b coefficients - a1 = (tdrk.rkm)[1].A # regular A coefficients - a2 = (tdrk.rkm)[2].A # derivative A coefficients - c1 = (tdrk.rkm)[1].c # regular c coefficients - c2 = (tdrk.rkm)[2].c # derivative c coefficients - - level_sequence = t.level_sequence - color_sequence = t.color_sequence - - # Choose root b vector based on root color - b_root = color_sequence[1] == 1 ? b1 : b2 - - # Recursive evaluation of the product of subtree contributions - function helper(t::ColoredRootedTree) - # Compute contribution of a subtree - function subtree_contribution(t::ColoredRootedTree) - A, c = t.color_sequence[1] == 1 ? (a1, c1) : (a2, c2) - - # Leaf node contributes only its c-value - if length(t.level_sequence) == 1 - return c - else - # Internal node: multiply A by recursive subtree product - return A * helper(t) - end - end - - # Initialize product vector (same length as b1) - product = ones(eltype(b1), length(b1)) - - # Get subtrees of the current node - subtrees_array = subtrees(t) - - # If no subtrees, return ones vector - if isempty(subtrees_array) - return product - else - # Multiply elementwise by each subtree’s contribution - for n in subtrees_array - product = product .* subtree_contribution(n) - end - end - - return product - end - - # Final weight = b_root ⋅ product of subtree contributions - product = helper(t) - result = LinearAlgebra.dot(b_root, product) - return result -end - -""" - elementary_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) -> Number - -Compute the elementary weight associated with the colored rooted tree `t` -for a two-derivative additive Runge–Kutta method `tdark`. - -This function first generates all possible collapsed forms of `t` using -[`collapse_tree`](@ref), then sums their contributions weighted by their -multiplicities. Each tree’s contribution is obtained from -[`tree_weight`](@ref). - -# Arguments -- `t`: A `ColoredRootedTree` representing the current term. -- `tdark`: The `TwoDerivativeAdditiveRungeKuttaMethod` whose coefficients define the weights. - -# Returns -A scalar weight equal to the sum over all collapsed trees. -""" -function elementary_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) - trees, multiplicities = collapse_tree(t) - - sum = 0 - # Sum over the tree and all its collapsed forms, weighted by multiplicity - for (i, tree) in enumerate(trees) - sum += multiplicities[i] * tree_weight(tree, tdark) - end - return sum -end - - -""" - tree_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) -> Number - -Compute the contribution (elementary weight) of a single colored rooted tree `t` -for the two-derivative additive Runge–Kutta method `tdark`. - -Each node color determines which coefficient tableau `(A, b, c)` is used: - -| Color | Meaning | Coefficients used | -|:------|:--------------|:-----------------| -| `0` | Additive `g` | `(a3, b3, c3)` | -| `1` | Regular `f` | `(a1, b1, c1)` | -| `2` | Derivative `f′` | `(a2, b2, c2)` | - -The tree is evaluated recursively by multiplying the contributions of its -subtrees and combining them according to the Butcher recursion. -With the proper A,b or c values according to the level sequence and colors of each according to the color sequence. - -# Arguments -- `t`: A `ColoredRootedTree` representing a single B-series term. -- `tdark`: A `TwoDerivativeAdditiveRungeKuttaMethod` defining all three tableau. - -# Returns -A scalar value representing the coefficient of the tree in the B-series expansion. -""" -function tree_weight(t::ColoredRootedTree, tdark::TwoDerivativeAdditiveRungeKuttaMethod) - # Key: 0 = additive (g), 1 = regular (f), 2 = derivative (f′) - b1 = (tdark.rkm)[1].b # regular f - b2 = (tdark.rkm)[2].b # derivative f′ - b3 = (tdark.rkm)[3].b # additive g - - a1 = (tdark.rkm)[1].A # regular f - a2 = (tdark.rkm)[2].A # derivative f′ - a3 = (tdark.rkm)[3].A # additive g - - c1 = (tdark.rkm)[1].c # regular f - c2 = (tdark.rkm)[2].c # derivative f′ - c3 = (tdark.rkm)[3].c # additive g - - level_sequence = t.level_sequence - color_sequence = t.color_sequence - - # Choose the root b-vector based on root color - if color_sequence[1] == 0 - b_root = b3 - elseif color_sequence[1] == 1 - b_root = b1 - else - b_root = b2 - end - - # Recursive computation of subtree products - function helper(t::ColoredRootedTree) - # Compute contribution of a subtree recursively - function subtree_contribution(t::ColoredRootedTree) - # Determine coefficient set based on subtree color - if t.color_sequence[1] == 0 - A, c = a3, c3 - elseif t.color_sequence[1] == 1 - A, c = a1, c1 - else - A, c = a2, c2 - end - - # A leaf contributes only its c-vector - if length(t.level_sequence) == 1 - return c - else - # Internal node: multiply A by the recursive subtree product - return A * helper(t) - end - end - - # Initialize with a ones-vector of appropriate length - product = ones(eltype(b1), length(b1)) - - # Get all subtrees of the current node - subtrees_array = subtrees(t) - - # No subtrees → return ones-vector - if isempty(subtrees_array) - return product - else - # Multiply elementwise by each subtree’s contribution - for n in subtrees_array - product = product .* subtree_contribution(n) - end - end - return product - end - - # Compute final result: b_root ⋅ (product of subtree contributions) - product = helper(t) - result = LinearAlgebra.dot(b_root, product) - return result -end - - - # Equation (2.6) of Miyatake & Butcher (2016) function _matrix_a(csrk::ContinuousStageRungeKuttaMethod) M = csrk.matrix @@ -1581,6 +1257,56 @@ end +""" + TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c = vec(sum(A1, dims=2))) + +Represent a two-derivative Runge-Kutta method with Butcher coefficients +`A1`, `b1`, and `c` for the first derivative and `A2`, `b2` for the second +derivative. +If `c` is not provided, the usual "row sum" requirement of consistency with +autonomous problems is applied. +``` + +""" +struct TwoDerivativeRungeKuttaMethod{T, + MatT <: AbstractMatrix{T}, + VecT <: AbstractVector{T}} <:RootedTrees.AbstractTimeIntegrationMethod + A1::MatT + b1::VecT + c1::VecT + A2::MatT + b2::VecT +end + +""" + TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2))) + +Construct a two-derivative Runge–Kutta method from coefficient arrays. +All inputs are promoted to a common numeric type. +`c1` defaults to the usual row-sum condition `sum(A1, dims=2)`. +""" +function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2))) + # promote all numeric types together + T = promote_type(eltype(A1), eltype(b1), eltype(A2), eltype(b2), eltype(c1)) + + A1T = T.(A1) + b1T = T.(b1) + c1T = T.(c1) + A2T = T.(A2) + b2T = T.(b2) + + return TwoDerivativeRungeKuttaMethod{T, typeof(A1T), typeof(b1T)}( + A1T, b1T, c1T, A2T, b2T + ) +end + + +""" + eltype(tdrk::TwoDerivativeRungeKuttaMethod) + +Return the element type of the coefficients. +""" +Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T, MatT, VecT}) where {T, MatT, VecT} = T """ bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) -> TruncatedBSeries @@ -1588,94 +1314,273 @@ end Construct the truncated B-series of a two-derivative Runge–Kutta method `tdrk` up to the specified `order`. -This implementation uses **bicolored rooted trees** to represent the method, -where color 1 corresponds to the the ability for all trees to collapse, which is essential. +Returns a `TruncatedBSeries{RootedTree, V}` where `V` is inferred from +the element type of `tdrk`. +""" +function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) + # determine coefficient type + V_tmp = eltype(tdrk) + V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp + + + series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}() + + # empty (order-0) tree + series[rootedtree(Int[])] = one(V) + + # iterate over orders + for o in 1:order + for t in RootedTreeIterator(o) + color_sequence = fill(1, length(t.level_sequence)) + gen_t = rootedtree(t.level_sequence, color_sequence) + # assign coefficient via elementary weight (only def for colored trees) + series[copy(t)] = collapse_elementary_weight(gen_t, tdrk) + end + end + + return series +end + +""" + collapse_elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + +Compute the elementary weight associated with the colored rooted tree `t` +for a two-derivative Runge–Kutta method `tdrk`. + +This function first generates all possible collapsed forms of `t` using +[`collapse_tree`](@ref), then sums their contributions weighted by their +multiplicities. Each tree’s contribution is obtained from [`tree_weight`](@ref). + + +This follows from the Butcher type order conditions exhibited in, +Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. +- Numer. Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 + + +# Arguments +- `t`: A `ColoredRootedTree` representing the current term. +- `tdrk`: The `TwoDerivativeRungeKuttaMethod` whose coefficients define the + weights. + +# Returns +A scalar weight equal to the sum over all collapsed trees. +""" +function collapse_elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + trees, multiplicities = collapse_tree(t) + + sum = 0 + for (i, tree) in enumerate(trees) # sum over tree and all collapsed forms + sum += multiplicities[i] * tree_weight(tree, tdrk) + end + return sum +end + + +""" + tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + +Compute the weight contribution of a single colored rooted tree `t` +for the two-derivative Runge–Kutta method `tdrk`. + +The function recursively evaluates the contribution of each subtree +using the appropriate `(A, b, c)` coefficients based on node color: + +- Color `1` → derivative part, F. (uses `a1`, `b1`, `c1`) +- Color `2` → second derivative part, F'. (uses `a2`, `b2`, `c2`) + +# Arguments +- `t`: A `ColoredRootedTree`. +- `tdrk`: A `TwoDerivativeRungeKuttaMethod` containing two embedded RK tableau. + +# Returns +A scalar weight equal to the B-series coefficient associated with the tree. +""" +function tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + b1 = tdrk.b1 + b2 = tdrk.b2 + a1 = tdrk.A1 + a2 = tdrk.A2 + c1 = tdrk.c1 + c2 = sum(a2, dims=2) + + level_sequence = t.level_sequence + color_sequence = t.color_sequence + + # Choose root b vector based on root color + b_root = color_sequence[1] == 1 ? b1 : b2 + + # Recursive evaluation of the product of subtree contributions + function helper(t::ColoredRootedTree) + # Compute contribution of a subtree + function subtree_contribution(t::ColoredRootedTree) + A, c = t.color_sequence[1] == 1 ? (a1, c1) : (a2, c2) + + # Leaf node contributes only its c-value + if length(t.level_sequence) == 1 + return c + else + # Internal node: multiply A by recursive subtree product + return A * helper(t) + end + end + + # Initialize product vector (same length as b1) + product = ones(eltype(b1), length(b1)) + + # Get subtrees of the current node + subtrees_array = subtrees(t) + + # If no subtrees, return ones vector + if isempty(subtrees_array) + return product + else + # Multiply elementwise by each subtree’s contribution + for n in subtrees_array + product = product .* subtree_contribution(n) + end + end + + return product + end + + # Final weight = b_root ⋅ product of subtree contributions + product = helper(t) + result = LinearAlgebra.dot(b_root, product) + return result +end + +# Multi-Derivative Features + +""" + collapse_tree_at_index(t::ColoredRootedTree, index::Int) -> ColoredRootedTree + +Collapse the node at position `index` in the colored rooted tree `t`. + +This operation is only valid for nodes of color `1`. The parent node of the +collapsed node must also have color `1`; otherwise, the tree is left unchanged. + +The collapse: +- Changes the color of the collapsed node to `-1`. +- Changes the color of its parent node to `2`. +- Decrements the level of all subsequent nodes at deeper levels until the + next node at the same level. +- Removes the collapsed node from the tree. -Returns a `TruncatedBSeries{ColoredRootedTree, V}` where `V` is inferred from -the element type of `tdrk`. +Returns the new `ColoredRootedTree` with the collapsed structure. +Throws an error if the node at `index` is not of color `1`. """ -function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) - # sanity check: we expect two Runge–Kutta components - if length(tdrk.rkm) != 2 - throw(ArgumentError("Only TwoDerivativeRungeKuttaMethod with a dual splitting is supported. Got $(length(tdrk.rkm)) components.")) + +function collapse_tree_at_index(t::ColoredRootedTree, index::Int) + if t.color_sequence[index] != 1 + error("Node at index $index is not of type 1 and cannot be collapsed.") end + level_of_node = t.level_sequence[index] - # determine coefficient type - V_tmp = eltype(tdrk) - V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp + #make a copy of the tree to work with + new_tree = deepcopy(t) + + #find the parent by looking for the last node with level one less than the current node + parent = findlast(x -> x == level_of_node - 1, new_tree.level_sequence[1:index-1]) - # construct an empty B-series with colored rooted trees - series = TruncatedBSeries{ColoredRootedTree{Int, Vector{Int}, Vector{Int}}, V}() + #if the parent node cant collapse then we return, that is if the color of the parent is not 1 + if new_tree.color_sequence[parent] != 1 + return + end + + + #main idea - # empty (order-0) tree - series[rootedtree(Int[], Int[])] = one(V) + #change the color of the node to -1 + new_tree.color_sequence[index] = -1 - # iterate over orders - for o in 1:order - for t in BicoloredRootedTreeIterator(o) - # make a colored copy (for now all nodes color=1) since these arnt addative - color_sequence = fill(1, length(t.level_sequence)) - colored_t = rootedtree(t.level_sequence, color_sequence) + #change the color of the parent to 2 + new_tree.color_sequence[parent] = 2 - # assign coefficient via elementary weight - series[copy(colored_t)] = elementary_weight(t, tdrk) + #decrement the level of each node after the index until the next node of the same level + for i in index+1:length(new_tree.level_sequence) + if new_tree.level_sequence[i] > level_of_node + new_tree.level_sequence[i] -= 1 + elseif new_tree.level_sequence[i] == level_of_node + break end end - return series + #take the part of the level and color sequence before and after index and concat them into a new tree + level_sequence = vcat(new_tree.level_sequence[1:index-1], new_tree.level_sequence[index+1:end]) + color_sequence = vcat(new_tree.color_sequence[1:index-1], new_tree.color_sequence[index+1:end]) + new_tree = rootedtree(level_sequence, color_sequence) + + return new_tree end +""" + collapse_tree(t::ColoredRootedTree) -> (trees, multiplicities) + +Generate all possible trees obtained by collapsing any combination of +collapsible nodes in `t`. That is we obtain the set of collapsed variants + +A node is collapsible if its color is `1`. +Each combination of collapses yields a new tree structure; identical trees +are merged and their multiplicities accumulated. +Returns a tuple `(trees, multiplicities)` where: +- `trees` is a vector of unique collapsed trees, +- `multiplicities` gives how many collapse combinations produce each tree, which is important for order conditions. """ - bseries(tdark::TwoDerivativeAdditiveRungeKuttaMethod, order, Flinear, Glinear) - -> TruncatedBSeries -Construct the truncated B-series for a two-derivative additive Runge–Kutta -method (`TwoDerivativeAdditiveRungeKuttaMethod`) up to the given `order`. -Again this is defined as an addative (think IMEX) method where just one component is paired with its second derivative, could be extended to both having second derivatives in future. +function collapse_tree(t::ColoredRootedTree) + trees = [] + multiplicities = [] + + # Get list of all collapsible node indices (excluding the root at index 1) + collapsible_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices + n = length(collapsible_nodes) -This version accounts for possible linearity in the functions `F` and `G`, where 'F' is the component with second derivative. + # Go through all 2^n combinations of collapses + for i in 0:(2^n - 1) + new_tree = deepcopy(t) + num_collapses = 0 + skip_combination = false -- `Flinear = true` skips all trees with nonlinear dependence on `F`, - i.e. trees whose root is color `1` with two or more children. -- `Glinear = true` skips all trees with nonlinear dependence on `G`, - i.e. trees whose root is color `0` with two or more children. + for j in 1:n + if ((i >> (j - 1)) & 1) == 1 + # we lose an index every time we collapse a node so we need to adjust the index if its after ours + adjusted_index = collapsible_nodes[j] - num_collapses -Returns a `TruncatedBSeries{ColoredRootedTree, V}` where `V` is inferred -from the coefficient type of `tdark`. -""" -function bseries(tdark::TwoDerivativeAdditiveRungeKuttaMethod, order, Flinear, Glinear) - if length(tdark.rkm) != 3 - throw(ArgumentError("Expected a TwoDerivativeAdditiveRungeKuttaMethod with three .")) - end + #check if the adjusted index is still valid + if new_tree.color_sequence[adjusted_index] != 1 + skip_combination = true + break + end - V_tmp = eltype(tdark) - V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp - # Initialize empty B-series over colored rooted trees - series = TruncatedBSeries{ColoredRootedTree{Int, Vector{Int}, Vector{Int}}, V}() + new_tree = collapse_tree_at_index(new_tree, adjusted_index) - # Add the empty (order-0) tree - series[rootedtree(Int[], Int[])] = one(V) + if new_tree === nothing + skip_combination = true + break + end - # Iterate over all bicolored trees up to the given order - for o in 1:order - for t in BicoloredRootedTreeIterator(o) - # Convert iterator output to a general colored rooted tree - t = rootedtree(t.level_sequence, Int.(t.color_sequence)) - - # Skip trees violating F/G linearity assumptions - if should_skip_tree(t, Flinear, Glinear) - continue + num_collapses += 1 end + end + + if skip_combination || new_tree === nothing + continue + end - # Assign the elementary weight corresponding to this tree - series[copy(t)] = elementary_weight(t, tdark) + # Check if tree is already in the list update tree list or multiplicity respectively + index = findfirst(t -> t == new_tree, trees) + if isnothing(index) + push!(trees, new_tree) + push!(multiplicities, 1) + else + multiplicities[index] += 1 end end - return series + return trees, multiplicities end """ @@ -3283,176 +3188,4 @@ function is_symplectic(series::TruncatedBSeries; kwargs...) order_of_symplecticity(series; kwargs...) == order(series) end -##################################################################### -# Multi-Derivative Features - -""" - collapse_tree_at_index(t::ColoredRootedTree, index::Int) -> ColoredRootedTree - -Collapse the node at position `index` in the colored rooted tree `t`. - -This operation is only valid for nodes of color `1`. The parent node of the -collapsed node must also have color `1`; otherwise, the tree is left unchanged. - -The collapse: -- Changes the color of the collapsed node to `-1`. -- Changes the color of its parent node to `2`. -- Decrements the level of all subsequent nodes at deeper levels until the - next node at the same level. -- Removes the collapsed node from the tree. - -Returns the new `ColoredRootedTree` with the collapsed structure. -Throws an error if the node at `index` is not of color `1`. -""" - -function collapse_tree_at_index(t::ColoredRootedTree, index::Int) - if t.color_sequence[index] != 1 - error("Node at index $index is not of type 1 and cannot be collapsed.") - end - level_of_node = t.level_sequence[index] - - #make a copy of the tree to work with - new_tree = deepcopy(t) - - #find the parent by looking for the last node with level one less than the current node - parent = findlast(x -> x == level_of_node - 1, new_tree.level_sequence[1:index-1]) - - #if the parent node cant collapse then we return, that is if the color of the parent is not 1 - if new_tree.color_sequence[parent] != 1 - return - end - - - #main idea - - #change the color of the node to -1 - new_tree.color_sequence[index] = -1 - - #change the color of the parent to 2 - new_tree.color_sequence[parent] = 2 - - #decrement the level of each node after the index until the next node of the same level - for i in index+1:length(new_tree.level_sequence) - if new_tree.level_sequence[i] > level_of_node - new_tree.level_sequence[i] -= 1 - elseif new_tree.level_sequence[i] == level_of_node - break - end - end - - #take the part of the level and color sequence before and after index and concat them into a new tree - level_sequence = vcat(new_tree.level_sequence[1:index-1], new_tree.level_sequence[index+1:end]) - color_sequence = vcat(new_tree.color_sequence[1:index-1], new_tree.color_sequence[index+1:end]) - new_tree = rootedtree(level_sequence, color_sequence) - - return new_tree -end - - -""" - collapse_tree(t::ColoredRootedTree) -> (trees, multiplicities) - -Generate all possible trees obtained by collapsing any combination of -collapsible nodes in `t`. That is we obtain the set of collapsed variants - -A node is collapsible if its color is `1`. -Each combination of collapses yields a new tree structure; identical trees -are merged and their multiplicities accumulated. - -Returns a tuple `(trees, multiplicities)` where: -- `trees` is a vector of unique collapsed trees, -- `multiplicities` gives how many collapse combinations produce each tree, which is important for order conditions. -""" - -function collapse_tree(t::ColoredRootedTree) - trees = [] - multiplicities = [] - - # Get list of all collapsible node indices (excluding the root at index 1) - collapsible_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices - n = length(collapsible_nodes) - - # Go through all 2^n combinations of collapses - for i in 0:(2^n - 1) - new_tree = deepcopy(t) - num_collapses = 0 - skip_combination = false - - for j in 1:n - if ((i >> (j - 1)) & 1) == 1 - # we lose an index every time we collapse a node so we need to adjust the index if its after ours - adjusted_index = collapsible_nodes[j] - num_collapses - - #check if the adjusted index is still valid - if new_tree.color_sequence[adjusted_index] != 1 - skip_combination = true - break - end - - - new_tree = collapse_tree_at_index(new_tree, adjusted_index) - - if new_tree === nothing - skip_combination = true - break - end - - num_collapses += 1 - end - end - - if skip_combination || new_tree === nothing - continue - end - - # Check if tree is already in the list update tree list or multiplicity respectively - index = findfirst(t -> t == new_tree, trees) - if isnothing(index) - push!(trees, new_tree) - push!(multiplicities, 1) - else - multiplicities[index] += 1 - end - end - - return trees, multiplicities -end - - -""" - should_skip_tree(t::ColoredRootedTree, Flinear::Bool, Glinear::Bool) -> Bool - -Determine whether the colored rooted tree `t` should be skipped when constructing -the B-series, based on the linearity assumptions for `F` and `G`. This allows for development of more methods when one of the functions is linear. - -If `F` (color 1) is linear, any tree whose root node is of color `1` and has -two or more child subtrees (corresponding to higher derivatives of `F`) -is skipped. - -If `G` (color 0) is linear, any tree whose root node is of color `0` and has -two or more child subtrees (corresponding to higher derivatives of `G`) -is skipped. - -Returns `true` if the tree should be ignored, and `false` otherwise. -""" -function should_skip_tree(t::ColoredRootedTree, Flinear::Bool, Glinear::Bool) - subtrees_array = subtrees(t) - - if Flinear - # If F is linear, skip trees where the F-root has more than one child - if t.color_sequence[1] == 1 && length(subtrees_array) > 1 - return true - end - end - - if Glinear - # If G is linear, skip trees where the G-root has more than one child - if t.color_sequence[1] == 0 && length(subtrees_array) > 1 - return true - end - end - - return false -end - end # module From 5e72636d539663c3b1b0179618743c9cd88352f3 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 24 Nov 2025 15:13:26 +0000 Subject: [PATCH 06/25] Update src/BSeries.jl add comments Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index b1fa90f7..a278cf1e 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1320,9 +1320,11 @@ the element type of `tdrk`. function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) # determine coefficient type V_tmp = eltype(tdrk) + # If people use integer coefficients, they will likely want to have results + # as exact as possible. However, general terms are not integers. Thus, we + # use rationals instead. V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp - series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}() # empty (order-0) tree From fe35120b6911f7491ed74ea35961f72e9a1495dc Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 24 Nov 2025 15:15:13 +0000 Subject: [PATCH 07/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index a278cf1e..306d0161 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1278,13 +1278,6 @@ struct TwoDerivativeRungeKuttaMethod{T, b2::VecT end -""" - TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2))) - -Construct a two-derivative Runge–Kutta method from coefficient arrays. -All inputs are promoted to a common numeric type. -`c1` defaults to the usual row-sum condition `sum(A1, dims=2)`. -""" function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2))) # promote all numeric types together T = promote_type(eltype(A1), eltype(b1), eltype(A2), eltype(b2), eltype(c1)) From bb3dcb03d55f347716f169dca37e22eb08bd8fef Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Tue, 9 Dec 2025 22:13:55 +0000 Subject: [PATCH 08/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 306d0161..b67ef931 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1294,12 +1294,7 @@ function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2)) end -""" - eltype(tdrk::TwoDerivativeRungeKuttaMethod) - -Return the element type of the coefficients. -""" -Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T, MatT, VecT}) where {T, MatT, VecT} = T +Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T}) where {T} = T """ bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) -> TruncatedBSeries From fef7a7d753d0b97b7d4adc5a4885c3b148652cf1 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Tue, 9 Dec 2025 22:14:04 +0000 Subject: [PATCH 09/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index b67ef931..7caa5dad 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1265,6 +1265,15 @@ Represent a two-derivative Runge-Kutta method with Butcher coefficients derivative. If `c` is not provided, the usual "row sum" requirement of consistency with autonomous problems is applied. + +Given an ODE ``u'(t) = f(t, u(t))`` with ``u''(t) = g(t, u(t))``, +one step from ``u^{n}`` to ``u^{n+1}`` is given by + +```math +\begin{aligned} + y^i &= u^n + \Delta t \sum_j a^{1}_{i,j} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_j a^{2}_{i,j} g(t^n + c_i \Delta t, y^i), \\ + u^{n+1} &= u^n + \Delta t \sum_i b^1_{i} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_i b^2_{i} g(t^n + c_i \Delta t, y^i), +\end{aligned} ``` """ From dd516e24fc30f785972c001cb41e29dcd8666b71 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Wed, 10 Dec 2025 01:10:23 +0000 Subject: [PATCH 10/25] Change to Recursive method of getting OC's Refactor order condition calculation to use the more straightforward recursive method introduced in the paper by Chan and Tsai (2010) following closely the implementation by Lisann Radmacher in her undergraduate thesis. This approach improves complexity compared to computing and storing all subtrees before calculating weights. --- src/BSeries.jl | 351 ++++++++++++++++--------------------------------- 1 file changed, 111 insertions(+), 240 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 7caa5dad..3b9f784c 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -19,7 +19,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 @@ -46,7 +46,7 @@ export is_energy_preserving, energy_preserving_order export order_of_symplecticity, is_symplectic -export TwoDerivativeRungeKuttaMethod, collapse_elementary_weight +export TwoDerivativeRungeKuttaMethod # Types used for traits # These traits may decide between different algorithms based on the @@ -1255,8 +1255,6 @@ end # TODO: bseries(mis::MultirateInfinitesimalSplitMethod) # should create a lazy version, optionally a memoized one - - """ TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c = vec(sum(A1, dims=2))) @@ -1269,12 +1267,12 @@ autonomous problems is applied. Given an ODE ``u'(t) = f(t, u(t))`` with ``u''(t) = g(t, u(t))``, one step from ``u^{n}`` to ``u^{n+1}`` is given by -```math -\begin{aligned} - y^i &= u^n + \Delta t \sum_j a^{1}_{i,j} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_j a^{2}_{i,j} g(t^n + c_i \Delta t, y^i), \\ - u^{n+1} &= u^n + \Delta t \sum_i b^1_{i} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_i b^2_{i} g(t^n + c_i \Delta t, y^i), -\end{aligned} -``` +'''math +\\begin{aligned} + y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_i \\Delta t, y^i), \\\ + u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i), +\\end{aligned} +''' """ struct TwoDerivativeRungeKuttaMethod{T, @@ -1285,9 +1283,10 @@ struct TwoDerivativeRungeKuttaMethod{T, c1::VecT A2::MatT b2::VecT + c2::VecT end -function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2))) +function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2)), c2 = vec(sum(A2, dims=2))) # promote all numeric types together T = promote_type(eltype(A1), eltype(b1), eltype(A2), eltype(b2), eltype(c1)) @@ -1296,9 +1295,10 @@ function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2)) c1T = T.(c1) A2T = T.(A2) b2T = T.(b2) + c2T = T.(c2) return TwoDerivativeRungeKuttaMethod{T, typeof(A1T), typeof(b1T)}( - A1T, b1T, c1T, A2T, b2T + A1T, b1T, c1T, A2T, b2T, c2T ) end @@ -1317,23 +1317,20 @@ the element type of `tdrk`. function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) # determine coefficient type V_tmp = eltype(tdrk) - # If people use integer coefficients, they will likely want to have results - # as exact as possible. However, general terms are not integers. Thus, we - # use rationals instead. - V = V_tmp <: Integer ? Rational{V_tmp} : V_tmp - + if V_tmp <: Integer + # If people use integer coefficients, they will likely want to have results + # as exact as possible. However, general terms are not integers. Thus, we + # use rationals instead. + V = Rational{V_tmp} + else + V = V_tmp + end series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}() - # empty (order-0) tree series[rootedtree(Int[])] = one(V) - - # iterate over orders for o in 1:order for t in RootedTreeIterator(o) - color_sequence = fill(1, length(t.level_sequence)) - gen_t = rootedtree(t.level_sequence, color_sequence) - # assign coefficient via elementary weight (only def for colored trees) - series[copy(t)] = collapse_elementary_weight(gen_t, tdrk) + series[copy(t)] = elementary_weight(t, tdrk) end end @@ -1341,275 +1338,149 @@ function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) end """ - collapse_elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number -Compute the elementary weight associated with the colored rooted tree `t` +Compute the elementary weight associated with the rooted tree `t` for a two-derivative Runge–Kutta method `tdrk`. -This function first generates all possible collapsed forms of `t` using -[`collapse_tree`](@ref), then sums their contributions weighted by their -multiplicities. Each tree’s contribution is obtained from [`tree_weight`](@ref). - -This follows from the Butcher type order conditions exhibited in, +This follows the recursive formula for the Butcher type order conditions exhibited in, Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. - Numer. Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 - + +#see formula 16 in the paper +# alpha(t) = b1*eta(subtrees(t)) +b2*eta(collapse_trees(nu)) # Arguments -- `t`: A `ColoredRootedTree` representing the current term. +- `t`: A `RootedTree` representing the current term. - `tdrk`: The `TwoDerivativeRungeKuttaMethod` whose coefficients define the weights. # Returns A scalar weight equal to the sum over all collapsed trees. """ -function collapse_elementary_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) - trees, multiplicities = collapse_tree(t) - - sum = 0 - for (i, tree) in enumerate(trees) # sum over tree and all collapsed forms - sum += multiplicities[i] * tree_weight(tree, tdrk) - end - return sum +function elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + b1 = tdrk.b1 + b2 = tdrk.b2 + # alpha(t) = b1*eta(subtrees(t)) +b2*eta(collapse_trees(nu)) + dot(b1, derivative_weight(t, tdrk)) + dot(b2, collapsed_derivative_weight(t, tdrk)) end - """ - tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -Compute the weight contribution of a single colored rooted tree `t` -for the two-derivative Runge–Kutta method `tdrk`. +Compute the derivative weight for the standered trees `t` in a two-derivative Runge–Kutta (TDRK) method. -The function recursively evaluates the contribution of each subtree -using the appropriate `(A, b, c)` coefficients based on node color: +this corresponds to formula 15 in Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. -- Color `1` → derivative part, F. (uses `a1`, `b1`, `c1`) -- Color `2` → second derivative part, F'. (uses `a2`, `b2`, `c2`) +eta(t) = A1*eta(t) + A2*eta(t/[1,2]) -# Arguments -- `t`: A `ColoredRootedTree`. -- `tdrk`: A `TwoDerivativeRungeKuttaMethod` containing two embedded RK tableau. +where we are evaluating eta(t) for the elementary weight of the tree t -# Returns -A scalar weight equal to the B-series coefficient associated with the tree. """ -function tree_weight(t::ColoredRootedTree, tdrk::TwoDerivativeRungeKuttaMethod) - b1 = tdrk.b1 - b2 = tdrk.b2 - a1 = tdrk.A1 - a2 = tdrk.A2 +function derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + A1 = tdrk.A1 c1 = tdrk.c1 - c2 = sum(a2, dims=2) - - level_sequence = t.level_sequence - color_sequence = t.color_sequence - - # Choose root b vector based on root color - b_root = color_sequence[1] == 1 ? b1 : b2 - - # Recursive evaluation of the product of subtree contributions - function helper(t::ColoredRootedTree) - # Compute contribution of a subtree - function subtree_contribution(t::ColoredRootedTree) - A, c = t.color_sequence[1] == 1 ? (a1, c1) : (a2, c2) - - # Leaf node contributes only its c-value - if length(t.level_sequence) == 1 - return c - else - # Internal node: multiply A by recursive subtree product - return A * helper(t) - end - end + A2 = tdrk.A2 + c2 = tdrk.c2 - # Initialize product vector (same length as b1) - product = ones(eltype(b1), length(b1)) + result1 = zero(c1) .+ one(eltype(c1)) - # Get subtrees of the current node - subtrees_array = subtrees(t) - - # If no subtrees, return ones vector - if isempty(subtrees_array) - return product - else - # Multiply elementwise by each subtree’s contribution - for n in subtrees_array - product = product .* subtree_contribution(n) - end + if t == rootedtree(Int64[]) || t == rootedtree([1]) + return zero(c1) .+ one(eltype(c1)) + else + subtrees_arr = subtrees(t) + l = 1 + for n in SubtreeIterator(t) + tmp = A1 * derivative_weight(subtrees_arr[l], tdrk) .+ + A2 * collapsed_derivative_weight(subtrees_arr[l], tdrk) + result1 = result1 .* tmp + l += 1 end - - return product + return result1 end - - # Final weight = b_root ⋅ product of subtree contributions - product = helper(t) - result = LinearAlgebra.dot(b_root, product) - return result end - -# Multi-Derivative Features - """ - collapse_tree_at_index(t::ColoredRootedTree, index::Int) -> ColoredRootedTree - -Collapse the node at position `index` in the colored rooted tree `t`. - -This operation is only valid for nodes of color `1`. The parent node of the -collapsed node must also have color `1`; otherwise, the tree is left unchanged. + collapsed_derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -The collapse: -- Changes the color of the collapsed node to `-1`. -- Changes the color of its parent node to `2`. -- Decrements the level of all subsequent nodes at deeper levels until the - next node at the same level. -- Removes the collapsed node from the tree. -Returns the new `ColoredRootedTree` with the collapsed structure. -Throws an error if the node at `index` is not of color `1`. -""" - -function collapse_tree_at_index(t::ColoredRootedTree, index::Int) - if t.color_sequence[index] != 1 - error("Node at index $index is not of type 1 and cannot be collapsed.") - end - level_of_node = t.level_sequence[index] +Compute the derivative weight for the standered trees `t` in a two-derivative Runge–Kutta (TDRK) method. - #make a copy of the tree to work with - new_tree = deepcopy(t) +this corresponds to formula 15 in Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. - #find the parent by looking for the last node with level one less than the current node - parent = findlast(x -> x == level_of_node - 1, new_tree.level_sequence[1:index-1]) +eta(t) = A1*eta(t) + A2*eta(t\\[1,2]) - #if the parent node cant collapse then we return, that is if the color of the parent is not 1 - if new_tree.color_sequence[parent] != 1 - return - end +where we are evaluating eta(t\\[1/2]) part for the elementary weight of the tree t +""" +function collapsed_derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) + A1 = tdrk.A1 + c1 = tdrk.c1 + A2 = tdrk.A2 + c2 = tdrk.c2 + result = zero(c2) - #main idea - - #change the color of the node to -1 - new_tree.color_sequence[index] = -1 - - #change the color of the parent to 2 - new_tree.color_sequence[parent] = 2 - - #decrement the level of each node after the index until the next node of the same level - for i in index+1:length(new_tree.level_sequence) - if new_tree.level_sequence[i] > level_of_node - new_tree.level_sequence[i] -= 1 - elseif new_tree.level_sequence[i] == level_of_node - break - end - end - - #take the part of the level and color sequence before and after index and concat them into a new tree - level_sequence = vcat(new_tree.level_sequence[1:index-1], new_tree.level_sequence[index+1:end]) - color_sequence = vcat(new_tree.color_sequence[1:index-1], new_tree.color_sequence[index+1:end]) - new_tree = rootedtree(level_sequence, color_sequence) - - return new_tree -end - - -""" - collapse_tree(t::ColoredRootedTree) -> (trees, multiplicities) - -Generate all possible trees obtained by collapsing any combination of -collapsible nodes in `t`. That is we obtain the set of collapsed variants - -A node is collapsible if its color is `1`. -Each combination of collapses yields a new tree structure; identical trees -are merged and their multiplicities accumulated. - -Returns a tuple `(trees, multiplicities)` where: -- `trees` is a vector of unique collapsed trees, -- `multiplicities` gives how many collapse combinations produce each tree, which is important for order conditions. -""" - -function collapse_tree(t::ColoredRootedTree) - trees = [] - multiplicities = [] - - # Get list of all collapsible node indices (excluding the root at index 1) - collapsible_nodes = findall(t.color_sequence[2:end] .== 1) .+ 1 #plus one is to adjust since [2:end] shifts indices - n = length(collapsible_nodes) - - # Go through all 2^n combinations of collapses - for i in 0:(2^n - 1) - new_tree = deepcopy(t) - num_collapses = 0 - skip_combination = false - - for j in 1:n - if ((i >> (j - 1)) & 1) == 1 - # we lose an index every time we collapse a node so we need to adjust the index if its after ours - adjusted_index = collapsible_nodes[j] - num_collapses - - #check if the adjusted index is still valid - if new_tree.color_sequence[adjusted_index] != 1 - skip_combination = true - break - end - - - new_tree = collapse_tree_at_index(new_tree, adjusted_index) - - if new_tree === nothing - skip_combination = true - break - end - - num_collapses += 1 + if t == rootedtree(Int64[]) + return zero(c1) .+ one(eltype(c1)) + else + collapsed_trees = collapse_tree(t) + number_of_trees = length(collapsed_trees) + + for k in 1:number_of_trees + treecombinations = collapsed_trees[k] + number2 = length(treecombinations) + sum = zero(c1) .+ one(eltype(c1)) + + for m in 1:number2 + step = A1 * derivative_weight(treecombinations[m], tdrk) .+ + A2 * collapsed_derivative_weight(treecombinations[m], tdrk) + sum = sum .* step end - end - if skip_combination || new_tree === nothing - continue + result = result .+ sum end - # Check if tree is already in the list update tree list or multiplicity respectively - index = findfirst(t -> t == new_tree, trees) - if isnothing(index) - push!(trees, new_tree) - push!(multiplicities, 1) - else - multiplicities[index] += 1 - end + return result end - - return trees, multiplicities end +# Multi-Derivative Features + """ - substitute(b, a, t::AbstractRootedTree) + collapse_tree(t::RootedTree) -Compute the coefficient corresponding to the tree `t` of the B-series that is -formed by substituting the B-series `b` into the B-series `a`. It is assumed -that the B-series `b` has the coefficient zero of the empty tree. +recursively collapse a rooted tree `t` by removing [1,2] type branches. -# References +A collapse groups the children of `t` into all possible merged subsets, +corresponding to the combinatorial partitions needed for the collapsed +derivative weights 'eta(t\\[1,2])` in two-derivative B-series. -Section 3.2 of -- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010) - Algebraic Structures of B-series. - Foundations of Computational Mathematics - [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1) +Each element of the returned array is one valid list of collapsed +subtrees of `t`. No modification of `t` is performed. """ -function substitute(b, a, t::AbstractRootedTree) - result = zero(first(values(a)) * first(values(b))) +function collapse_tree(t::RootedTree) + CollapsedArray = [] - for (forest, skeleton) in PartitionIterator(t) - update = a[skeleton] - update isa Rational && iszero(update) && continue - for tree in forest - update *= b[tree] + subtrees_arr = subtrees(t) + subtrees_multiplicity = length(subtrees_arr) + + # Recustive approach to create all the possibilities of subtrees + for i in 1:subtrees_multiplicity + subsubtrees = subtrees(subtrees_arr[i]) + numberofsubsubtrees = length(subsubtrees) + + for j in 1:subtrees_multiplicity + if j == i + elseif j > i + push!(subsubtrees, subtrees_arr[j]) + else j < i + pushfirst!(subsubtrees, subtrees_arr[i-j]) + end end - result += update + push!(CollapsedArray, subsubtrees) end - return result + return CollapsedArray end """ From 40e8f9ee9ebd7af2fa71bbeb0cfc79dff289ccc9 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Wed, 10 Dec 2025 01:12:03 +0000 Subject: [PATCH 11/25] Update tests for TwoDerivativeRungeKuttaMethod Add a test for multi-derivative Runge-Kutta methods, specifically a 5s7p method from Chan and Tsai (2010). --- test/runtests.jl | 50 +++++++++++++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6c03428a..59282992 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3278,24 +3278,38 @@ using Aqua: Aqua #should be 4th order @test @inferred(order_of_accuracy(tdrk_series)) == 4 - #test a collapsing tree as well - tree = ColoredRootedTree([1,2,3,4], [1,1,1,1]) - - # Compute elementary weight - w = collapse_elementary_weight(tree, tdrk) - - # Manual expected weight - c = vec(sum(A1, dims=2)) - cp = vec(sum(A2, dims=2)) - b = b1 - bp = b2 - A = A1 - Ap = A2 - - w_expected = dot(b, A * (A * c)) + dot(bp, A * c) + dot(b, Ap * c) + dot(b, A * cp) + dot(bp, cp) - - @test w == w_expected - + #now test with 5 stage 7th order method from Chan and Tsai (2010) + + A_1 = [ 0 0 0 0 0; + 2//7 0 0 0 0; + 2//5 0 0 0 0; + 4//7 0 0 0 0; + 1 0 0 0 0 + ] + + b_1 = [1,0,0,0,0] + + A_2 = [ + 0 0 0 0 0; + 2//49 0 0 0 0; + 2//25 0 0 0 0; + 4//49 4//49 0 0 0; + -159//832 1715//832 -1875//832 735//832 0 + ] + + b_2 = [ + 71//960, + 2401//4800, + -625//1728, + 2401//8640, + 13//1350 + ] + + tdrk_2 = TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2) + # bseries + tdrk_series_2 = @inferred bseries(tdrk_2, 8) + #should be 7th order + @test @inferred(order_of_accuracy(tdrk_series_2)) == 7 end end # @testset "BSeries" From 12f92b0940d893ac1410fcea02cfb6d0487e8353 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 15:41:09 +0000 Subject: [PATCH 12/25] Update test/runtests.jl Co-authored-by: Hendrik Ranocha --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 59282992..333fd7a6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3261,7 +3261,7 @@ using Aqua: Aqua #constructor - tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2) + tdrk = @inferred TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2) @test tdrk isa TwoDerivativeRungeKuttaMethod @test size(tdrk.A1) == (2, 2) From e96e367924365aef8d696fade77b974a3128613b Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 15:41:18 +0000 Subject: [PATCH 13/25] Update test/runtests.jl Co-authored-by: Hendrik Ranocha --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 333fd7a6..7e444f81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3305,7 +3305,7 @@ using Aqua: Aqua 13//1350 ] - tdrk_2 = TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2) + tdrk_2 = @inferred TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2) # bseries tdrk_series_2 = @inferred bseries(tdrk_2, 8) #should be 7th order From 07c117ef65583832bb9cb5e1f01c8e1a3cb223a3 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 15:51:32 +0000 Subject: [PATCH 14/25] Add back substitute function Substitute function got deleted, not sure how all tests ran. --- src/BSeries.jl | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index 3b9f784c..f931f1b5 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1483,6 +1483,37 @@ function collapse_tree(t::RootedTree) return CollapsedArray end + +""" + substitute(b, a, t::AbstractRootedTree) + +Compute the coefficient corresponding to the tree `t` of the B-series that is +formed by substituting the B-series `b` into the B-series `a`. It is assumed +that the B-series `b` has the coefficient zero of the empty tree. + +# References + +Section 3.2 of +- Philippe Chartier, Ernst Hairer, Gilles Vilmart (2010) + Algebraic Structures of B-series. + Foundations of Computational Mathematics + [DOI: 10.1007/s10208-010-9065-1](https://doi.org/10.1007/s10208-010-9065-1) +""" +function substitute(b, a, t::AbstractRootedTree) + result = zero(first(values(a)) * first(values(b))) + + for (forest, skeleton) in PartitionIterator(t) + update = a[skeleton] + update isa Rational && iszero(update) && continue + for tree in forest + update *= b[tree] + end + result += update + end + + return result +end + """ substitute(b, a) From b8a07f48f3ba6fe6a5c9b94cca1ed0891020571d Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 15:52:19 +0000 Subject: [PATCH 15/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index f931f1b5..7852d8eb 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1272,7 +1272,7 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_i \\Delta t, y^i), \\\ u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i), \\end{aligned} -''' +``` """ struct TwoDerivativeRungeKuttaMethod{T, From 867e870c836810df08f67f2a755f084c2bc27869 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 15:52:30 +0000 Subject: [PATCH 16/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 7852d8eb..a3ac73aa 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1267,7 +1267,7 @@ autonomous problems is applied. Given an ODE ``u'(t) = f(t, u(t))`` with ``u''(t) = g(t, u(t))``, one step from ``u^{n}`` to ``u^{n+1}`` is given by -'''math +```math \\begin{aligned} y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_i \\Delta t, y^i), \\\ u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i), From fc998a39c9415f140a2e7011abadf973be43e571 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:09:53 +0000 Subject: [PATCH 17/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index a3ac73aa..bb7b6006 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1309,10 +1309,16 @@ Base.eltype(tdrk::TwoDerivativeRungeKuttaMethod{T}) where {T} = T bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) -> TruncatedBSeries Construct the truncated B-series of a two-derivative Runge–Kutta method `tdrk` -up to the specified `order`. +up to the specified integer `order`. -Returns a `TruncatedBSeries{RootedTree, V}` where `V` is inferred from -the element type of `tdrk`. +See also [`TwoDerivativeRungeKuttaMethod`](@ref). + +!!! note "Normalization by elementary differentials" + The coefficients of the B-series returned by this method need to be + multiplied by a power of the time step divided by the `symmetry` of the + rooted tree and multiplied by the corresponding elementary differential + of the input vector field ``f``. + See also [`evaluate`](@ref). """ function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) # determine coefficient type From 9de47b28f1d4833549bf75f3261ec1c78de1bdc5 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 29 Dec 2025 16:19:30 +0000 Subject: [PATCH 18/25] Add Reference for Two Derivative Method's --- src/BSeries.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/BSeries.jl b/src/BSeries.jl index bb7b6006..5eeba33b 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1274,6 +1274,13 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by \\end{aligned} ``` +# References + +-Chan, R.P.K., Tsai, A.Y.J. +"On explicit two-derivative Runge-Kutta methods." +Numer Algor 53, 171–194 (2010): +[DOI: 10.1007/s11075-009-9349-1](https://doi.org/10.1007/s11075-009-9349-1) + """ struct TwoDerivativeRungeKuttaMethod{T, MatT <: AbstractMatrix{T}, From f86592477b39f48374ca657d9cd3a0484bd7a297 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Mon, 12 Jan 2026 19:30:09 +0000 Subject: [PATCH 19/25] Refactor TwoDerivativeRungeKuttaMethod to use one c simplify code by excluding the c2 which is implicitly being computed as A_2 * e. See comment thread, for my thoughts about this. --- src/BSeries.jl | 28 ++++++++++++---------------- 1 file changed, 12 insertions(+), 16 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 5eeba33b..8f00fc93 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1287,25 +1287,23 @@ struct TwoDerivativeRungeKuttaMethod{T, VecT <: AbstractVector{T}} <:RootedTrees.AbstractTimeIntegrationMethod A1::MatT b1::VecT - c1::VecT A2::MatT b2::VecT - c2::VecT + c::VecT end -function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c1 = vec(sum(A1, dims=2)), c2 = vec(sum(A2, dims=2))) +function TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c = vec(sum(A1, dims=2))) # promote all numeric types together - T = promote_type(eltype(A1), eltype(b1), eltype(A2), eltype(b2), eltype(c1)) + T = promote_type(eltype(A1), eltype(b1), eltype(A2), eltype(b2), eltype(c)) A1T = T.(A1) b1T = T.(b1) - c1T = T.(c1) A2T = T.(A2) b2T = T.(b2) - c2T = T.(c2) + cT = T.(c) return TwoDerivativeRungeKuttaMethod{T, typeof(A1T), typeof(b1T)}( - A1T, b1T, c1T, A2T, b2T, c2T + A1T, b1T, A2T, b2T, cT ) end @@ -1393,14 +1391,13 @@ where we are evaluating eta(t) for the elementary weight of the tree t """ function derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) A1 = tdrk.A1 - c1 = tdrk.c1 A2 = tdrk.A2 - c2 = tdrk.c2 + c = tdrk.c - result1 = zero(c1) .+ one(eltype(c1)) + result1 = zero(c) .+ one(eltype(c)) if t == rootedtree(Int64[]) || t == rootedtree([1]) - return zero(c1) .+ one(eltype(c1)) + return zero(c) .+ one(eltype(c)) else subtrees_arr = subtrees(t) l = 1 @@ -1427,14 +1424,13 @@ where we are evaluating eta(t\\[1/2]) part for the elementary weight of the tree """ function collapsed_derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) A1 = tdrk.A1 - c1 = tdrk.c1 A2 = tdrk.A2 - c2 = tdrk.c2 + c = tdrk.c - result = zero(c2) + result = zero(c) if t == rootedtree(Int64[]) - return zero(c1) .+ one(eltype(c1)) + return zero(c) .+ one(eltype(c)) else collapsed_trees = collapse_tree(t) number_of_trees = length(collapsed_trees) @@ -1442,7 +1438,7 @@ function collapsed_derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKutt for k in 1:number_of_trees treecombinations = collapsed_trees[k] number2 = length(treecombinations) - sum = zero(c1) .+ one(eltype(c1)) + sum = zero(c) .+ one(eltype(c)) for m in 1:number2 step = A1 * derivative_weight(treecombinations[m], tdrk) .+ From 9f161b393fd420dd8b26bd0c1e13ee93f05a5c6d Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Sun, 18 Jan 2026 15:29:13 +0000 Subject: [PATCH 20/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 8f00fc93..df529027 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1276,10 +1276,10 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by # References --Chan, R.P.K., Tsai, A.Y.J. -"On explicit two-derivative Runge-Kutta methods." -Numer Algor 53, 171–194 (2010): -[DOI: 10.1007/s11075-009-9349-1](https://doi.org/10.1007/s11075-009-9349-1) +- Chan, R.P.K., Tsai, A.Y.J. + "On explicit two-derivative Runge-Kutta methods." + Numer Algor 53, 171–194 (2010): + [DOI: 10.1007/s11075-009-9349-1](https://doi.org/10.1007/s11075-009-9349-1) """ struct TwoDerivativeRungeKuttaMethod{T, From 52ba9fe31dbfe229dbd7556c84fab37fa3197a44 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Sun, 18 Jan 2026 16:02:06 +0000 Subject: [PATCH 21/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index df529027..260c78c6 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1269,7 +1269,7 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by ```math \\begin{aligned} - y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_i \\Delta t, y^i), \\\ + y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_j \\Delta t, y^j) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_j \\Delta t, y^j), \\\ u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i), \\end{aligned} ``` From 7f598ac425781c5df6f3ce3e27989d63326b5a8d Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Sun, 18 Jan 2026 16:02:14 +0000 Subject: [PATCH 22/25] Update src/BSeries.jl Co-authored-by: Hendrik Ranocha --- src/BSeries.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 260c78c6..9785e6c8 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1270,7 +1270,7 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by ```math \\begin{aligned} y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_j \\Delta t, y^j) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_j \\Delta t, y^j), \\\ - u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i), + u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i). \\end{aligned} ``` From 9c30d0408559de9674e906524082cfb8460898b0 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Sun, 18 Jan 2026 16:26:49 +0000 Subject: [PATCH 23/25] clean/reformat doc-strings Updated mathematical notation and descriptions for clarity. --- src/BSeries.jl | 40 +++++++++++----------------------------- 1 file changed, 11 insertions(+), 29 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 9785e6c8..f9bae247 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1268,10 +1268,10 @@ Given an ODE ``u'(t) = f(t, u(t))`` with ``u''(t) = g(t, u(t))``, one step from ``u^{n}`` to ``u^{n+1}`` is given by ```math -\\begin{aligned} - y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_j \\Delta t, y^j) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_j \\Delta t, y^j), \\\ - u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i). -\\end{aligned} +\begin{aligned} + y^i &= u^n + \Delta t \sum_j a^{1}_{i,j} f(t^n + c_j \Delta t, y^j) + \Delta t^2 \sum_j a^{2}_{i,j} g(t^n + c_j \Delta t, y^j), \\ + u^{n+1} &= u^n + \Delta t \sum_i b^1_{i} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_i b^2_{i} g(t^n + c_i \Delta t, y^i). +\end{aligned} ``` # References @@ -1349,26 +1349,12 @@ function bseries(tdrk::TwoDerivativeRungeKuttaMethod, order) end """ - elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -> Number + elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -Compute the elementary weight associated with the rooted tree `t` -for a two-derivative Runge–Kutta method `tdrk`. +Compute the elementary weight \$\Phi(t)\$ for a TDRK method. - -This follows the recursive formula for the Butcher type order conditions exhibited in, -Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. -- Numer. Algor 53, 171–194 (2010). https://doi.org/10.1007/s11075-009-9349-1 - -#see formula 16 in the paper -# alpha(t) = b1*eta(subtrees(t)) +b2*eta(collapse_trees(nu)) - -# Arguments -- `t`: A `RootedTree` representing the current term. -- `tdrk`: The `TwoDerivativeRungeKuttaMethod` whose coefficients define the - weights. - -# Returns -A scalar weight equal to the sum over all collapsed trees. +Following Chan & Tsai (2010, Eq. 16), the weight is: +\$\$ \Phi(t) = b_1 \cdot \eta(t) + b_2 \cdot \bar{\eta}(t) \$\$ """ function elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) b1 = tdrk.b1 @@ -1380,14 +1366,10 @@ end """ derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) -Compute the derivative weight for the standered trees `t` in a two-derivative Runge–Kutta (TDRK) method. - -this corresponds to formula 15 in Chan, R.P.K., Tsai, A.Y.J. On explicit two-derivative Runge-Kutta methods. - -eta(t) = A1*eta(t) + A2*eta(t/[1,2]) - -where we are evaluating eta(t) for the elementary weight of the tree t +Compute the derivative weight \$\eta(t)\$ for a TDRK method. +Following Chan & Tsai (2010, Eq. 15): +\$\$ \eta(t) = A_1 \cdot \prod \eta(t_i) + A_2 \cdot \prod \bar{\eta}(t_j) \$\$ """ function derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) A1 = tdrk.A1 From 0a057c7e458f8117c76c44b1b63ebfc5f70eaca3 Mon Sep 17 00:00:00 2001 From: John Driscoll <168301593+JohnDriscollAcademic@users.noreply.github.com> Date: Sun, 18 Jan 2026 17:37:26 +0000 Subject: [PATCH 24/25] Fix LaTeX formatting in BSeries.jl documentation This seems to work --- src/BSeries.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index f9bae247..3ebfdfd7 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1268,10 +1268,10 @@ Given an ODE ``u'(t) = f(t, u(t))`` with ``u''(t) = g(t, u(t))``, one step from ``u^{n}`` to ``u^{n+1}`` is given by ```math -\begin{aligned} - y^i &= u^n + \Delta t \sum_j a^{1}_{i,j} f(t^n + c_j \Delta t, y^j) + \Delta t^2 \sum_j a^{2}_{i,j} g(t^n + c_j \Delta t, y^j), \\ - u^{n+1} &= u^n + \Delta t \sum_i b^1_{i} f(t^n + c_i \Delta t, y^i) + \Delta t^2 \sum_i b^2_{i} g(t^n + c_i \Delta t, y^i). -\end{aligned} +\\begin{aligned} + y^i &= u^n + \\Delta t \\sum_j a^{1}_{i,j} f(t^n + c_j \\Delta t, y^j) + \\Delta t^2 \\sum_j a^{2}_{i,j} g(t^n + c_j \\Delta t, y^j), \\\\ + u^{n+1} &= u^n + \\Delta t \\sum_i b^1_{i} f(t^n + c_i \\Delta t, y^i) + \\Delta t^2 \\sum_i b^2_{i} g(t^n + c_i \\Delta t, y^i). +\\end{aligned} ``` # References From d8c7ffd64f07614ce9e23128edd38cf4bd1bcaf4 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 7 Feb 2026 14:46:01 +0100 Subject: [PATCH 25/25] fix --- src/BSeries.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/BSeries.jl b/src/BSeries.jl index 3ebfdfd7..3ed722b2 100644 --- a/src/BSeries.jl +++ b/src/BSeries.jl @@ -1276,7 +1276,7 @@ one step from ``u^{n}`` to ``u^{n+1}`` is given by # References -- Chan, R.P.K., Tsai, A.Y.J. +- Chan, R.P.K., Tsai, A.Y.J. "On explicit two-derivative Runge-Kutta methods." Numer Algor 53, 171–194 (2010): [DOI: 10.1007/s11075-009-9349-1](https://doi.org/10.1007/s11075-009-9349-1) @@ -1356,7 +1356,7 @@ Compute the elementary weight \$\Phi(t)\$ for a TDRK method. Following Chan & Tsai (2010, Eq. 16), the weight is: \$\$ \Phi(t) = b_1 \cdot \eta(t) + b_2 \cdot \bar{\eta}(t) \$\$ """ -function elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) +function RootedTrees.elementary_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) b1 = tdrk.b1 b2 = tdrk.b2 # alpha(t) = b1*eta(subtrees(t)) +b2*eta(collapse_trees(nu)) @@ -1371,7 +1371,7 @@ Compute the derivative weight \$\eta(t)\$ for a TDRK method. Following Chan & Tsai (2010, Eq. 15): \$\$ \eta(t) = A_1 \cdot \prod \eta(t_i) + A_2 \cdot \prod \bar{\eta}(t_j) \$\$ """ -function derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) +function RootedTrees.derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKuttaMethod) A1 = tdrk.A1 A2 = tdrk.A2 c = tdrk.c @@ -1408,9 +1408,9 @@ function collapsed_derivative_weight(t::RootedTree, tdrk::TwoDerivativeRungeKutt A1 = tdrk.A1 A2 = tdrk.A2 c = tdrk.c - + result = zero(c) - + if t == rootedtree(Int64[]) return zero(c) .+ one(eltype(c)) else @@ -1454,7 +1454,7 @@ function collapse_tree(t::RootedTree) subtrees_arr = subtrees(t) subtrees_multiplicity = length(subtrees_arr) - + # Recustive approach to create all the possibilities of subtrees for i in 1:subtrees_multiplicity subsubtrees = subtrees(subtrees_arr[i])