8
$\begingroup$

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 demonstrate the effects of different time steps on the global accuracy of methods with different orders, and I wrote and ran a Python/Pylab script to that effect. For comparison I chose:

The strange thing is, whatever timestep I choose, Ruth's 3rd-order method seems 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 momenta by the masses they are applied to. Thus forces turn into accelerations, and momenta turn into velocities.

The symplectic integrators come with special (given, constant) coefficients which 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 form

  • For $i=1,\ldots,n$:

    1. Compute vector $g$ of all accelerations, given vector $q$ of all positions
    2. Change vector $v$ of all velocities by $b_i\,g\,\delta t$
    3. 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$ chosen somewhere between $10$ and $100$. Taking leap2´s 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$:

N=12N=12, zoomed

N=36N=36, zoomed

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 np
import matplotlib.pyplot as plt

def 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,t

def 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,v

c = 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: -q
qvt0 = (1.0, 0.0, 0.0)
tmax = 2.0 * np.math.pi
N = 36

fig, ax = plt.subplots(1, figsize=(6, 6))
ax.axis([-1.3, 1.3, -1.3, 1.3])
ax.set_aspect('equal')
ax.set_title(r"Phase 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 that sequencing 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 smaller implied 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 method behave like a higher-order method?
$\endgroup$
  • $\begingroup$ Can you give numerical values of errors? It's a little hard to tell from the plot. How do the errors scale with changing $N$? Do they scale as expected from the methods' orders? Usually one would plot errors against $N$ on a log-log plot to check this. $\endgroup$ – Kirill Aug 22 '15 at 19:48
  • $\begingroup$ @Kirill: Working on that. Will edit soon. $\endgroup$ – ccorn Aug 22 '15 at 19:58
  • $\begingroup$ One thing I'm suspicious about is the choice of a linear r.h.s.: truncation errors of the methods usually depend on some high derivative of the r.h.s., so if all high derivatives of r.h.s. vanish, you might observe some weird convergence behaviour. It's probably worth trying a more unusual r.h.s.. $\endgroup$ – Kirill Aug 22 '15 at 20:09
8
$\begingroup$

Following Kirill's suggestion, I ran the test with $N$ from a list of roughly geometrically increasing values, and for each $N$ computed the error as $$\epsilon(N) = \|\tilde z(2\pi) - \tilde z(0)\|_2 \quad\text{where}\quad \tilde z(t) = \left(\tilde{y}(t),\tilde y'(t)\right)$$ where $\tilde z$ represents the approximation obtained by numerical integration. Here is the result in a log-log plot:

enter image description here

So ruth3 indeed has the same order $4$ as ruth4 for that test case and implied constants of only $\frac{1}{100}$ the magnitude.

Interesting. I will have to investigate further, perhaps trying other tests.

By the way, here are the additions to the Python script for the error plot:

def int_error(qvt0, accel, qvt1, Ns, coeffs):
    e = np.empty((len(Ns),))
    for i,N in enumerate(Ns):
        t = np.linspace(qvt0[2], qvt1[2], N+1)
        q,v = symplectic_integrate(qvt0, accel, t, coeffs)
        e[i] = np.math.sqrt((q[-1]-qvt1[0])**2+(v[-1]-qvt1[1])**2)
    return e

qvt1 = (1.0, 0.0, tmax)
Ns = [12,16,20,24,32,40,48,64,80,96,128,160,192,
      256,320,384,512,640,768,1024,1280,1536,2048,2560,3072]

fig, ax = plt.subplots(1)
ax.set_xscale('log')
ax.set_xlabel(r"$N$")
ax.set_yscale('log')
ax.set_ylabel(r"$\|z(2\pi)-z(0)\|$")
ax.set_title(r"Error after 1 period vs #steps")
ax.grid(True)
e = int_error(qvt0, accel, qvt1, Ns, leap2)
ax.plot(Ns, e, label='leap2', color='black')
e = int_error(qvt0, accel, qvt1, Ns, ruth3)
ax.plot(Ns, e, label='ruth3', color='red')
e = int_error(qvt0, accel, qvt1, Ns, ruth4)
ax.plot(Ns, e, label='ruth4', color='blue')
ax.legend(loc='upper right')
fig.show()
$\endgroup$
  • $\begingroup$ Not relevant to the question, but can you please put changes and updates into the question itself, instead of posting as a separate answer? This wold maintain the convention that answers should answer the question. $\endgroup$ – Kirill Aug 22 '15 at 21:46
  • 1
    $\begingroup$ @Kirill: It is an answer. ruth3 indeed has higher order and lesser constants here. Discovered because of your suggestion of making a log-log error plot. The question is therefore answered, and I will decidedly not change the point of a question once it has been answered, not even if the answer has been composed by me. $\endgroup$ – ccorn Aug 22 '15 at 21:53
  • $\begingroup$ That said, I'd be happy to accept further analysis. (Self-answered questions get accepted automatically, but one can change that I suppose.) $\endgroup$ – ccorn Aug 22 '15 at 21:57
  • $\begingroup$ Ah, okay, I misunderstood you. I thought the point of the question was to explain why the order of convergence was wrong, not just demonstrate that the code is working. $\endgroup$ – Kirill Aug 22 '15 at 22:02
  • $\begingroup$ @Kirill: I did not know whether the order was better than expected. (I thought that improbable enough to not consider that possibility seriously.) Thanks to your suggestion, that is clear now. If furthermore it can be explained why that is the case, that would be surely relevant and constitute a very useful answer, $\endgroup$ – ccorn Aug 22 '15 at 22:10

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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