Test des symplektischen Integrators 3. Ordnung gegen 4. Ordnung mit seltsamem Ergebnis

10

In meiner Antwort auf eine Frage zu MSE bezüglich einer 2D-Hamiltonschen Physiksimulation habe ich vorgeschlagen, einen symplektischen Integrator höherer Ordnung zu verwenden .

Dann hielt ich es für eine gute Idee, die Auswirkungen verschiedener Zeitschritte auf die globale Genauigkeit von Methoden mit unterschiedlichen Ordnungen zu demonstrieren, und schrieb und führte ein Python / Pylab-Skript zu diesem Zweck aus. Zum Vergleich habe ich gewählt:

Das Seltsame ist, egal welchen Zeitschritt ich wähle, Ruths Methode 3. Ordnung scheint in meinem Test genauer zu sein als Ruths Methode 4. Ordnung, selbst um eine Größenordnung.

Meine Frage ist daher: Was mache ich hier falsch? Details unten.

Die Methoden entfalten ihre Stärke in Systemen mit trennbaren Hamiltonianern, dh solchen, die als wobei alle Positionskoordinaten umfasst , die konjugierten Impulse umfasst, die Kinetik darstellt Energie und potentielle Energie.

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

In unserem Setup können wir Kräfte und Impulse durch die Massen normalisieren, auf die sie angewendet werden. So werden Kräfte zu Beschleunigungen und Impulse zu Geschwindigkeiten.

Die symplektischen Integratoren kommen mit speziellen (gegebenen, konstanten) Koeffizienten, die ich mit und . Mit diesen Koeffizienten nimmt ein Schritt zum Entwickeln des Systems vom Zeitpunkt zum Zeitpunkt die Form ana1,,anb1,,bntt+δt

  • Für :i=1,,n

    1. Berechnen Sie den Vektor aller Beschleunigungen unter Berücksichtigung des Vektors aller Positionengq
    2. Ändern Sie den Vektor aller Geschwindigkeiten umvbigδt
    3. Ändern Sie den Vektor aller Positionen umqaivδt

Die Weisheit liegt nun in den Koeffizienten. Dies sind

[a1a2b1b2]=[121201](leap2)[a1a2a3b1b2b3]=[2323172434124](ruth3)[a1a2a3a4b1b2b3b4]=1223[12123212321201231](ruth4)

Zum Testen habe ich das 1D-Anfangswertproblem das einen trennbaren Hamilton-Operator hat. Hier werden mit identifiziert .

y+y=0y(0)=1y(0)=0
(y(t),y(t))=(cost,sint)
(q,v)(y,y)

Ich habe das IVP mit den obigen Methoden über mit einer Schrittgröße von mit einer ganzen Zahl die irgendwo zwischen und . Unter Berücksichtigung der Geschwindigkeit von leap2 habe ich für diese Methode verdreifacht . Ich habe dann die resultierenden Kurven im Phasenraum aufgezeichnet und auf gezoomt, wo die Kurven idealerweise wieder bei ankommen sollten .t[0,2π]δt=2πNN10100N(1,0)t=2πN(1,0)t=2π

Hier sind Diagramme und Zooms für und :N=12N=36

N = 12N = 12, gezoomt

N = 36N = 36, gezoomt

Für kommt Sprung2 mit der Schrittgröße näher zu Hause als ruth4 mit der Schrittgröße . Für , ruth4 gewinnt über leap2 . Allerdings kommt ruth3 mit der gleichen Schrittgröße wie ruth4 in allen Einstellungen, die ich bisher getestet habe , viel näher zu Hause an als die beiden anderen.N=122 π2π3N 2π2πNN=36

Hier ist das Python / Pylab-Skript:

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

Ich habe bereits nach einfachen Fehlern gesucht:

  • Kein Wikipedia-Tippfehler. Ich habe insbesondere die Referenzen überprüft ( 1 , 2 , 3 ).
  • Ich habe die richtige Koeffizientenfolge. Wenn Sie mit der Bestellung von Wikipedia vergleichen, beachten Sie, dass die Sequenzierung der Operatoranwendung von rechts nach links funktioniert. Meine Nummerierung stimmt mit Candy / Rozmus überein . Und wenn ich trotzdem eine andere Bestellung versuche, werden die Ergebnisse schlechter.

Mein Verdacht:

  • Falsche Schrittgrößenreihenfolge: Vielleicht hat Ruths Schema 3. Ordnung irgendwie viel kleinere implizite Konstanten, und wenn die Schrittgröße wirklich klein gemacht würde, würde die Methode 4. Ordnung gewinnen? Aber ich habe sogar ausprobiert , und die Methode 3. Ordnung ist immer noch überlegen.N=360
  • Falscher Test: Etwas Besonderes an meinem Test lässt Ruths Methode 3. Ordnung wie eine Methode höherer Ordnung verhalten?
ccorn
quelle
Können Sie numerische Fehlerwerte angeben? Es ist ein wenig schwer von der Handlung zu unterscheiden. Wie skalieren die Fehler mit der Änderung von ? Skalieren sie wie erwartet aus den Bestellungen der Methoden? Normalerweise würde man Fehler gegen in einem Log-Log-Plot zeichnen, um dies zu überprüfen. N.NN
Kirill
@Kirill: Daran arbeiten. Wird bald bearbeitet.
Ccorn
1
Eine Sache, bei der ich misstrauisch bin, ist die Wahl eines linearen rhs: Kürzungsfehler der Methoden hängen normalerweise von einer hohen Ableitung der rhs ab. Wenn also alle hohen Ableitungen von rhs verschwinden, können Sie ein seltsames Konvergenzverhalten beobachten. Es lohnt sich wahrscheinlich, eine ungewöhnlichere Rhs zu probieren.
Kirill

Antworten:

9

Auf Kirills Vorschlag hin führte ich den Test mit aus einer Liste von ungefähr geometrisch ansteigenden Werten durch und berechnete für jedes den Fehler als wobei die Approximation darstellt erhalten durch numerische Integration. Hier ist das Ergebnis in einem Log-Log-Plot:N ϵ ( N ) = ˜ z ( 2 π ) - ˜ z ( 0 ) 2NN ˜ z

ϵ(N)=z~(2π)z~(0)2wherez~(t)=(y~(t),y~(t))
z~

Geben Sie hier die Bildbeschreibung ein

Also hat ruth3 tatsächlich die gleiche Ordnung wie ruth4 für diesen Testfall und impliziert Konstanten von nur der Größe.141100

Interessant. Ich muss weiter nachforschen und vielleicht andere Tests versuchen.

Übrigens, hier sind die Ergänzungen zum Python-Skript für den Fehlerplot:

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()
ccorn
quelle
Nicht relevant für die Frage, aber können Sie bitte Änderungen und Aktualisierungen in die Frage selbst einfügen, anstatt sie als separate Antwort zu veröffentlichen? Diese Wold die Konvention halten , dass die Antworten sollten beantwortet die Frage.
Kirill
1
@ Kirill: Es ist eine Antwort. ruth3 hat hier tatsächlich höhere Ordnung und niedrigere Konstanten. Entdeckt aufgrund Ihres Vorschlags, ein Protokoll-Protokoll-Fehlerdiagramm zu erstellen. Die Frage wird daher beantwortet, und ich werde den Punkt einer Frage nach ihrer Beantwortung entschieden nicht ändern, auch wenn die Antwort von mir verfasst wurde.
Ccorn
Trotzdem würde ich gerne weitere Analysen akzeptieren. (Selbst beantwortete Fragen werden automatisch akzeptiert, aber ich nehme an, dass man das ändern kann.)
ccorn
2
Ich sah es mir ein wenig an und konnte keine überzeugende Erklärung finden. Diese Konvergenz 4. Ordnung von ruth3 hat viel mit den Anfangsbedingungen zu tun: Versuchen Sie, neq0 einzustellen, und versuchen Sie, nicht für eine volle Periode (noch eine halbe Periode) zu integrieren. Eine Sache, die leicht passieren kann, ist, dass der Kürzungsfehler eine "Null-Mittelwert" -Komponente enthält, die abgebrochen wird, wenn Sie eine vollständige Periode integrieren. Ich habe auch versucht , um zu überprüfen, ob dies etwas mit mit wenigen hohen Ableitungen zu tun hat , aber in meinen Tests scheinen Anfangsbedingungen und Periodizität mehr damit zu tun zu haben. V ( q ) = 1 / q + log q V.p00V(q)=1/q+logqV
Kirill
2
Dies ist eine Anzeige der Superkonvergenz. Einfache Testprobleme wie dieses haben in vielen Fällen dieses Problem. Die Verwendung einer linearen Gleichung kann zu diesem Verhalten führen, und oft können sich die ungeraden Terme einer Taylor-Reihe aufheben, wenn dies geschieht. Bei einem nichtlinearen Testproblem ohne analytische Lösung ist dies viel weniger wahrscheinlich.
Chris Rackauckas
2

Das Zeichnen des Fehlers von , über das gesamte Intervall, skaliert mit der Potenz der Schrittgröße, die durch die erwartete Reihenfolge gegeben ist, ergibt die Diagrammeq¨=qq(0)=1,q˙(0)=0

Geben Sie hier die Bildbeschreibung ein

Wie erwartet nähern sich die Graphen für eine zunehmende Anzahl von Teilintervallen zunehmend einer Grenzkurve an, die der führende Fehlerkoeffizient ist. In allen außer einem Diagramm ist diese Konvergenz sichtbar schnell, es gibt fast keine Divergenz. Dies bedeutet, dass selbst bei relativ großen Schrittgrößen der führende Fehlerterm alle anderen Terme dominiert.

Bei der Ruth-Methode 3. Ordnung scheint der führende Koeffizient der Komponente Null zu sein, die sichtbare Grenzkurve liegt nahe oder gleich der horizontalen Achse. Die sichtbaren Grafiken zeigen deutlich die Dominanz des Fehlerterms 4. Ordnung. Die Skalierung für einen Fehler 4. Ordnung ergibt eine Darstellung ähnlich den anderen.p

Wie man sehen kann, ist in allen 3 Fällen der Fehlerkoeffizient der führenden Ordnung in der Komponente bei nach einer vollen Periode Null . Bei der Kombination der Fehler beider Komponenten dominiert somit das Verhalten der Komponente, was fälschlicherweise den Eindruck einer Methode 4. Ordnung im Loglog-Plot erweckt.qt=2πp

Ein maximaler Koeffizient in der Komponente kann um . Das Loglog-Diagramm dort sollte die korrekten globalen Fehlerreihenfolgen widerspiegeln.q3π/2


Dass das Verschwinden des Fehlerterms 3. Grades in Ruth3p ein Artefakt der Einfachheit der linearen Gleichung ist, zeigt das nichtlineare Beispiel , mit den entsprechenden Darstellungenq¨=sin(q)q(0)=1.3, q˙(0)=0

Geben Sie hier die Bildbeschreibung ein

Lutz Lehmann
quelle