In my answer to a question on MSE regarding a 2D Hamiltonian physics simulation, I have suggested using a higher-order symplectic integrator.Then I thought it might be a good idea to demonstratethe effects of different time steps on the global accuracy of methodswith different orders,and I wrote and ran a Python/Pylab script to that effect.For comparison I chose:(leap2) Wikipedia's 2nd-order examplewith which I am familiar, although I know it under the name leapfrog,(ruth3) Ruth's 3rd-order symplectic integrator,(ruth4) Ruth's 4th-order symplectic integrator.The strange thing is, whatever timestep I choose, Ruth's 3rd-order methodseems to be be more accurate in my test than Ruth's 4th-order method,even by an order of magnitude.My question is therefore: What am I doing wrong here? Details below.The methods unfold their strength in systems with separable Hamiltonians,i.e. those that can be written as$$H(q,p) = T(p) + V(q)$$where $q$ comprises all position coordinates,$p$ comprises the conjugate momenta,$T$ represents kinetic energy and $V$ potential energy.In our setup, we can normalize forces and momentaby the masses they are applied to. Thus forces turn into accelerations,and momenta turn into velocities.The symplectic integrators come with special (given, constant) coefficientswhich I will label $a_1,\ldots,a_n$ and $b_1,\ldots,b_n$.With those coefficients, one step for evolving the system from time $t$to time $t+\delta t$ takes the formFor $i=1,\ldots,n$:Compute vector $g$ of all accelerations, given vector $q$ of all positionsChange vector $v$ of all velocities by $b_i\,g\,\delta t$Change vector $q$ of all positions by $a_i\,v\,\delta t$The wisdom now lies in the coefficients. These are$$\begin{align}\begin{bmatrix}a_1 & a_2 \\b_1 & b_2\end{bmatrix}&= \begin{bmatrix}\frac{1}{2} & \frac{1}{2} \\ 0 & 1\end{bmatrix}&&\textsf{(leap2)} \\\begin{bmatrix}a_1 & a_2 & a_3 \\b_1 & b_2 & b_3\end{bmatrix}&= \begin{bmatrix}\frac{2}{3} & -\frac{2}{3} & 1 \\\frac{7}{24} & \frac{3}{4} & -\frac{1}{24}\end{bmatrix}&&\textsf{(ruth3)} \\\begin{bmatrix}a_1 & a_2 & a_3 & a_4 \\b_1 & b_2 & b_3 & b_4\end{bmatrix}&= \frac{1}{2-\sqrt[3]{2}}\begin{bmatrix}\frac{1}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1-\sqrt[3]{2}}{2} & \frac{1}{2} \\0 & 1 & -\sqrt[3]{2} & 1\end{bmatrix}&&\textsf{(ruth4)}\end{align}$$For testing, I have chosen the 1D initial value problem$$\begin{align}y''+y &= 0 & y(0) &= 1 & y'(0) &= 0\end{align}$$$$\therefore\quad \left(y(t),y'(t)\right) = (\cos t, -\sin t)$$which has a separable Hamiltonian.Here, $(q,v)$ are identified with $(y,y')$.I have integrated the IVP with the above methods over $t\in[0,2\pi]$with a stepsize of $\delta t=\frac{2\pi}{N}$ with an integer $N$ chosensomewhere between $10$ and $100$.Taking leap2s speed into account, I tripled $N$ for that method.I then plotted the resulting curves in phase space and zoomed in at $(1,0)$where the curves should ideally arrive again at $t=2\pi$.Here are plots and zooms for $N=12$ and $N=36$:For $N=12$, leap2 with step size $\frac{2\pi}{3N}$happens to arrive closer to home than ruth4  with step size $\frac{2\pi}{N}$.For $N=36$, ruth4 wins over leap2.However, ruth3, with the same step size as ruth4,arrives much closer to home than both the others,in all settings I have tested so far.Here is the Python/Pylab script:import numpy as npimport matplotlib.pyplot as pltdef symplectic_integrate_step(qvt0, accel, dt, coeffs):    q,v,t = qvt0    for ai,bi in coeffs.T:        v += bi * accel(q,v,t) * dt        q += ai * v * dt        t += ai * dt    return q,v,tdef symplectic_integrate(qvt0, accel, t, coeffs):    q = np.empty_like(t)    v = np.empty_like(t)    qvt = qvt0    q[0] = qvt[0]    v[0] = qvt[1]    for i in xrange(1, len(t)):        qvt = symplectic_integrate_step(qvt, accel, t[i]-t[i-1], coeffs)        q[i] = qvt[0]        v[i] = qvt[1]    return q,vc = np.math.pow(2.0, 1.0/3.0)ruth4 = np.array([[0.5, 0.5*(1.0-c), 0.5*(1.0-c), 0.5],                  [0.0,         1.0,          -c, 1.0]]) / (2.0 - c)ruth3 = np.array([[2.0/3.0, -2.0/3.0, 1.0], [7.0/24.0, 0.75, -1.0/24.0]])leap2 = np.array([[0.5, 0.5], [0.0, 1.0]])accel = lambda q,v,t: -qqvt0 = (1.0, 0.0, 0.0)tmax = 2.0 * np.math.piN = 36fig, ax = plt.subplots(1, figsize=(6, 6))ax.axis([-1.3, 1.3, -1.3, 1.3])ax.set_aspect('equal')ax.set_title(rPhase plot $(y(t),y'(t))$)ax.grid(True)t = np.linspace(0.0, tmax, 3*N+1)q,v = symplectic_integrate(qvt0, accel, t, leap2)ax.plot(q, v, label='leap2 (%d steps)' % (3*N), color='black')t = np.linspace(0.0, tmax, N+1)q,v = symplectic_integrate(qvt0, accel, t, ruth3)ax.plot(q, v, label='ruth3 (%d steps)' % N, color='red')q,v = symplectic_integrate(qvt0, accel, t, ruth4)ax.plot(q, v, label='ruth4 (%d steps)' % N, color='blue')ax.legend(loc='center')fig.show()I have checked for simple errors already:No Wikipedia typo. I have checked the references, in particular(1,2,3).I have got the coefficient sequence right.If you compare with Wikipedia's ordering, note thatsequencing of operator application works from right to left.My numbering agrees with Candy/Rozmus.And if I try another ordering nevertheless, results get worse.My suspicions:Wrong stepsize order: Maybe Ruth's 3rd-order scheme has somehow much smallerimplied constants, and if the step size were made really small,then the 4th-order method would win?But I have even tried $N=360$, and the 3rd-order method still is superior.Wrong test: Something special about my test lets Ruth's 3rd-order methodbehave like a higher-order method?
Test of 3rd-order vs 4th-order symplectic integrator with strange result
python;time integration