Test integratora symetrycznego trzeciego rzędu vs czwartego rzędu z dziwnym wynikiem

W mojej odpowiedzi na pytanie dotyczące MSE dotyczące symulacji fizyki Hamiltoniana 2D zasugerowałem użycie integratora symplektycznego wyższego rzędu .

Potem pomyślałem, że dobrym pomysłem może być wykazanie wpływu różnych kroków czasowych na globalną dokładność metod z różnymi zamówieniami. Napisałem i uruchomiłem skrypt Python / Pylab w tym celu. Do porównania wybrałem:

• ( leap2 ) Przykład drugiego rzędu Wikipedii, który znam, chociaż znam go pod nazwą leapfrog ,
• ( ruth3 ) Integrator symplektyczny trzeciego rzędu Rutha ,
• ( ruth4 ) Integrator symplektyczny czwartego rzędu Rutha .

Dziwną rzeczą jest to, że niezależnie od wybranego przeze mnie przedziału czasu metoda trzeciego rzędu Rutha wydaje się być bardziej dokładna w moim teście niż metoda czwartego rzędu Rutha, nawet o rząd wielkości.

Moje pytanie brzmi zatem: co tutaj robię źle? Szczegóły poniżej.

Metody rozwijają swoją siłę w układach z separowalnymi hamiltonianami, tj. W tych, które można zapisać jako gdzie obejmuje wszystkie współrzędne pozycji, zawiera sprzężoną pęd, reprezentuje kinetykę energia i energia potencjalna

$$H(q,p) = T(p) + V(q)$$ $q$$p$$T$$V$

W naszym ustawieniu możemy znormalizować siły i momenty według mas, do których są przyłożone. W ten sposób siły zamieniają się w przyspieszenia, a momenty w prędkości.

Integratory symplektyczne mają specjalne (podane, stałe) współczynniki, które i . Przy tych współczynnikach przybiera postać jeden etap ewolucji systemu od czasu do czasu$a_1,\ldots,a_n$$b_1,\ldots,b_n$$t$$t+\delta t$

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

  1. Oblicz wektor $g$ wszystkich przyspieszeń, biorąc pod uwagę wektor $q$ wszystkich pozycji
  2. Zmień wektor $v$ wszystkich prędkości o $b_i\,g\,\delta t$
  3. Zmiana wektora $q$ wszystkich pozycjach przez $a_i\,v\,\delta t$

Mądrość leży teraz w współczynnikach. Są to

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

Do testowania wybrałem problem wartości początkowej 1D

$$\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')$

Zintegrowałem IVP z powyższymi metodami w stosunku do ze skokiem wielkości z liczbą całkowitą wybraną gdzieś pomiędzy a . Biorąc pod uwagę szybkość skoku 2 , potroiłem dla tej metody. Następnie wykreśliłem otrzymane krzywe w przestrzeni fazowej i powiększyłem w gdzie krzywe powinny idealnie dotrzeć ponownie przy .$t\in[0,2\pi]$$\delta t=\frac{2\pi}{N}$$N$$10$$100$N(1,0)t=2π$N$$(1,0)$$t=2\pi$

Oto wykresy i powiększenia dla i :$N=12$$N=36$

Dla , skok 2 z krokiem wielkości zdarza się zbliżać do domu niż ruth4 z krokiem wielkości . Dla , ruth4 wygranych ponad leap2 . Jednak ruth3 , z tym samym krokiem co ruth4 , zbliża się znacznie do domu niż oba pozostałe, we wszystkich testowanych do tej pory ustawieniach.$N=12$2 π$\frac{2\pi}{3N}$ 2π$\frac{2\pi}{N}$$N=36$

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

Sprawdziłem już proste błędy:

• Brak literówki w Wikipedii. Sprawdziłem referencje, w szczególności ( 1 , 2 , 3 ).
• Mam prawidłową sekwencję współczynników. Jeśli porównasz z zamówieniami Wikipedii, zauważ, że sekwencjonowanie aplikacji operatora działa od prawej do lewej. Moja numeracja zgadza się z Candy / Rozmus . A jeśli spróbuję jeszcze raz zamówić, wyniki będą gorsze.

Moje podejrzenia:

• Zła kolejność stopni: Może schemat Rutha trzeciego rzędu ma w jakiś sposób znacznie mniejsze implikowane stałe, a jeśli rozmiar kroku byłby naprawdę mały, to metoda czwartego rzędu wygrałaby? Ale nawet wypróbowałem , a metoda trzeciego rzędu wciąż jest lepsza.$N=360$
• Błędny test: Coś specjalnego w moim teście pozwala metodzie Rutha trzeciego rzędu zachowywać się jak metoda wyższego rzędu?

Następujące Cyryla „s sugestii uruchomiony test z z listy przybliżeniu geometrycznie rosnących wartości, a dla każdego z obliczony błędu jako gdzie reprezentuje przybliżenie uzyskane przez całkowanie numeryczne. Oto wynik w dzienniku dziennika:$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$

Zatem ruth3 rzeczywiście ma tę samą kolejność co ruth4 dla tego przypadku testowego i implikuje stałe tylko wielkości.$4$$\frac{1}{100}$

Ciekawy. Będę musiał zbadać dalej, być może próbując innych testów.

Nawiasem mówiąc, oto dodatki do skryptu Python dla wykresu błędów:

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

Wykreślenie błędu , w pełnym przedziale, skalowane według potęgi wielkości kroku podanej w oczekiwanym porządku, daje wykresy$\ddot q=-q$$q(0)=1,\dot q(0)=0$

Zgodnie z oczekiwaniami wykresy rosnącej liczby pod-przedziałów coraz bardziej zbliżają się do krzywej granicznej, która jest wiodącym współczynnikiem błędu. We wszystkich wątkach z wyjątkiem jednego zbieżność ta jest widocznie szybka, prawie nie ma rozbieżności. Oznacza to, że nawet przy względnie dużych rozmiarach kroku wiodący termin błędu dominuje nad wszystkimi innymi terminami.

W metodzie Rutego trzeciego rzędu współczynnik wiodący komponentu wydaje się wynosić zero, widoczna krzywa graniczna jest bliska lub równa osi poziomej. Widoczne wykresy wyraźnie pokazują dominację terminu błędu czwartego rzędu. Skalowanie dla błędu czwartego rzędu daje wykres podobny do innych.$p$

Jak widać, we wszystkich 3 przypadkach współczynnik błędu rzędu wiodącego w składniku wynosi zero przy po pełnym okresie. Łącząc błędy obu komponentów dominuje zatem zachowanie komponentu , co fałszywie daje wrażenie metody czwartego rzędu na wykresie dziennika.$q$$t=2\pi$$p$

Maksymalny współczynnik w składniku można znaleźć w okolicach . Wykres dziennika powinien odzwierciedlać prawidłowe globalne zamówienia na błędy.$q$$3\pi/2$

---

To, że zniknięcie terminu błędu 3 stopnia w Ruth3p jest artefaktem prostoty równania liniowego, pokazuje nieliniowy przykład$\ddot q=-\sin(q)$$q(0)=1.3,~\dot q(0)=0$