Prueba de integrador simpléctico de tercer orden versus cuarto orden con resultado extraño

En mi respuesta a una pregunta sobre MSE con respecto a una simulación de física hamiltoniana en 2D, he sugerido usar un integrador simpléctico de orden superior .

Entonces pensé que podría ser una buena idea demostrar los efectos de diferentes pasos de tiempo en la precisión global de los métodos con diferentes órdenes, y escribí y ejecuté un script Python / Pylab para ese efecto. Para comparar, elegí:

• ( leap2 ) El ejemplo de segundo orden de Wikipedia con el que estoy familiarizado, aunque lo conozco bajo el nombre de leapfrog ,
• ( ruth3 ) Integrador simpléctico de tercer orden de Ruth ,
• ( ruth4 ) Integrador simpléctico de cuarto orden de Ruth .

Lo extraño es que, cualquiera que sea el paso de tiempo que elija, el método de tercer orden de Ruth parece ser más preciso en mi prueba que el método de cuarto orden de Ruth, incluso por un orden de magnitud.

Mi pregunta es por lo tanto: ¿qué estoy haciendo mal aquí? Detalles abajo.

Los métodos despliegan su fuerza en sistemas con Hamiltonianos separables , es decir, aquellos que pueden escribirse como

$$H(q,p) = T(p) + V(q)$$ donde $q$ comprende todas las coordenadas de posición, $p$ comprende el momento conjugado, $T$ representa cinética energía y energía potencial $V$

En nuestra configuración, podemos normalizar fuerzas y momentos por las masas a las que se aplican. Así, las fuerzas se convierten en aceleraciones, y los momentos se convierten en velocidades.

Los integradores simplécticos vienen con coeficientes especiales (dados, constantes) que denominaré $a_1,\ldots,a_n$ y $b_1,\ldots,b_n$ . Con esos coeficientes, un paso para evolucionar el sistema de tiempo $t$ a tiempo $t+\delta t$ toma la forma

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

  1. Calcule el vector $g$ de todas las aceleraciones, dado el vector $q$ de todas las posiciones
  2. Cambiar el vector $v$ de todas las velocidades por $b_i\,g\,\delta t$
  3. Cambiar el vector $q$ de todas las posiciones por $a_i\,v\,\delta t$

La sabiduría ahora reside en los coeficientes. Estos son

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

$$\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)$$ $(q,v)$$(y,y')$

He integrado el IVP con los métodos anteriores sobre con un tamaño de paso de con un entero elegido en algún lugar entre y . Teniendo en cuenta la velocidad de leap2 , tripliqué para ese método. Luego tracé las curvas resultantes en el espacio de fase y amplié la imagen en donde las curvas deberían llegar idealmente de nuevo a .$t\in[0,2\pi]$$\delta t=\frac{2\pi}{N}$$N$$10$$100$$N$$(1,0)$$t=2\pi$

Aquí hay gráficos y zooms para y :$N=12$$N=36$

Para , leap2 con tamaño de paso llega más cerca de casa que ruth4 con tamaño de paso . Para , ruth4 gana sobre leap2 . Sin embargo, ruth3 , con el mismo tamaño de paso que ruth4 , llega mucho más cerca de casa que los otros, en todos los entornos que he probado hasta ahora.$N=12$$\frac{2\pi}{3N}$$\frac{2\pi}{N}$$N=36$

Aquí está el script Python / Pylab:

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

Ya he verificado si hay errores simples:

• No hay error tipográfico en Wikipedia. He comprobado las referencias, en particular ( 1 , 2 , 3 ).
• Tengo la secuencia de coeficientes correcta. Si compara con los pedidos de Wikipedia, tenga en cuenta que la secuencia de la aplicación del operador funciona de derecha a izquierda. Mi numeración está de acuerdo con Candy / Rozmus . Y si intento otro pedido, los resultados empeoran.

Mis sospechas:

• Orden de tamaño incorrecto: tal vez el esquema de tercer orden de Ruth tiene constantes implícitas mucho más pequeñas, y si el tamaño del paso se hiciera realmente pequeño, ¿ganaría el método de cuarto orden? Pero incluso he intentado , y el método de tercer orden sigue siendo superior.$N=360$
• Prueba incorrecta: ¿Algo especial sobre mi prueba permite que el método de tercer orden de Ruth se comporte como un método de orden superior?

Siguiendo la sugerencia de Kirill , ejecuté la prueba con de una lista de valores que aumentan aproximadamente geométricamente, y para cada calculé el error como donde representa la aproximación obtenido por integración numérica. Aquí está el resultado en un gráfico log-log:$N$$N$

$$\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)$$ $\tilde z$

Entonces, ruth3 tiene el mismo orden que ruth4 para ese caso de prueba y constantes implícitas de solo la magnitud.$4$$\frac{1}{100}$

Interesante. Tendré que investigar más, quizás probando otras pruebas.

Por cierto, aquí están las adiciones al script de Python para el diagrama de error:

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

Trazar el error de , en todo el intervalo, escalado por la potencia del tamaño de paso dado por el orden esperado, da las gráficas$\ddot q=-q$$q(0)=1,\dot q(0)=0$

Como era de esperar, los gráficos para aumentar el número de subintervalos se acercan cada vez más a una curva límite, que es el coeficiente de error principal. En todas las parcelas menos una, esta convergencia es visiblemente rápida, casi no hay divergencia. Esto significa que incluso para tamaños de paso relativamente grandes, el término de error principal domina todos los demás términos.

En el método de Ruth de tercer orden, el coeficiente principal del componente parece ser cero, la curva límite visible es cercana o igual al eje horizontal. Los gráficos visibles muestran claramente el dominio del término de error de cuarto orden. Escalar para un error de 4to orden da una gráfica similar a las otras.$p$

Como se puede ver, en los 3 casos el coeficiente de error del orden inicial en el componente es cero en después de un período completo. Al combinar los errores de ambos componentes, el comportamiento del componente domina, dando falsamente la impresión de un método de cuarto orden en el diagrama de registro.$q$$t=2\pi$$p$

Se puede encontrar un coeficiente máximo en el componente alrededor de . El diagrama de registro debe reflejar los pedidos de error global correctos.$q$$3\pi/2$

---

Que la desaparición del término de error de tercer grado en Ruth3p es un artefacto de la simplicidad de la ecuación lineal muestra el ejemplo no lineal , con las parcelas correspondientes$\ddot q=-\sin(q)$$q(0)=1.3,~\dot q(0)=0$