Implementing CSRK code in BSeries.jl #144

Sondar74 · 2023-06-26T12:27:02Z

This code is based on the paper "A characterization of energy-preserving methods and the construction of parallel integrators for Hamiltonian systems" (https://doi.org/10.48550/arXiv.1505.02537).

Added a struct in BSeries.jl called ContinuousStageRungeKutta and a function CSRK().

Example:

The energy-preserving 4x4 matrix given by Miyatake & Butcher (2015) is

M = [-6//5   72//5  -36//1  24//1;
        72//5  -144//5 -48//1  72//1;
        -36//1 -48//1   720//1 -720//1;
        24//1   72//1  -720//1  720//1]

Then, we calculate the bseries with the following code:

csrk = ContinuousStageRungeKuttaMethod(M)
series = bseries(csrk, 4)
TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entries:
  RootedTree{Int64}: Int64[]      => 1//1
  RootedTree{Int64}: [1]          => 1//1
  RootedTree{Int64}: [1, 2]       => 1//2
  RootedTree{Int64}: [1, 2, 3]    => 1//6
  RootedTree{Int64}: [1, 2, 2]    => 1//3
  RootedTree{Int64}: [1, 2, 3, 4] => 1//24
  RootedTree{Int64}: [1, 2, 3, 3] => 1//12
  RootedTree{Int64}: [1, 2, 3, 2] => 1//8
  RootedTree{Int64}: [1, 2, 2, 2] => 1//4

…should work better than the others.

ranocha · 2023-06-26T13:19:08Z

The failing tests come from a deprecation warning fixed in JuliaPy/PyCall.jl#1042. It will be fixed by #145

ranocha

Thanks a lot for the PR! Now, everything seems to work. Please close the obsolete PRs.

Please add references and explain why you are doing what you are doing.
Please add tests.
Please check type stability (also in the tests).
Could you please describe an example of the new functionality in your opening post?
How is the performance behavior compared to the other implementations generating B-series we have so far?

src/BSeries.jl

ranocha · 2023-06-26T13:37:33Z

src/BSeries.jl

+    # we need variables to work with
+    variable1 = Sym(t)
+    variable2 = Sym(z)


I do not really like that this introduces a hard dependency on SymPy.jl here (which is not covered in the import statements at the top of this file or Project.toml). If I remember correctly, @ketch reported some installation issues with some Julia packages build on Python packages. If possible, I would like to avoid such issues.

Since we work only with polynomials and just need to integrate, in principle everything can be done by working directly with the coefficients. This will make the code a bit harder to read but may also be faster.

…differentials_csrk'

src/BSeries.jl

ranocha · 2023-06-27T12:43:37Z

src/BSeries.jl

+function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
+    csrk = csrk.matrix
+    bseries(o) do t, series
+        if order(t) in (0, 1)
+            return one(csrk)
+        else
+            return elementary_differentials_csrk(csrk, t)
+        end
+    end
+end


Please add tests. This does not work correctly right now.

Okey, I'm working on this.

When working with the function one(), I have the problem that eltype(csrk) returns Any. How can I fix it?

What is the csrk that you use here?

It is the csrk obtained from csrk = ContinuousStageRungeKuttaMethod(M), for a matrix M.

Use eltype(csrk.matrix) instead or define eltype(csrk::ContinuousStageRungeKuttaMethod) appropriately. The latter would be my preferred option, e.g.,

struct ContinuousStageRungeKuttaMethod{T, MatT <: AbstractMatrix{T}} <: RootedTrees.AbstractTimeIntegrationMethod A::MatT end Base.eltype(::ContinuousStageRungeKuttaMethod{T}) where {T} = T

ranocha

Please add tests with floating point coefficients.
Please add tests with symbolic coefficients using SymPy.jl, SymEngine.jl, and Symbolics.jl.

src/BSeries.jl

ranocha · 2023-06-28T04:53:33Z

src/BSeries.jl

+  RootedTree{Int64}: [1, 2, 3]    => 6004799503160661//36028797018963968
+  RootedTree{Int64}: [1, 2, 2]    => 6004799503160661//18014398509481984
+  RootedTree{Int64}: [1, 2, 3, 4] => 6004799503160661//144115188075855872
+  RootedTree{Int64}: [1, 2, 3, 3] => 6004799503160661//72057594037927936


Is this really correct?

No, sorry: seems like there is something happening with the definition of V in

series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}() series[rootedtree(Int[])] = one(V)

so that, although the output of elementary_differentials_csrk(csrk, t) is correct, for some values it will throw an approximation to the original result instead of saving the value of elementary_differentials_csrk(csrk, t) Do you know how can I fix it?

Seems like the numerical function N() solves this issue, but I don't know if it has limitations for future applications of this code (since this is the first time I'm working with this function).

Let's see. First, please add more tests etc. - and fix the CI problems. Then, possible problem can show up already

src/BSeries.jl

test/runtests.jl

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

…kupdated

Sondar74 · 2023-06-28T10:31:28Z

Please add tests with floating point coefficients.

Please add tests with symbolic coefficients using SymPy.jl, SymEngine.jl, and Symbolics.jl.

Seems like a symbolic matrix with inputs from SymEngine is not compatible with the PolynomialsA function.

…and SymEngine

ranocha · 2023-06-28T13:09:30Z

Please add tests with floating point coefficients.

Please add tests with symbolic coefficients using SymPy.jl, SymEngine.jl, and Symbolics.jl.

Seems like a symbolic matrix with inputs from SymEngine is not compatible with the PolynomialsA function.

That doesn't surprise me since SymPy.jl is probably not compatible with SymEngine.jl etc.

test/runtests.jl

ranocha · 2023-06-28T13:20:13Z

test/runtests.jl

+            # Create symbolic matrix
+            Minv = [a 1//2 b 1//6;
+            1//2 1//4 1//6 1//8;
+            b 1//6 c 1//10;
+            1//6 1//8 1//10 1//12]


Please check your comments everywhere. This is not a symbolic matrix since you set the coefficients above to fixed floating point values. Where did you get these special choices from? Please add some references

Okey, I'll fix the comments. As for the values, I took them from the project I am working on with @ketch, I calculated them.
Since it is the same matrix from Miyatake & Butcher, I will add this reference too.

ranocha · 2023-06-28T13:20:42Z

test/runtests.jl

+            # Generate the bseries up to order 4
+            order = 4
+            series = bseries(csrk, order)
+            expected_coefficients = [1.0 ,1.0, 0.5000000000000002, 0.16666666666666666, 0.3333333333333336, 0.04166666666666661, 0.0833333333333332, 0.12500000000000003, 0.2500000000000003]
+            @test collect(values(series)) == expected_coefficients


Same comment about the order, energy-preserving property etc. as above

I think the is_energy_preserving function will not necesary work accurately for Floating coefficients: for example, if we wanted to express 1/3 in floating we may need to cut it in some point at some decimal place, so there would be loss of information and is_energy_preserving would return false

Did you test it?

test/runtests.jl

ranocha · 2023-06-28T13:22:21Z

test/runtests.jl

+            # Generate the bseries up to order 4
+            order = 3
+            series = bseries(csrk, order)
+            expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]


Did you compute this by hand? How did you arrive at these coefficients?

No, I just copied from the output of the working code

Hmm... So it is just a regression test - assuming that the current implementation is correct. Could you please also add real tests using examples with known behavior? IIRC, Miyatake & Butcher (2015) have some examples with free coefficients in their paper. You could express these free coefficients as symbolic variables and check whether the resulting B-series of the CSRK method

has the correct order_of_accuracy

is_energy_preserving

test/runtests.jl

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

…kupdated

Gabriel and others added 6 commits May 17, 2023 11:13

comentario

7853dcb

Merge branch 'main' into Energy_Preserving

aaf5940

Merge branch 'main' of github.com:Sondar74/BSeries.jl into main

afe7373

updated 'main' branch and created new branch 'csrkupdated'. This one …

e4d79cc

…should work better than the others.

deleted extra file

823684d

fixing 'Spell Check'

0df736b

ranocha requested changes Jun 26, 2023

View reviewed changes

Sondar74 added 3 commits June 27, 2023 11:15

added docstrings and removed elementary_differentials_csrk from export

09d4629

modified the comments in the functions 'PolynomialA' and 'elementary_…

335265f

…differentials_csrk'

Modified the function 'bseries' for CSRK

8c4e0b1

ranocha requested changes Jun 27, 2023

View reviewed changes

Sondar74 added 3 commits June 27, 2023 15:47

exported 'ContinuousStageRungeKuttaMethod'

7ad5325

added one test, and modified the function 'bseries' for CSRK

5d82c87

Spell Check

c0a724d

ranocha requested changes Jun 28, 2023

View reviewed changes

Sondar74 and others added 8 commits June 28, 2023 10:33

Update src/BSeries.jl

c45de1c

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

Update src/BSeries.jl

6616b9c

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

Update src/BSeries.jl

1227e1a

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

Update src/BSeries.jl

3691f69

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

Update test/runtests.jl

68cecb0

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

updating docstring

5d30865

Merge branch 'csrkupdated' of github.com:Sondar74/BSeries.jl into csr…

7387ffc

…kupdated

fixed the values of bseries()

53fd493

Sondar74 added 2 commits June 28, 2023 13:34

added test for floating-point and symbolic coefficientes using SymPy …

d6a23d6

…and SymEngine

added test for Symbolics.jl: same as SymEngine

c0ba86f

ranocha requested changes Jun 28, 2023

View reviewed changes

Update test/runtests.jl

b5813d2

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

Sondar74 and others added 5 commits June 28, 2023 17:06

Update test/runtests.jl

691507f

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

Update test/runtests.jl

02acff4

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

added Energy_Preserving test

fd1a23e

Merge branch 'csrkupdated' of github.com:Sondar74/BSeries.jl into csr…

0c29499

…kupdated

fixed the funtion bseries(csrk, order)

c8113e6

ranocha mentioned this pull request Jul 22, 2023

first version of continuous stage Runge-Kutta methods #151

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implementing CSRK code in BSeries.jl #144

Implementing CSRK code in BSeries.jl #144

Sondar74 commented Jun 26, 2023 •

edited

Loading

ranocha commented Jun 26, 2023 •

edited

Loading

ranocha left a comment

ranocha Jun 26, 2023

ketch Jun 26, 2023

ranocha Jun 27, 2023

Sondar74 Jun 27, 2023

Sondar74 Jun 27, 2023

ranocha Jun 28, 2023

Sondar74 Jun 28, 2023

ranocha Jun 28, 2023 •

edited

Loading

ranocha left a comment

ranocha Jun 28, 2023

Sondar74 Jun 28, 2023

Sondar74 Jun 28, 2023

ranocha Jun 28, 2023 •

edited

Loading

Sondar74 commented Jun 28, 2023

ranocha commented Jun 28, 2023

ranocha Jun 28, 2023

Sondar74 Jun 28, 2023 •

edited

Loading

ranocha Jun 28, 2023

Sondar74 Jun 28, 2023

ranocha Jun 28, 2023

ranocha Jun 28, 2023

Sondar74 Jun 28, 2023

ranocha Jun 28, 2023

Implementing CSRK code in BSeries.jl #144

Are you sure you want to change the base?

Implementing CSRK code in BSeries.jl #144

Conversation

Sondar74 commented Jun 26, 2023 • edited Loading

Example:

ranocha commented Jun 26, 2023 • edited Loading

ranocha left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ranocha Jun 28, 2023 • edited Loading

Choose a reason for hiding this comment

ranocha left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ranocha Jun 28, 2023 • edited Loading

Choose a reason for hiding this comment

Sondar74 commented Jun 28, 2023

ranocha commented Jun 28, 2023

Choose a reason for hiding this comment

Sondar74 Jun 28, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Sondar74 commented Jun 26, 2023 •

edited

Loading

ranocha commented Jun 26, 2023 •

edited

Loading

ranocha Jun 28, 2023 •

edited

Loading

ranocha Jun 28, 2023 •

edited

Loading

Sondar74 Jun 28, 2023 •

edited

Loading