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:

• 6th order method from here
• Kahan summation
• Tests
  • Add pendulum
• 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

❌ Patch coverage is 85.54217% with 12 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.38%. Comparing base (6dea961) to head (645aa1a).
⚠️ Report is 7 commits behind head on develop.

Files with missing lines	Patch %	Lines
include/boost/math/quadrature/symplectic.hpp	80.35%	11 Missing ⚠️
test/test_symplectic.cpp	96.29%	1 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #1411      +/-   ##
===========================================
- Coverage    95.39%   95.38%   -0.01%     
===========================================
  Files          826      828       +2     
  Lines        68901    69002     +101     
===========================================
+ Hits         65726    65820      +94     
- Misses        3175     3182       +7     

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

... 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 6dea961...645aa1a. 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]

@JacobHass8 : It appears to me that if we wanted to return both objects of the same type, we would do a std::array<ReturnType,2>. But it also appears that these two objects have different units. Could it work with dimensioned types like mpunitz or boost::units? That would indeed force you into a std::pair.

[JacobHass8]

I don't think the return object necessarily need to be the same type. It is just returning the input types of p0 and q0.

[JacobHass8]

However, there is a constraint on these types. Either we need to be able to add the types of p0 and q0 or we need the function dHdp to take in a type of p0 and return a type of q0 (for dHdq we need q0 -> p0).

This is because of the lines

    ReturnType p = p0;
    ReturnType q = q0;

    q = q + dt / 2 * dHdp(p);
    p = p - dt * dHdq(q);

I'm not familiar with mpunitz or boost::units. I was only thinking about singletons versus array types.

[NAThompson]

I don't think you used boost unit test framework here right? You just used the boost::math floating point comparisons.

[JacobHass8]

I'm not entirely sure what the unit test framework is. I just copied it from test_trapezoidal...

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

@JacobHass8 : I believe that std::valarray has been essentially deprecated in the C++ community. @mborland : Is this correct?

In general, we like to use template typename <RandomAccessContainer>.

[JacobHass8]

I only used std::valarray because it supports vectorized arithmetic. I think Eigen vectors would be okay too.

Maybe this shouldn't be a requirement though? I would need to either convert to an std::valarray or just write for loops. The for loops option isn't bad but I thought it would obfuscate the algorithm.

[NAThompson]

Yeah, I'd personally just leave it generic and do a for loop.

[JacobHass8]

Hmmm, I can't figure out how to allow the function to accept both numeric types and iterators. I might just require the user to always input an iterator. I can't figure out how to use for loops for iterators and none for numeric types.

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