Adding Symplectic Integrators

[JacobHass8]

Adds symplectic solvers for ODE systems with a conserved quantity (i.e. energy). This was a requested scipy feature (see #303 and scipy/scipy#12690). I'd ultimately like to merge this into scipy using cython. I'm still working on this and have a couple of features I'd like to add:

• 11th order method from here
• Tests
  • Implement Hénon–Heiles potential
  • Add floating point tests
• Add multiprecision support
• Update documentation for singleton types
• Driver function which allows users to run a ode to a specific time with a specific method
• Helper function for 1d integrator
• Examples
• Make sure the function signature is similar to other methods in boost

Does this seem like a good plan? I'd appreciate any input.

[NAThompson]

@JacobHass8 : I think a good test is to identify the conserved quantity you wish to be conserved, and show that quantity is much better preserved with a symplectic integrator than with (say) RK4.

I would also recommend an API for event detection and return a solution skeleton that can be interpolated as a Hermite spline, i.e., return ${t_k, y_k, dot{y}k}{k=0}^{n-1}$ rather than the typical solution skeleton {t_k, y_k, }_{k=0}^{n-1}.

[JacobHass8]

@JacobHass8 : I think a good test is to identify the conserved quantity you wish to be conserved, and show that quantity is much better preserved with a symplectic integrator than with (say) RK4.

I've tried this on a harmonic oscillator and the energy fluctuations seem to be of order 1e-11 (for the 6th order method). I haven't checked RK4 (or any other method though!).

I would also recommend an API for event detection and return a solution skeleton that can be interpolated as a Hermite spline, i.e., return ${t_k, y_k, dot{y}k}{k=0}^{n-1}$ rather than the typical solution skeleton {t_k, y_k, }_{k=0}^{n-1}.

I'm not sure I totally understand what you mean. Are you saying return a third object that could be used to interpolate between different $t_k$?

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.41%. Comparing base (8ee12a5) to head (4d11732).
⚠️ Report is 19 commits behind head on develop.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1411      +/-   ##
===========================================
+ Coverage    95.40%   95.41%   +0.01%     
===========================================
  Files          828      831       +3     
  Lines        69072    69405     +333     
===========================================
+ Hits         65896    66226     +330     
- Misses        3176     3179       +3     

Files with missing lines	Coverage Δ
include/boost/math/quadrature/symplectic.hpp	100.00% <100.00%> (ø)
test/test_symplectic.cpp	100.00% <100.00%> (ø)

... and 4 files with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 8ee12a5...4d11732. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[NAThompson]

Are you saying return a third object that could be used to interpolate between different t_k?

Exactly. Say you have an $\mathcal{O}(\delta t^6)$ accurate ODE stepper. That means you only need to store a very sparse set of $(t_k, y_k)$ pairs. And each ODE stepper does come with a "natural interpolant", so naively this is no loss of information. But the user gets this data back and has no clue how to interpolate it, or what the "natural interpolant" even is. So they resort to just using linear interpolation. This makes the display and subsequent use of the data worse than if you used a less accurate ODE stepper.

However, you have to compute $\dot{y}_k = f(t_k, y_k)$ during the course of the computation anyway, so if you just store it you can at least give the user a way to display the data as a hermite spline. I've been trying to get this built up in H5Web-see here.

Note that this is still inferior to the "natural interpolant", but given there seems to be no hope to get every ODE stepper's interpolant into the standard graphics packages, I think this is a reasonable compromise.

[JacobHass8]

The methods implemented here do appear better than those in Scipy. I test a simple harmonic oscillator where the Hamiltonian is given by $H = p^2 + x^2$. Starting with the initial condition $x=1$ and $p=0$, the total energy of the system is $1$. I ran the system from time 0 to 1 with time steps of 0.05. Here's a plot of the position as a function of time. Both graphs match up pretty nicely so we know they are both performing as expected.

However, if we look $p^2 + x^2 -1$, which I call the energy loss, it's clear that the symplectic integrators implemented here outperform scipy's Runge-Kutta method.

[NAThompson]

Should we do an enum of available methods rather than a vector of strings?

[JacobHass8]

I wasn't quite sure how to implement the enum of available methods. I'm just not very familiar with enum types. Instead, I used a map. Is this okay or is there an easier way with enum?

[NAThompson]

Actually, you're targeting a scipy wrap correct? So this question should be answered in the context of "what is the most expressive way to do this in python via nanobind?"

[JacobHass8]

Scipy actually uses cython to bind c++. I think strings are probably okay for now. This is the same behavior as scipy.integrate.solve_ivp.

[mborland]

Also highly recommend going the enum route instead of a vector of strings and then a map.

enum class available_methods : unsigned 
{
    Y6 = 0,
    Y4 = 1, 
    Y2 = 2 
}

And then you update the signature:

template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian_imp(const RandomAccessContainer p0,
                                                                                                             const RandomAccessContainer q0,
                                                                                                             const RealType dt,
                                                                                                             const unsigned steps,
                                                                                                             Func dHdp,
                                                                                                             Func dHdq,
                                                                                                             available_methods method,
                                                                                                             const Policy& pol)

[JacobHass8]

@NAThompson I created a .qbk file to add documentation for the functions I added. I couldn't figure out how to build the documentation though. I kept getting errors when trying to follow the boost math guide.

I wanted to include an equation so I created a .svg and .png using latex. However, I noticed that all other equations have a .mml file. Do you know how to convert these?

[NAThompson]

I wanted to include an equation so I created a .svg and .png using latex. However, I noticed that all other equations have a .mml file. Do you know how to convert these?

I just use this website and drop pngs in it.

I also experienced errors while building the docs-I think @mborland may have to help based on your particular error.

[jzmaddock]

I created a .qbk file to add documentation for the functions I added. I couldn't figure out how to build the documentation though. I kept getting errors when trying to follow the boost math guide.

What are the errors you are getting?

Unfortunately there is a rather long toolchain required for building the docs, we did have it all documented on the old wiki, but I have a hunch that's bitten the dust :(

I wanted to include an equation so I created a .svg and .png using latex. However, I noticed that all other equations have a .mml file. Do you know how to convert these?

Don't worry about .mml, these are all a "historical artefact": when we started there was no easy way to convert LaTex to SVG or something web friendly (other than blocky pngs) but MML was the new kid on the block and there was a nice GUI editor and some scripts to convert the mml to SVG so that was what we used. Nowadays those tools are all unmaintained, and LaTex is so much better... please just make sure that the original LaTex is archived somewhere (a comment in the docs next to the link including the SVG would be good) that way if someone later spots an error it's so much easier/quicker to fix!

[JacobHass8]

What are the errors you are getting?

Unfortunately there is a rather long toolchain required for building the docs, we did have it all documented on the old wiki, but I have a hunch that's bitten the dust :(

Actually, no errors. I found the .xml file but don't know how to convert this to html. I just ran b2 in the doc directory.

Don't worry about .mml, these are all a "historical artefact": when we started there was no easy way to convert LaTex to SVG or something web friendly (other than blocky pngs) but MML was the new kid on the block and there was a nice GUI editor and some scripts to convert the mml to SVG so that was what we used. Nowadays those tools are all unmaintained, and LaTex is so much better... please just make sure that the original LaTex is archived somewhere (a comment in the docs next to the link including the SVG would be good) that way if someone later spots an error it's so much easier/quicker to fix!

Sounds good! I added the latex command.

[JacobHass8]

@jzmaddock or @NAThompson do you think Kahan summation might help the accuracy of these algorithms? Essentially, the algorithm is just summing together a bunch of small steps so I thought it might help.

Edit: I implemented kahan summation, but it didn't do anything. I could have done it wrong, but I don't think that's the case. I'm only testing against the simple harmonic oscillator example. Thus, the error term just oscillates from a small negative to a small positive value.

[NAThompson]

@JacobHass8 : I personally don't see any reason Kahan summation is going to help here. I would try to compute the condition number of the ODE using the techniques in Corless, A Graduate Introduction to Numerical Methods. This will tell you whether it's the ODE itself that forces the inaccuracy, or it's the method that forces inaccuracy on the ODE.

[JacobHass8]

@JacobHass8 : I personally don't see any reason Kahan summation is going to help here. I would try to compute the condition number of the ODE using the techniques in Corless, A Graduate Introduction to Numerical Methods. This will tell you whether it's the ODE itself that forces the inaccuracy, or it's the method that forces inaccuracy on the ODE.

Unfortunately, I haven't been able to find a copy of that book. I think I'll leave Kahan summation out for the first implementation. After this ports to scipy maybe someone will find need for it but I'll leave it out for now.

[JacobHass8]

On another note, is there any need to do policy based up/down casting on the input/result? From my understanding, boost policy casting is only implemented numeric types and not containers.

[JacobHass8]

@NAThompson, @jzmaddock, @mborland if all the tests come back passing, I think this PR is ready for review. There are only two small issues that I can think of .

First, if the timestep, dt, is less than 0 an error is thrown according to the policy. However, if the policy wants to silently return an nan, this isn't handled. The algorithm just proceeds as normally. I don't know why anyone would want to run this without throwing an error though.

[JacobHass8]

@NAThompson, @jzmaddock or @mborland gentle ping to see if anyone would be interested in reviewing. I can also make plots of more systems with my algorithm vs scipy if that would be helpful.

[mborland]

Your accuracy is clearly good based on the plots in the main thread. I think reducing re-allocation or heavy objects will help significantly with performance.

[mborland]

Redefinition of void_t?

[JacobHass8]

It is a redefinition but it was only defined in C++17. I duplicated it because I need it for C++14 test. Is there a better way to do this?

[mborland]

No I meant it as you have a definition here at line 48-49, but there's also another one at 21-22.

[mborland]

This is a copy of a copy, but p0 is never used. Same for q0.

[JacobHass8]

I think p0 is just copied once now when integrate_hamiltonian_imp is called. Should I pass everything by reference when I can? Most function signatures have the parameters const RealType dt, Func dHdp, Func dHdq. The timestep dt is easy to pass by reference. Should I also pass the functions dHdp/dq?

[mborland]

Yes, to the extent possible either pass by reference either const or non-const to the extent you can when moving containers around. Like the in-place add I mentioned I think there are other calculations you can do in-place since you're not mutating the users data.

[mborland]

Can these be competed to higher precision like you did above?

[JacobHass8]

I'll have to reread the paper to see how they did it. Even for the above (simpler) case it was quite complicated. I'll get a better idea and then report how difficult it would be.

[mborland]

Also highly recommend going the enum route instead of a vector of strings and then a map.

enum class available_methods : unsigned 
{
    Y6 = 0,
    Y4 = 1, 
    Y2 = 2 
}

And then you update the signature:

template <typename RandomAccessContainer, typename RealType, class Func, class Policy>
std::pair<std::vector<RandomAccessContainer>, std::vector<RandomAccessContainer> > integrate_hamiltonian_imp(const RandomAccessContainer p0,
                                                                                                             const RandomAccessContainer q0,
                                                                                                             const RealType dt,
                                                                                                             const unsigned steps,
                                                                                                             Func dHdp,
                                                                                                             Func dHdq,
                                                                                                             available_methods method,
                                                                                                             const Policy& pol)

[mborland]

The one other change now is to mark your public APIs with BOOST_MATH_EXPORT. Details should not be exported so they'll stay out of the purview of the consumer.