What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

Question

In this comment I wrote:

...default SciPy integrator, which I'm assuming only uses symplectic methods.

in which I am refering to SciPy's odeint, which uses either a "non-stiff (Adams) method" or a "stiff (BDF) method". According to the source:

def odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0,
           ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0,
           hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12,
           mxords=5, printmessg=0):
    """
    Integrate a system of ordinary differential equations.

    Solve a system of ordinary differential equations using lsoda from the
    FORTRAN library odepack.

    Solves the initial value problem for stiff or non-stiff systems
    of first order ode-s::
        dy/dt = func(y, t0, ...)
    where y can be a vector.
    """

Here is an example where I propagate a satellite's orbit around the earth for three months just to show that it precesses as expected.

I believe that non-symplectic integrators have the undesirable property that they will tend not to conserve energy (or other quantities) and so are undesirable in orbital mechanics for example. But I'm not exactly sure what it is that makes a symplectic integrator symplectic.

Is it possible to explain what the property is (that makes a symplectic integrator symplectic) in a straightforward and (fairly) easy to understand but not inaccurate way? I'm asking from the point of view of how the integrator functions internally, rather than how it performs in testing.

And is my suspicion correct that odeint does use only symplectic integrators?

Answer 1

Let me start off with corrections. No, odeint doesn't have any symplectic integrators. No, symplectic integration doesn't mean conservation of energy.

What does symplectic mean and when should you use it?

First of all, what does symplectic mean? Symplectic means that the solution exists on a symplectic manifold. A symplectic manifold is a solution set which is defined by a 2-form. The details of symplectic manifolds probably sound like mathematical nonsense, so instead the gist of it is there is a direct relation between two sets of variables on such a manifold. The reason why this is important for physics is because Hamiltonian's equations naturally have that the solutions reside on a symplectic manifold in phase space, with the natural splitting being the position and momentum components. For the true Hamiltonian solution, that phase space path is constant energy.

A symplectic integrator is an integrator whose solution resides on a symplectic manifold. Because of discretization error, when it is solving a Hamiltonian system it doesn't get exactly the correct trajectory on the manifold. Instead, that trajectory itself is perturbed $\mathcal{O}(\Delta t^n)$ for the order $n$ from the true trajectory. Then there's a linear drift due to numerical error of this trajectory over time. Normal integrators tend to have a quadratic (or more) drift, and do not have any good global guarantees about this phase space path (just local).

What this tends to mean is that symplectic integrators tend to capture the long-time patterns better than normal integrators because of this lack of drift and this almost guarantee of periodicity. This notebook displays those properties well on the Kepler problem. The first image shows what I'm talking about with the periodic nature of the solution.

This was solved using the 6th order symplectic integrator from Kahan and Li from DifferentialEquations.jl. You can see that the energy isn't exactly conserved, but its variation is dependent on how far the perturbed solution manifold is from the true manifold. But since the numerical solution itself resides on a symplectic manifold, it tends to be almost exactly periodic (with some linear numerical drift that you can see), making it do very nicely for long term integration. If you do the same with RK4, you can get disaster:

You can see that the issue is that there's no true periodicity in the numerical solution and therefore overtime it tends to drift.

This highlights the true reason to choose symplectic integrators: symplectic integrators are good on long-time integrations on problems that have the symplectic property (Hamiltonian systems). So let's walk through a few things. Note that you don't always need symplectic integrators even on a symplectic problem. For this case, an adaptive 5th order Runge-Kutta method can do fine. Here's Tsit5:

Notice two things. One, it gets a good enough accuracy that you cannot see the actual drift in the phase space plot. However, on the right side you can see that there is this energy drift, and so if you are doing a long enough integration this method will not do as well as the solution method with the periodic properties. But that raises the question, how does it fare efficiency-wise vs just integrating extremely accurately? Well, this is a bit less certain. In SciMLBenchmarks.jl you can find some benchmarks investigating this question. For example, this notebook looks at the energy error vs runtime on a Hamiltonian equation system from a quadruple Boson model and shows that if you want really high accuracy, then even for quite long integration times it's more efficient to just use a high order RK or Runge-Kutta Nystrom (RKN) method. This makes sense because to satisfy the symplectic property the integrators give up some efficiency and pretty much have to be fixed time step (there is some research making headway into the latter but it's not very far along).

In addition, notice from both of these notebooks that you can also just take a standard method and project it back to the solution manifold each step (or every few steps). This is what the examples using the DifferentialEquations.jl ManifoldProjection callback are doing. You see that guarantees conservation laws are upheld but with an added cost of solving an implicit system each step. You can also use a fully-implicit ODE solver or singular mass matrices to add on conservation equations, but the end result is that these methods are more computationally-costly as a tradeoff.

So to summarize, the class of problems where you want to reach for a symplectic integrator are those that have a solution on a symplectic manifold (Hamiltonian systems) where you don't want to invest the computational resources to have a very exact (tolerance <1e-12) solution and don't need exact energy/etc. conservation. This highlights that it's all about long-term integration properties, so you shouldn't just flock to them all willy-nilly like some of the literature suggests. But they are still a very important tool in many fields like Astrophysics where you do have long time integrations that you need to solve sufficiently fast without having absurd accuracy.

Where do I find symplectic integrators? What kind of symplectic integrators exist?

There are generally two classes of symplectic integrators. There are the symplectic Runge-Kutta integrators (which are the ones shown in the above examples) and there are implicit Runge-Kutta methods which have the symplectic property. As @origimbo mentions, the symplectic Runge-Kutta integrators require that you provide them with a partitioned structure so they can handle the position and momentum parts separately. However, counter to the comment, the implicit Runge-Kutta methods are symplectic without requiring this, but instead require solving a nonlinear system. This isn't too bad because if the system is non-stiff this nonlinear system can be solved with functional iteration or Anderson acceleration, but the symplectic RK methods should still probably be preferred for efficiency (it's a general rule that the more information you provide to an integrator, the more efficient it is).

That said, odeint does not have methods from either of these families, so it is not a good choice if you're looking for symplectic integrators. In Fortran, Hairer's site has a small set you can use. Mathematica has a few built in. The GSL ODE solvers have implicit RK Gaussian point integrators which IIRC are symplectic, but that's about the only reason to use the GSL methods.

But the most comprehensive set of symplectic integrators can be found in DifferentialEquations.jl in Julia (recall this was used for the notebooks above). The list of available symplectic Runge-Kutta methods is found on this page and you'll notice that the implicit midpoint method is also symplectic (the implicit Runge-Kutta Trapezoid method is considered "almost symplectic" because it's reversible). Not only does it have the largest set of methods, but it's also open-source (you can see the code and its tests in a high-level language) and has a lot of benchmarks. A good introductory notebook for using it to solve physical problems is this tutorial notebook. But of course it's recommended you get started with the package through the first ODE tutorial.

In general you can find a detailed analysis of numerical differential equation suites at this blog post. It's quite detailed but since it has to cover a lot of topics it does each at less detail than this, so feel free to ask for it to be expanded in any way.

Answer 2

To complement Chris Rackauckas answer, to state some of the mathematical nonsense as well as some stuff you almost certainly know, a dynamical system is Hamiltonian if there is a description with coordinates $\mathbf{p}$ and $\mathbf{q}$ and a functional, $\mathcal{H(\mathbf{p},\mathbf{q})}$ such that $$\frac{d\mathbf{q}}{dt}=+\frac{\partial \mathcal{H}}{\partial\mathbf{p}}$$ and $$\frac{d\mathbf{p}}{dt}=-\frac{\partial \mathcal{H}}{\partial\mathbf{q}}.$$ This motion conserves the value of $\mathcal{H}$ along trajectories, but it also has an additional property, namely that if we define a mapping $$ \mathbf{p}(t),\mathbf{q}(t)= \phi_t(\mathbf{p}(t_0),\mathbf{q}(t_0))$$ then this mapping conserves the two form $d\mathbf{p}\wedge d\mathbf{q}$. For a problem in which $p$ and $q$ are one dimensional you can think of this as saying that the area inside closed curves on the phase space are conserved. This ensures all kinds of nice stability properties, since "balls" of trajectories have to stay "close" to each other.

In terms of numerics, a symplectic integrator acts in the same way, also conserving this area/two form. In turn this means that there is a conserved "numerical Hamiltonian" (which may not be [read 'is not'] the same as the exact one). Note that stability is not the same as accuracy, so that most of the advantages of symplectic methods come when integrating for very long times (e.g. your method may rapidly place a satellite on the wrong side of the Earth, while never allowing it to decay into it).