Symplectic Integrators

Dynamical systems can be divided into two basic types: conservative and dissipative. In biology, we almost always model dissipative systems and thus if we want to computationally simulate the system almost any numerical solver will do the job (unless the problem is stiff, which I’ll leave to another post). However, when simulating a conservative system, we must take care to conserve the conserved quantities. Here, I will give a very elementary review of symplectic integrators for numerically solving conservative systems.

The prime example is a Hamiltonian system, where energy is conserved. More precisely, the symplectic form, which includes orientation as well as magnitude, is preserved. Consider the simple harmonic oscillator with the Hamiltonian:

$H = \frac{1}{2} p^2 + \frac{1}{2}q^2$

Hamilton’s equations are

$\dot{q}=\frac{\partial H}{\partial p}=p$

$\dot{p}=-\frac{\partial H}{\partial q}=-q$

If we were to blindly use Euler’s method to numerically integrate these equations we would get

$q^{t+1}=q^{t}+ h p^t$

$p^{t+1}=p^{t}-h q^t$

where h is the step size. We can rewrite this as

$\left(\begin{array}{c} q^{t+1} \\ p^{t+1}\end{array}\right)=\left(\begin{array}{cc} 1 & h\\ -h & 1\end{array}\right)\left(\begin{array}{c} q^{t} \\ p^{t}\end{array}\right)$

We assess stability of this numerical scheme by setting

$(q^t,p^t) =\lambda^t (v,w)$

, in other words, compute the eigenvalues of the matrix system. The characteristic equation,

$(\lambda-1)^2 +h^2$

, has solutions

$\lambda = 1\pm i h$

. The magnitude of the largest eigenvalue is

$\sqrt{1+h^2}$

, which is greater than one so the system is unstable. Hence, the Euler scheme will never work for a Hamiltonian system no matter how small you make h. Eventually, it will blow up.

This can be fixed with an extremely simple tweak to preserve the symplectic form and naturally enough this is called a symplectic integrator:

$q^{t+1}=q^{t}+ h p^t$

$p^{t+1}=p^{t}-h q^{t+1}$

The only difference is that you update the p equation using the updated value of q. In fact, if you were to absentmindedly code this system you could easily accidentally discover the symplectic integrator. However, you wouldn’t know why it worked or even that it may not always work and that is a very dangerous thing.

If we substitute

$(q^t,p^t) =\lambda^t (v,w)$

, we now get the system

$\lambda v = v-hw$

$\lambda(w-hv)=w$

which has solutions

$\lambda = (2-h^2)\pm h \sqrt{h^2-4})/2$

, which are always magnitude one. Hence, the symplectic integrator is stable and the energy will remain bounded. The symplectic scheme generalizes to nonlinear Hamiltonian equations as well. The trick is just that you update each equation sequentially using all the information up to that point in time.

This idea can generalize to PDEs as well. Consider the advection equation

$\partial_t v -\partial_x v = 0$