Add support for Two-Derivative RK Methods #299

JohnDriscollAcademic · 2025-11-12T02:20:17Z

Adds support for two-derivative and additive two-derivative Runge-Kutta methods, including the necessary structures and functions to work with B-series. Note I defined the structures directly where the existing methods seem to import them from RootedTrees.jl.

Implement TwoDerivativeRungeKuttaMethod and TwoDerivativeAdditiveRungeKuttaMethod structures with associated functions for B-series computations.

Removed unnecessary blank lines and comments in BSeries.jl.

ranocha · 2025-11-21T10:15:42Z

src/BSeries.jl

+#####################################################################
+# defining the new methods, not sure where best place for this code is.


Thanks a lot for your contribution! I have a couple of suggestions/questions to improve the code before we merge it into the main branch.

Could you please focus first on one of the two new types, e.g., only on the standard two-derivative RK methods and not the additive ones? This allows us to improve the code of this version first so that we don't have the same discussions twice.

Could you please keep the code for the TwoDerivativeRungeKuttaMethods together? Right now, you inserted some parts between the code belonging to the ContinuousStageRungeKuttaMethods. This is a bit confusing for people reading the code (and code is always read more often than it is written). A good place to insert this code could be, e.g.,

BSeries.jl/src/BSeries.jl

Line 1122 in 5183752

or

BSeries.jl/src/BSeries.jl

Line 1255 in 5183752

Could you please add some tests (to https://github.com/ranocha/BSeries.jl/blob/main/test/runtests.jl) to ensure that everything works well and keeps working well when the code develops?

Thanks for the suggestions. I have cut back to just focusing on the standard two-derivative methods for this pr, put all the new code in one section, and added a test. I am not very experienced with contributing on GitHub, so I have changed this to a "Draft" Pull Request in case more iteration is needed.

Regarding the test, so far I have just had it use a known method from the paper by Robert Chan, and had it check that the proper TwoDerivativeRungeKuttaMethod struct was created, that the B-series computed the correct order of accuracy for the method, and checked a specific tree weight. This seems to be in line with what other tests have done, but please let me know if anything else should be included.

src/BSeries.jl

Added tests for the TwoDerivativeRungeKuttaMethod including constructor validation and accuracy checks.

Placed all new functionality in one section. This adds support for multi-derivative RK methods and defines these methods explicitly from two tableau

ranocha

Thanks! It looks like your approach will already allocate some temporary arrays. Could you please compare the performance and results with an approach creating the B-series using compose as in https://ranocha.de/BSeries.jl/stable/tutorials/bseries_creation/#tutorial-bseries-creation-RK-compose?

This should work directly for an explicit method. For an implicit method, it should work when used as a fixed point map with enough iterations so that the result does not change anymore up to the desired order. This approach is used for Runge-Kutta methods in Section 313 of Butcher's book "Numerical Methods for Ordinary Differential Equations" to derive the order conditions of RK methods.

src/BSeries.jl

ranocha · 2025-11-24T11:01:43Z

src/BSeries.jl

+up to the specified `order`.
+
+Returns a `TruncatedBSeries{RootedTree, V}` where `V` is inferred from
+the element type of `tdrk`.


Suggested change

up to the specified `order`.

Returns a `TruncatedBSeries{RootedTree, V}` where `V` is inferred from

the element type of `tdrk`.

up to the specified integer `order`.

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).

src/BSeries.jl

ranocha · 2025-11-24T11:09:05Z

src/BSeries.jl

+        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


I am not sure whether I like this approach. A two-derivative RK method should have a B-series expansion in terms of standard rooted trees, not colored trees.

I agree. In this commit the series itself is using only standard rooted trees; I'm converting to colored rooted trees only to compute the weights. This is because when I was writing the collapsing function I used for weight computation, I was planning ahead for when we want to also collapse trees in an additive method. I could presumably create a slightly simpler function to handle only the multi-derivative case, but then I wouldn’t be able to reuse the function for the additive multi-derivative trees.

I think this will again boil down to whether the way I'm computing the tree weights is significantly more computationally expensive than the recursive or compose methods, which I plan to test soon.

ranocha · 2025-11-24T11:12:39Z

test/runtests.jl

+    @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;


Can you please add more methods to the tests, in particular some higher-order methods from their paper? Some special cases of trees start at order 5 and higher.

I thought I had checked my logic and all the trees many times, but alas, I found a mistake in one of my fifth-order trees that causes a miscount and leads to an error in the order. Thank you for the suggestion to check the fifth order!! I’m going to see if there’s a quick fix for this before continuing with anything else.

I knew my implementation was far from memory optimized but I did not think it was this bad. I did some quick benchmarking and here are the results. For my method, collapse all trees, then find order conditions from there, I get the following where I just test the time it takes to generate the B-series for a given order:

Computing the B-series coefficients up to order: 7 0.009580 seconds (291.81 k allocations: 14.426 MiB) 0.018783 seconds (291.64 k allocations: 14.417 MiB, 49.28% gc time) Computing the B-series coefficients up to order: 8 0.061796 seconds (1.44 M allocations: 72.963 MiB, 15.03% gc time) 0.054118 seconds (1.44 M allocations: 72.963 MiB, 13.11% gc time) Computing the B-series coefficients up to order: 9 0.295879 seconds (7.27 M allocations: 379.771 MiB, 18.92% gc time) 0.290079 seconds (7.27 M allocations: 379.771 MiB, 16.50% gc time)

where I’ve run twice to account for Julia's compilation time.

Calculating differentials the more streamlined way, that is as recursively defined in the paper and implemented in the thesis you shared with me, we get far nicer results, polynomial rather than exponential!!

Computing the B-series coefficients up to order: 7 0.005005 seconds (76.89 k allocations: 3.987 MiB) 0.004866 seconds (76.86 k allocations: 3.985 MiB) Computing the B-series coefficients up to order: 8 0.020274 seconds (272.40 k allocations: 14.138 MiB) 0.019488 seconds (272.40 k allocations: 14.138 MiB) Computing the B-series coefficients up to order: 9 0.087649 seconds (968.06 k allocations: 50.192 MiB, 28.02% gc time) 0.063105 seconds (968.06 k allocations: 50.192 MiB)

This shown, I think the recursive method is obviously what I should have implemented. Why I originally veered away from it was my goal was generating order conditions for the low-stage IMEX multi-derivative method (which now I know is possible only since we went to order 4 :). It's not immediately clear to me how to generate the weights recursively for the additive case, however with some thinking I’m sure we could figure it out. Additionally, doing things by generating trees also allowed me to quickly reduce the order conditions in the case that either component was linear; it's also not clear to me how easy this would be to implement in the recursive case.

All that said, if you would like to move forward with multi-derivative methods, I could push a version using the recursive definition, which I just tested as agreeing with the 5s7p method. That said, my implementation is almost line for line with that of the thesis, so if we were to merge something like that I would hope that student would be properly credited. As for multi-derivative additive methods, I would need to do some more thinking before reimplementing that in the memory-friendly manner.

Thanks a lot for performing these benchmarks! Since we aim for an efficient method that is able to handle trees also up to order 9 and so, I would prefer the more efficient method you mentioned. Could you please update the PR accordingly?

And sorry for the late reply...

Additionally, doing things by generating trees also allowed me to quickly reduce the order conditions in the case that either component was linear; it's also not clear to me how easy this would be to implement in the recursive case.

You could first generate the full B-series and then discard all trees where the elementary differentials are zero for linear problems. That's how I would approach this problem.

Apologies for the late reply, things got busy with our Thanksgiving holiday followed by the start of our exam period. I have updated the PR to include the more efficient method and added a test for a method up to order 5 from the referenced paper. I'll try benchmarking the current implementation again; while it certainly outperformed my original implementation, it will be nice to see how it compares to something like the standard or additive methods existing.

That would be nice indeed 👍

add comments Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

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.

Add a test for multi-derivative Runge-Kutta methods, specifically a 5s7p method from Chan and Tsai (2010).

ranocha · 2025-12-19T13:30:20Z

src/BSeries.jl

-function substitute(b, a, t::AbstractRootedTree)
-    result = zero(first(values(a)) * first(values(b)))
+struct TwoDerivativeRungeKuttaMethod{T,
+                                     MatT <: AbstractMatrix{T},
+                                     VecT <: AbstractVector{T}} <:RootedTrees.AbstractTimeIntegrationMethod


It looks like the method substitute(b, a, t::AbstractRootedTree) has been deleted here?

ranocha · 2025-12-19T13:31:26Z

src/BSeries.jl

+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


Suggested change

'''math

```math

ranocha · 2025-12-19T13:31:38Z

src/BSeries.jl

+  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}
+'''


Suggested change

'''

```

ranocha · 2025-12-19T13:32:15Z

src/BSeries.jl

-    substitute(b, a, t::AbstractRootedTree)
+    TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2, c = vec(sum(A1, dims=2)))

-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.
+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.


Can you please add some references (following the format we have used before) to the end of the docstring?

ranocha · 2025-12-19T13:32:51Z

src/BSeries.jl

+    b1::VecT
+    c1::VecT
+    A2::MatT
+    b2::VecT
+    c2::VecT


The documentation says we have only

Suggested change

b1::VecT

c1::VecT

A2::MatT

b2::VecT

c2::VecT

b1::VecT

c::VecT

A2::MatT

b2::VecT

ranocha · 2025-12-19T13:34:40Z

test/runtests.jl

+
+        #constructor
+
+        tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)


Suggested change

tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)

tdrk = @inferred TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)

To check type stability

ranocha · 2025-12-19T13:35:05Z

test/runtests.jl

+            13//1350
+        ]
+
+        tdrk_2 = TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2)


Suggested change

tdrk_2 = TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2)

tdrk_2 = @inferred TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2)

JohnDriscollAcademic added 3 commits November 12, 2025 02:11

Add support for Two-Derivative and Addative Two-Derivative RK Methods

4448537

Implement TwoDerivativeRungeKuttaMethod and TwoDerivativeAdditiveRungeKuttaMethod structures with associated functions for B-series computations.

Fix typos in BSeries.jl

0f00454

Clean up BSeries.jl by removing blank lines

c0d4d78

Removed unnecessary blank lines and comments in BSeries.jl.

ranocha requested changes Nov 21, 2025

View reviewed changes

JohnDriscollAcademic marked this pull request as draft November 23, 2025 18:32

JohnDriscollAcademic changed the title ~~Add support for Two-Derivative and Addative Two-Derivative RK Methods~~ Add support for Two-Derivative RK Methods Nov 24, 2025

JohnDriscollAcademic added 2 commits November 24, 2025 00:12

Implement tests for TwoDerivativeRungeKuttaMethod

e5886fc

Added tests for the TwoDerivativeRungeKuttaMethod including constructor validation and accuracy checks.

Condensed all multi-derivative functionality into a single section

9ea2fcb

Placed all new functionality in one section. This adds support for multi-derivative RK methods and defines these methods explicitly from two tableau

ranocha reviewed Nov 24, 2025

View reviewed changes

JohnDriscollAcademic and others added 6 commits November 24, 2025 15:13

Update src/BSeries.jl

5e72636

add comments Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Update src/BSeries.jl

fe35120

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Update src/BSeries.jl

bb3dcb0

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Update src/BSeries.jl

fef7a7d

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

Update tests for TwoDerivativeRungeKuttaMethod

40e8f9e

Add a test for multi-derivative Runge-Kutta methods, specifically a 5s7p method from Chan and Tsai (2010).

ranocha requested changes Dec 19, 2025

View reviewed changes

		#####################################################################
		# defining the new methods, not sure where best place for this code is.

-up to the specified `order`.
-Returns a `TruncatedBSeries{RootedTree, V}` where `V` is inferred from
-the element type of `tdrk`.
+up to the specified integer `order`.
+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).


		#constructor

		tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)

	tdrk = TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)
	tdrk = @inferred TwoDerivativeRungeKuttaMethod(A1, b1, A2, b2)

	tdrk_2 = TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2)
	tdrk_2 = @inferred TwoDerivativeRungeKuttaMethod(A_1, b_1, A_2, b_2)

Add support for Two-Derivative RK Methods #299

Are you sure you want to change the base?

Add support for Two-Derivative RK Methods #299

Uh oh!

Conversation

JohnDriscollAcademic commented Nov 12, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

ranocha left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants