Behavior of integration method

Question

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):

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.

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}$.

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)$.