Add the capability to check whether a B-series is energy-preserving

[Sondar74]

This code is based on the Theorem 2 of the paper "Energy-Preserving Integrators and the Structure of B-series" (link: https://link.springer.com/article/10.1007/s10208-010-9073-1).
In this code we want to know if a method is Energy-Preserving or not up to some order.
The function 'EnergyPreserving(A,b,s)' receives as input a square Matrix which is the inner part of the Butcher Array of a numerical method, a vector b (which length is the same as the amount of rows of A) and an intenger s, which will be the order up to which we are testing this method.

For example, the classical Runge-Kutta method:

A = [0//1  0//1  0//1  0//1
 1//2  0//1  0//1  0//1
 0//1  1//2  0//1  0//1
 0//1  0//1  1//1  0//1]
 b =  [1//6, 1//3, 1//3, 1//6]

the function 'EnergyPreserving(A,b,s)' will return 'True' for s = 3,4 but will return 'False' for s = 5, i.e. 'EnergyPreserving(A,b,5)'.

Furthermore, for any method described by the square matrix A and vector b, the function 'OrderMethod' will say up to which order it is energy-preserving. For example, for the same (A,b) as before this function will return "Energy Preserving for order < 4".

Finally, if you have already generated a BSeries serie following the next steps:

 rka = RungeKuttaMethod(A, b)
#generate bseries 
 series = bseries(rka, s)

the function is_energy_preserving(serie) will return True for s = 4, and False for s = 5.

[ketch]

Thanks, @Sondar74 ! Could you edit the box above (where it says "no description provided") and add a description of the contents of this PR? Specifically:

• Describe what functionality is added
• Give example(s) of how to use it

[ketch]

Also, it would be great to make a notebook like https://github.com/ranocha/BSeries.jl/blob/main/docs/src/tutorials/bseries_basics.ipynb but explaining and demonstrating the functionality of the code in this PR. From the Jupyter notebook we can generate a page for the online docs that will look similar to https://ranocha.de/BSeries.jl/dev/tutorials/bseries_basics/.

[ranocha]

Thanks for the PR! As @ketch noted, we have a bit of work to do before we can get this into BSeries.jl. Some general comments (besides the aspects already mentioned by @ketch)

• Please follow standard Julia style, e.g., snake_case for functions and variables. See https://docs.julialang.org/en/v1/manual/style-guide/
• Please add comments explaining why you do something and not only what you do
• Please add references to the basic algorithms that you apply
• Please add docstrings

[Sondar74]

The Style Guide site recommends not to use the from x -> f(x). Should I rewrite this lines?
last_j_occurrence = findlast(x -> x == j, a)

[ranocha]

No, that's fine as it is

[ranocha]

I haven't looked at the algorithm in details yet. How does it compare to the approach of Sundklakk (2015)?

[Sondar74]

I haven't looked at the algorithm in details yet. How does it compare to the approach of Sundklakk (2015)?

I was trying to download the paper from this website, is this the reference you mentioned? And also, can you share the paper with me?

https://www.semanticscholar.org/paper/A-Library-for-Computing-with-Trees-and-B-Series-Sundklakk/4c5c94924f63ed3c61997c5972955d20f83a2145

[ranocha]

Done on Slack

[ranocha]

Thanks! The current version is already much more efficient 👍

Next, there are still some type instabilities and further possibilities for improvements. You can get an impression of the overall performance profile via

julia> using ProfileView, BSeries

julia> series = bseries(10) do t, series
           if order(t) in (0, 1)
               return 1 // 1
           else
               v = 1 // 1
               n = 0
               for subtree in SubtreeIterator(t)
                   v *= series[subtree]
                   n += 1
               end
               return v / (n + 1)
           end
       end

julia> @profview is_energy_preserving(series) # ignore first profile including compilation

julia> @profview is_energy_preserving(series) # use second profile after compilation

This yields

for me. Some takeaways are

1. A substantial amount of time (roughly 50%) is spent in modified_equation - which is fine.
2. You can see that there are type instabilities - the red bars. They are for example cause by lines such as same_order_trees = [] and same_order_coeffs = [], where you create Vector{Any}, fill them, and pass them to another function. Something like this should be avoided by specifying the types concretely (based on the coefficients you get).
3. You use dense linear algebra. I think it may be worth a try to use sparse linear algebra instead. Then, you should also avoid collecting all the column vectors by just collecting the relevant values and indices while assembling the sparse matrix.
4. You call rank twice. Depending on whether this is a bottleneck at the end, I wonder whether there is a clever way to avoid this using the SVD or something similar directly.

[ranocha]

I will make some fixes to this PR and perform a final review before merging the first working version

[ranocha]

Thanks a lot!