I was playing with N-body simulations of a game called Kerbal Space Program, which itself uses the patched conics approximation. I have read that for long term stability it is best to use symplectic integrators and one of the simplest would be Störmer-Verlet method: $$ \ddot{\vec{x}}_n = F(\vec{x}_n) $$ $$ \vec{x}_{n+1} = 2\vec{x}_n - \vec{x}_{n-1} + F(\vec{x}_n){\Delta t}^2 $$ $$ t_{n+1} = t_n + \Delta t $$ I thought that I might add to this method using a similar approach as how this method can be derived, starting with the Taylor expansions: $$ \vec{x}(t+\Delta t) = \sum^N_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^{N+1}) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = \sum^5_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + \sum^5_{n=0} \frac{(-\Delta t)^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^6) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = 2\vec{x}(t) + \Delta t^2 \frac{\partial^2}{\partial t^2}\vec{x}(t) + \frac{\Delta t^4}{12} \frac{\partial^4}{\partial t^2}\vec{x}(t) + O(\Delta t^6) $$ Since the accelerations can be calculated with $F(\vec{x}(t))$ I only had to find an expression for the forth derivative (jounce or snap). For this I used finite difference: $$ \frac{\partial^4}{\partial t^4}\vec{x}(t) = \frac{\ddot{\vec{x}}(t+\Delta t) - 2\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{\Delta t^2} + O(\Delta t^2) $$ $$ \vec{x}(t+\Delta t) = 2\vec{x}(t) - \vec{x}(t-\Delta t) + \Delta t^2 \frac{\ddot{\vec{x}}(t+\Delta t) + 10\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{12} + O(\Delta t^6) $$ This method would be implicit, since for the next position you would need to know the acceleration at that point, however I thought that higher order precision might give it some advantage over Störmer-Verlet. Similar to Störmer-Verlet I assumed that the global error would only drop two orders, such that this method would be a forth order integrator.

Kepler Orbit

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint): Error of integrator versus time-step size

The actual slopes are $4.98\pm 0.01$ for my method and $1.98\pm 0.01$ for Verlet, where the errors are a 95% confidence bound.

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below. enter image description here

The actual slopes are $8.00\pm 0.01$ and $3.998\pm 0.003$ for my method and $3.992\pm 0.005$ and $1.997\pm 0.002$ for Verlet, where the errors are again a 95% confidence bound.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?

PS: In a similar way an eighth order method can be found, by increasing the Taylor expansion by two orders and by using a forth order approximation of the second derivative of the acceleration and a second order approximation of the forth derivative of the acceleration: $$ \vec{x}_{n+1} = 2\vec{x}_{n} - \vec{x}_{n-1} + \Delta t^2 \frac{3\ddot{\vec{x}}_{n+1} + 48\ddot{\vec{x}}_{n+1/2} + 78\ddot{\vec{x}}_{n} + 48\ddot{\vec{x}}_{n-1/2} + 3\ddot{\vec{x}}_{n-1}}{180} + O(\Delta t^8) $$ However I have not yet figured out a good way to find $\ddot{\vec{x}}_{n+1/2}$.

  • $\begingroup$ Typically, you prove that a method is symplectic. Methods for doing that can be found in a "Geometric Numerical Integration" by Hairer, Wanner, and Lubich. (In general, the books by Hairer and Wanner are good reference texts for numerical integration of ODEs.) Numerical results are then usually used for illustration and to provide a sanity check for the proof (i.e., if you obviously fail to conserve energy for long time intervals, your method is probably not symplectic, and your proof is probably wrong). $\endgroup$ Sep 25, 2014 at 1:57
  • $\begingroup$ Symplecticity is the property of conserving area in phase space; a symplectic method would map $\Phi:(x_n,v_n)\mapsto (x_{n+1},v_{n+1})$ and preserve areas: $\mathrm{vol}(\Phi(\Omega))=\mathrm{vol}(\Omega)$ for some set $\Omega=\{(x,v)\}\subset\mathbb{R}^2$. Your method only maps $(x_n,x_{n-1})\mapsto (x_{n+1},x_n)$, without any mention of velocities. I don't know what it would mean for such a method to be symplectic. What is your definition of a symplectic method? $\endgroup$
    – Kirill
    Sep 25, 2014 at 4:12
  • $\begingroup$ Also, when I implemented this method, and tried it on $\ddot x=-x$ for one period, my estimate of its order was $5$, which matches the local truncation error of order $6$. Are you sure it's sixth-order? $\endgroup$
    – Kirill
    Sep 25, 2014 at 4:14
  • $\begingroup$ @Kirill how did you determine the second position and how did you solve the implicitly? Since if you approximate these with an lower order, than the end result would have a lower order too. I only based my estimate of the order on the slope of the error of my figure. I will try and see what results I get in terms of error for this ODE. $\endgroup$
    – fibonatic
    Sep 25, 2014 at 10:59
  • $\begingroup$ @Kirill and isn't (position) Störmer-Verlet method symplectic, which also only uses positions? Also I have to admit that the slope in my figure is a close call between 6 and 5. $\endgroup$
    – fibonatic
    Sep 25, 2014 at 11:24

1 Answer 1


Your method is an example of a multistep linear method for a second-order ODE. Have a look at Hairer, Wanner, Norsett, "Solving Ordinary Differential Equations I", they discuss this method on p. 464, and call it the Implicit Stormer method.

Multistep linear methods for a first-order ODE are those methods that have the form $$ \alpha_0 y_n + \alpha_1 y_{n+1} + \cdots + \alpha_s y_{n+s} = \beta_0 f_n + \beta_1 f_{n+1} + \cdots + \beta_s f_{n+s}, $$ where $\alpha$'s and $\beta$'s are chosen for consistency and numerical stability. In your case, one way to generalize is to use $$ y_{n+1}-2y_{n}+y_{n-1} = \Delta t^2 \big(\beta_{-1}f_{n+1} + \beta_0 f_n + \beta_1 f_{n-1} + \beta_2 f_{n-2} + \cdots), $$ and choose the coefficients appropriately. Hairer, Wanner, Norsett list other related methods as well.

Order. I tested your method on $$\ddot x = x, \qquad x(0) = 1, \quad x(-h) = e^{-h}, \quad h = 1/n $$ with $n=2^{\{3,\ldots,13\}}$, integrating to $t=1$. The error of $x(1)$ has order $4$ as expected.

For a method for a second-order ODE with local error of order $p+2$ the global error has order $p$. (So some of my comments are dead wrong.) The easy way to see this is to note that the inhomogeneous recurrence equation $$ y_{n+1}-2y_n+y_{n-1}=d $$ has the particular solution $y_n=\frac12 dn^2$, which is $O(h^p)$ when $d=O(h^{p+2})$.

Symplecticity. When a system of ODEs comes from a Hamiltonian, i.e., when it is $p'=-H_q$, $q'=H_p$, it is symplectic: the area (or volume) $dp\,dq$ in $\mathbb{R}^2$ is preserved by the flow defined by the equation. Here $p,q$ are generalized momentum and generalized coordinate; the space of all $\{(p,q)\}$ is called the phase space and all this is standard in physics textbooks.

"Geometric Numerical Integration" by Hairer, Wanner, and Lubich, mentioned by Geoff Oxberry in a comment above, is a good reference for symplectic methods. The definition of symplecticity for an ODE method is the same as in physics: a method that defines a map $(x_n,v_n)\mapsto (x_{n+1}, v_{n+1})$ between consecutive timesteps is symplectic if it preserves area in $\mathbb{R}^2$.

As Geoff Oxberry said, to establish symplecticity it is not enough to numerically establish that it conserves energy of a pendulum for a reasonably large number of periods.

As you have written it, I would say your method is not symplectic because it does not define a flow, let alone a symplectic flow. This says nothing about whether it conserves or exponentially decays energy, just that it cannot be said to be symplectic in this form.

Note that the Stormer-Verlet method $$ y_{n+1}-2y_n+y_{n-1} = h^2 f_n $$ does not define a symplectic flow, but it can be rewritten as $$ v_{n+\frac12} = v_n+\frac h2f_n, \quad y_{n+1} = y_n+h v_{n+\frac12}, \quad v_{n+1} = v_{n+\frac12} + \frac h2f_{n+1}, $$ which does define a map $(y_n,v_n)\mapsto (y_{n+1},v_{n+1})$, and it is symplectic because $dy_{n+1}\wedge dv_{n+1} = dy_n\wedge dv_n$ ("$\wedge$" is the anti-symmetric product defined on differential forms; it is the product preserved by symplectic maps).

When I tried to write your method in a form applied to $(y_n, v_n)$, I got (for an arbitrary parameter $\theta$): $$ y_{n+1} = y_n + hv_n + \tfrac1{12}h^2f_{n+1} + h^2(\theta-\tfrac1{12})f_{n}, \\ v_{n+1} = v_n + h\theta f_{n} + h(1-\theta)f_{n+1}, $$ for which I got $$ dy_{n+1}\wedge dv_{n+1} = \frac{1-\frac{h^2}{12}f'_n}{1-\frac{h^2}{12}f'_{n+1}} dy_{n}\wedge dv_n. $$ So unless I made a mistake here, your method is only symplectic when the r.h.s. of $y''=f(y)$ is linear in $y$. When it's not linear, it will scale volume in phase space by a factor $1+O(h^2)$.

  • $\begingroup$ You say that this method is a fifth order, however I also got results (with the same simple harmonic oscillator) where it seemed to be forth order at your equivalent end position of $x(\frac{n}{2}\Delta t)$. Did I do something wrong, or could it actually be only forth order? $\endgroup$
    – fibonatic
    Sep 26, 2014 at 15:11
  • $\begingroup$ @fibonatic I made a mistake there. If the local error has order $p+2$, then the method has order $p$; I forgot it was a 2nd-order ODE, not a 1st-order one. The equation $\ddot x=-x$ is not very good for checking this, I think $\ddot x=x$ works better. $\endgroup$
    – Kirill
    Sep 26, 2014 at 22:19
  • $\begingroup$ I am not familiar with the anti-symmetric product. When subtracting two consecutive the Taylor expansions of the velocity you can obtain a fifth order expression for the velocity: $$ \vec{v}_{n+1} = \vec{v}_{n-1} + \frac{\Delta t}{3}\left(\ddot{\vec{x}}_{n+1} + 4\ddot{\vec{x}}_{n} + \ddot{\vec{x}}_{n-1}\right) $$ Would this conserve the area? $\endgroup$
    – fibonatic
    Oct 8, 2014 at 16:04
  • $\begingroup$ @fibonatic It must be in the form $(y_n,v_n)\mapsto (y_{n+1}, v_{n+1})$; $v_{n-1}$ can be a part of it. $\endgroup$
    – Kirill
    Oct 8, 2014 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.