12
$\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?
Improve this question
$\endgroup$
3
  • $\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, 2015 at 19:48
  • $\begingroup$ @Kirill: Working on that. Will edit soon. $\endgroup$
    – ccorn
    Aug 22, 2015 at 19:58
  • 1
    $\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, 2015 at 20:09

2 Answers 2

Reset to default
9
$\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()
Improve this answer
$\endgroup$
9
  • $\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, 2015 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, 2015 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, 2015 at 21:57
  • 2
    $\begingroup$ I looked at it a little, and I couldn't find a convincing explanation. This 4th-order convergence of ruth3 has a lot to do with initial conditions: try setting $p_0\neq0$ and try not integrating for a full period (nor half-period). One thing that can happen easily is that the truncation error has some "zero-mean" component that cancels out when you integrate to a full period. I also tried $V(q)=1/q+\log q$ to check if this has anything to do with $V$ having few high derivatives, but in my tests it seems initial conditions and periodicity have more to do with this. $\endgroup$
    – Kirill
    Aug 23, 2015 at 1:25
  • 3
    $\begingroup$ This is a display of superconvergence. Simple test problems like this end up having this issue in many cases. Using a linear equation can give this behavior, and many times the odd terms of a Taylor series can cancel out when that happens. A nonlinear test problem without an analytical solution is much less likely to have this occur. $\endgroup$
    – Chris Rackauckas
    Apr 7, 2019 at 16:21
5
$\begingroup$

Plotting the error of $\ddot q=-q$, $q(0)=1,\dot q(0)=0$ over the full interval, scaled by the power of the step size given by the expected order, gives the plots

enter image description here

As expected, the graphs for increasing number of sub-intervals increasingly approach a limit curve which is the leading error coefficient. In all but one plot this convergence is visibly rapid, there is almost no divergence. This means that even for relatively large step sizes the leading error term dominates all the other terms.

In the 3rd order Ruth method, the leading coefficient of the $p$ component is appears to be zero, the visible limit curve is close or equal to the horizontal axis. The visible graphs clearly show the dominance of the 4th order error term. Scaling for an 4th order error gives a plot similar to the others.

As one can see, in all 3 cases the error coefficient of the leading order in the $q$ component is zero at $t=2\pi$ after a full period. In combining the errors of both components, the behavior of the $p$ component thus dominates, falsely giving the impression of a 4th order method in the loglog plot.

A maximal coefficient in the $q$ component can be found around $3\pi/2$. The loglog plot there should reflect the correct global error orders.

That the vanishing of the 3rd degree error term in Ruth3p is an artefact of the simplicity of the linear equation shows the non-linear example $\ddot q=-\sin(q)$, $q(0)=1.3,~\dot q(0)=0$ with the corresponding plots

enter image description here

Improve this answer
$\endgroup$

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.