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 (64a8d75) to head (c63df4a).

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1411      +/-   ##
===========================================
+ Coverage    95.39%   95.41%   +0.01%     
===========================================
  Files          829      831       +2     
  Lines        69186    69422     +236     
===========================================
+ Hits         66003    66241     +238     
+ Misses        3183     3181       -2     

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

... and 1 file 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 64a8d75...c63df4a. 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.

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

[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]

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]

Was the answer here that these are overly complicated to calculate?

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

[JacobHass8]

@mborland I think this is ready to merge if you don't have any further comments. Thanks for the help!

The drone failure looks spurious. Here's the error message

+ git clone -b develop --depth 1 https://github.com/boostorg/boost.git boost-root
Cloning into 'boost-root'...
error: RPC failed; curl 16 Error in the HTTP2 framing layer
fatal: error reading section header 'shallow-info'

[mborland]

Looks much better. Comments are much smaller in size and scope this go around.

[mborland]

Was the answer here that these are overly complicated to calculate?

[mborland]

Adding a default: BOOST_MATH_UNREACHABLE; will help make sure you don't ever miss a method, for example if you add additional options.

[mborland]

We should probably have some kind of assertion that both p0 and q0 have at least steps amount of elements inside of them otherwise we will end up with a SEGFAULT in the stepper below.

[JacobHass8]

Looks much better. Comments are much smaller in size and scope this go around.

All the tests are passing now. I think this is ready for another look! Thanks for all your help with this.