Störmer-Verlet is an order 2 symplectic method. This means it will preserve the Hamiltonian for the system the differential equations describe, for example, important for modelling planetary motion; the application of something like the much-loved Runge-Kutta 4th order method would either model the planet spiralling toward or away from the Sun!

Here's a diagram showing the orbit of Jupiter around the Sun.

To create this, consider the \(n\)-body problem. The Hamiltonian is

\[ {\mathbb H} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} - \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|} \]

Apply Hamilton's equations will gives \(2n\) first order equations. To use stormerVerlet2 this needs to be \(n\) second order equations. In this case, the Lagrangian is easy

\[ {\mathcal{L}} = \frac{1}{2}\sum_{i=0}^n \frac{p_i^\top p_i}{m_i} + \frac{G}{2}\sum_{i=0}^n\sum_{j \neq i} \frac{m_i m_j}{\|q_i - q_j\|} \]

Applying Lagrange's equation

\[ \frac{\mathrm{d}}{\mathrm{d}t}\bigg(\frac{\partial{\mathcal{L}}}{\partial\dot{q}_j}\bigg) = \frac{\partial{\mathcal{L}}}{\partial{q}_j} \]

gives

\[ m_j\ddot{q}_j = G\sum_{k \neq j}m_k m_j \frac{q_k - q_j}{\|q_k - q_j\|^3} \]

For \(n = 2\) this gives

\[ \begin{aligned} \ddot{q}_1 &= m_2G\frac{q_1 - q_2}{\|q_1 - q_2\|^3} \\ \ddot{q}_2 &= m_1G\frac{q_2 - q_1}{\|q_2 - q_1\|^3} \end{aligned} \]

{-# LANGUAGE NegativeLiterals      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}

import qualified Data.Vector as V
import Control.Monad.ST

import Math.Integrators.StormerVerlet
import Math.Integrators

import qualified Linear as L
import Linear.V
import Data.Maybe ( fromJust )

import Diagrams.Prelude

import Control.Monad
import Control.Monad.State.Class

import Plots

-- First some constants describing the system

gConst :: Double
gConst = 6.67384e-11

nStepsTwoPlanets :: Int
nStepsTwoPlanets = 44

-- A step size of 100 days!

stepTwoPlanets :: Double
stepTwoPlanets = 24 * 60 * 60 * 100

sunMass, jupiterMass :: Double
sunMass     = 1.9889e30
jupiterMass = 1.8986e27

jupiterPerihelion :: Double
jupiterPerihelion = 7.405736e11

jupiterV :: L.V3 Double
jupiterV = L.V3 (-1.0965244901087316e02) (-1.3710001990210707e04) 0.0

jupiterQ :: L.V3 Double
jupiterQ = L.V3 (-jupiterPerihelion) 0.0 0.0

sunV :: L.V3 Double
sunV = L.V3 0.0 0.0 0.0

sunQ :: L.V3 Double
sunQ =  L.V3 0.0 0.0 0.0

-- The right hand side of the second order differential equation system.

kepler :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double)
kepler (L.V2 q1 q2) =
    let r  = q2 L.^-^ q1
        ri = r `L.dot` r
        rr = ri * (sqrt ri)
        q1' = pure gConst * r / pure rr
        q2' = negate q1'
        q1'' = q1' * pure sunMass
        q2'' = q2' * pure jupiterMass
    in L.V2 q1'' q2''

-- Initial values

initPQs :: L.V2 (L.V2 (L.V3 Double))
initPQs = L.V2 (L.V2 jupiterV sunV) (L.V2 jupiterQ sunQ)

-- Steps at which to evolve the system

tm :: V.Vector Double
tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets

-- The results

result1 :: V.Vector (L.V2 (L.V2 (L.V3 Double)))
result1 = runST $ integrateV (\h -> stormerVerlet2 kepler (pure h)) initPQs tm

preMorePts :: [(Double, Double)]
preMorePts = map (\(L.V2 _ (L.V2 (L.V3 x y _z) _)) -> (x,y))  (V.toList result1)

morePts :: [P2 Double]
morePts = map p2 $ preMorePts

-- Finally plot the results

addPoint :: (Plotable (Diagram B) b, MonadState (Axis b V2 Double) m) =>
            Double -> (Double, Double) -> m ()
addPoint o (x, y) = addPlotable'
                    ((circle 1e11 :: Diagram B) #
                     fc brown #
                     opacity o #
                     translate (r2 (x, y)))

jSaxis :: Axis B V2 Double
jSaxis = r2Axis &~ do
  addPlotable' ((circle 1e11 :: Diagram B) # fc yellow)
  let l = length preMorePts
  let os = [0.05,0.1..]
  let ps = take (l `div` 4) [0,4..]
  zipWithM_ addPoint os (map (preMorePts!!) ps)
  linePlot' $ map unp2 $ take 200 morePts

jupiterOrbit = renderAxis jSaxis # bg ivory