Add the Störmer-Verlet method + symplectic test changes

[packquickly]

Changes:

• Add the Störmer-Verlet method
• Add all_symplectic_solvers to test/helpers.
• Change test_semi_implicit_euler to a parameterised test test_symplectic_solvers which includes Störmer-Verlet
• Add an order test test_symplectic_ode_order similar to test_ode_order using a simple harmonic oscillator as the base problem instead of the linear equation in test_ode_order.

Störmer-Verlet is implemented in the general partitioned form updating p and q with generic vector fields f(p) and g(q) rather than the more common version which assumes g(q) = q. This is consistent with the Diffrax implementation of SemiImplicitEuler.

[patrick-kidger]

Isn't this just semi-implicit Euler written in kick-drift-kick form? (I.e. offset by half a step.) Justification: it looks to me like y1_2 on this step is y0_2 on the next step, so the term_1.vf_prod(...) evaluations happen at the same point twice. (Which also means that this is increasing runtime/compiletime.)

[packquickly]

This is pretty subtle actually. long story short: the half step difference has an impact, and they are indeed different (you can check numerically, Störmer-Verlet is order 2, symplectic Euler is order 1.)

Störmer-Verlet is the composition of the symplectic euler method and it's adjoint (reverse method,) ie. it is both variants of symplectic Euler stacked with step-size $h/2$:

$$ \begin{aligned} q_{n + 1/2} &= q_n + \frac{h}{2} f(p_n) \\ p_{n + 1/2} &= p_n + \frac{h}{2} g(q_{n + 1/2}) \\ p_{n + 1} &= p_{n + 1/2} + \frac{h}{2} f(q_{n + 1/2}) \\ q_{n + 1} &= q_{n + 1/2} + \frac{h}{2} g(p_{n + 1}) \end{aligned} $$

The implementation in non kick-drift-kick form, ie.

$$ \begin{aligned} q_{n + 1/2} &= q_{n - 1/2} + h f(p_n)\\ p_{n + 1} &= p_n + h g(q_{n + 1/2}) \end{aligned} $$

is distinct from symplectic Euler primarily because of the initialization. Looking at the second-order diffeq case ($g(q) = q$) the initial $p_1$ for symplectic Euler is

$$ p_1 = p_0 + h q_0 + h^2 f(p_0)$$

and for Störmer-Verlet it's:

$$ p_1 = p_0 + hq_0 + \frac{h^2}{2} f(p_0).$$

This is pretty much the only difference for these two though, and the non kick-drift-kick (often called the leapfrog-Verlet implementation, ugh) is definitely the better one when we don't care about knowing $q_{n+1}$. However, I assumed we needed to be able to return the tuple $(p_{n+1}, q_{n+1})$ in each call to step, hence the more expensive kick-drift-kick implementation.

If you see a way to switch to the leapfrog-Verlet implementation though by all means I'll do that instead.

[packquickly]

Also, you might find it interesting that Hairer wrote a long article about Störmer-Verlet as a precursor to "Geometric Numerical Integration," where he used the method to demonstrate a bunch of the ideas later expanded on in the book

[adrn]

I was just about to open a PR with a very similar implementation that I'm using in a project, and then stumbled on this! What's left here to resolve the discussion?