Die N-Körper-Simulation verliert immer noch an Genauigkeit, obwohl Arithmetik mit beliebiger Genauigkeit und symplektischer Integrator verwendet werden

Es gibt verschiedene mögliche Ursachen dafür, dass sich ein Umlaufbahnsimulator schlecht verhält, selbst wenn ein hervorragendes Paket mit beliebiger Genauigkeit wie mpmath verwendet wird . Wenn Sie numerische Werte bereitstellen, müssen Sie sie als Zeichenfolgen oder Ganzzahlen und nicht als Floats an mpmath übergeben, damit sie ordnungsgemäß konvertiert werden können. Bitte lesen Sie die mpmath-Dokumentation für weitere Details.

Eine weitere mögliche Ursache für Probleme ist die katastrophale Stornierung , die auftreten kann, wenn Zahlen mit sehr unterschiedlichen Größen kombiniert werden. Sie können dieses Problem oft minimieren, indem Sie die Reihenfolge der Berechnungen sorgfältig organisieren und in einer Skala arbeiten, in der der Bereich der Größenordnungen Ihrer Zahlen nicht zu groß ist. Wenn Sie eine Bibliothek wie mpmath verwenden, können Sie natürlich einfach die Genauigkeit erhöhen, obwohl dies die Dinge verlangsamt.

Wenn eine korrekt codierte Orbit-Simulation unannehmbar hohe Fehler liefert, ist es üblich, den Zeitschritt kleiner zu machen. Die Verwendung einer großen Anzahl von Zeitschritten erhöht jedoch die kumulativen Rundungsfehler. Wenn Sie also den Zeitschritt zu klein machen, wird das Ergebnis tatsächlich schlechter, wie David Hammen bereits erwähnt hat. Wenn Sie mpmath verwenden, können Sie dieses Problem natürlich umgehen, indem Sie die Genauigkeit erhöhen.

Als Referenz finden Sie hier einen Python-Code, der eine einfache 1-Körper-Orbit-Simulation unter Verwendung der synchronisierten Form der Leapfrog-Integration durchführt, die eigentlich mit dem Geschwindigkeits-Verlet- Algorithmus identisch ist . Es verwendet einfache Python-Gleitkommaarithmetik, nicht mpmath. Der Python- floatTyp verwendet das IEEE-754-Binary64- Format mit einer Genauigkeit von 53 Bit (fast 16 Dezimalstellen).

Dieses Programm verwendet natürliche Einheiten, dh der Gravitationsparameter des Zentralkörpers ist$GM=\mu=4\pi^2$, also hat eine Kreisbahn mit Radius 1 Raumeinheit eine Periode von 1 Zeiteinheit.

Bei einer einfachen kreisförmigen Umlaufbahn mit 10.000 Zeitschritten pro Umlaufbahn liegt der maximale Fehler im Radius unter 2 Teilen von 10.000.000. Für exzentrische Umlaufbahnen ist ein kleinerer Zeitschritt erforderlich.

''' Simple orbit sim using the
    synchronised Leapfrog algorithm
    Written by PM 2Ring 2017.05.20
'''

from math import pi

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # Time step
        self.delta_time = delta_time

        # Current acceleration due to central gravity
        self.ax, self.ay = acc(x, y)

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the synchronised leapfrog algorithm
        '''
        dt = self.delta_time
        dt2 = dt * dt

        # Update position
        x, y = self.x, self.y
        x += self.vx * dt + 0.5 * self.ax * dt2
        y += self.vy * dt + 0.5 * self.ay * dt2

        # Update velocity using mean acceleration
        ax, ay = acc(x, y)
        self.vx += 0.5 * (self.ax + ax) * dt
        self.vy += 0.5 * (self.ay + ay) * dt

        self.ax, self.ay = ax, ay
        self.x, self.y = x, y
        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 10000
delta_time = period / steps
istep = steps // 100

body = Body(radius,0, 0,v, delta_time)

print("step radius x y")
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

# Plot
P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

Das ist größtenteils einfacher Python-Code, abgesehen von den letzten paar Zeilen, die Sage zum Plotten verwenden (über matplotlib).

Hier ist die Ausgabe. Wie Sie sehen können, tritt der maximale Fehler auf halber Strecke auf.

Semi-major axis: 100.0
Period: 1000.0
step radius x y
100 100.00000001947542 99.8026728363567 6.279052365941193
200 100.00000007782481 99.2114701058467 12.533324180026883
300 100.00000017481786 98.22872501631164 18.738132688008204
400 100.00000031007181 96.85831601489033 24.86899034488533
500 100.0000004830529 95.10565148152806 30.901701456144753
600 100.00000069307836 92.97764838452808 36.812457667191076
700 100.0000009393194 90.48270498238338 42.577931924118325
800 100.00000122080417 87.63066767962191 48.17537053499955
900 100.00000153642179 84.43279216747251 53.582682968369106
1000 100.00000188492673 80.90169900271201 58.7785290345032
1100 100.0000022649435 77.05132380000519 63.74240310543257
1200 100.00000267497244 72.89686223430738 68.45471504130965
1300 100.00000311339528 68.4547100703815 72.89686750374857
1400 100.0000035784818 63.74239845610735 77.05132935101909
1500 100.00000406839649 58.77852473495271 80.90170482543533
1600 100.00000458120591 53.582679050660076 84.43279825988925
1700 100.00000511488614 48.17536703380712 87.63067404816206
1800 100.00000566733111 42.577928875365075 90.48271164233866
1900 100.00000623636048 36.812455106635035 92.97765536027522
2000 100.0000068197286 30.90169941794281 95.10565880655513
2100 100.00000741513318 24.86898886015765 96.85832373162384
2200 100.00000802022444 18.73813178343009 98.22873317574528
2300 100.00000863261437 12.53332387647376 99.21147876697717
2400 100.00000924988619 6.279052677211943 99.8026820654341
2500 100.00000986960366 9.316423505156823e-07 100.00000986960366
2600 100.00001048932121 -6.279050817680861 99.80268342431145
2700 100.00001110659296 -12.533322029108925 99.21148149399528
2800 100.00001171898288 -18.738129959342203 98.22873728916358
2900 100.00001232407413 -24.868987074894257 96.8583292581682
3000 100.00001291947869 -30.901697693138058 95.10566578063481
3100 100.00001350284674 -36.81245347153712 92.97766382296042
3200 100.00001407187602 -42.577927368209636 90.48272164011684
3300 100.0000146243209 -48.17536570300814 87.63068563149452
3400 100.00001515800102 -53.582677955791624 84.43281148159379
3500 100.00001567081028 -58.7785239475003 80.9017197389242
3600 100.00001616072485 -63.74239805996919 77.05134600842264
3700 100.00001662581123 -68.45471016210631 72.89688595395924
3800 100.0000170642339 -72.8968629230579 68.45473532798135
3900 100.00001747426262 -77.0513252072431 63.742425264971644
4000 100.00001785427928 -80.90170126160325 58.77855309408762
4100 100.00001820278406 -84.43279542201651 53.58270894403392
4200 100.00001851840162 -87.6306720835148 48.17539842982273
4300 100.00001879988628 -90.48271069762764 42.577961726543805
4400 100.00001904612719 -92.97765557980723 36.81248934952026
4500 100.00001925615258 -95.10566033036771 30.9017349732278
4600 100.00001942913354 -96.85832669363184 24.869025633039456
4700 100.00001956438729 -98.22873770192699 18.738169664187293
4800 100.00001966138026 -99.21148497362778 12.533362741262588
4900 100.00001971972958 -99.80269005751201 6.279092389071799
5000 100.00001973920497 -100.00001973919642 4.13416988284121e-05
5100 100.00001971972962 -99.80269524924384 -6.279009868830847
5200 100.00001966138022 -99.21149533660193 -12.533280709847839
5300 100.0000195643871 -98.2287531952455 -18.738088445339077
5400 100.00001942913315 -96.85834725614961 -24.868945547291286
5500 100.00001925615217 -95.10568588093415 -30.901656336641377
5600 100.00001904612671 -92.97768601758592 -36.812412472438055
5700 100.00001879988581 -90.48274590249471 -42.577886912364384
5800 100.00001851840106 -87.6307119165327 -48.17532597380332
5900 100.00001820278348 -84.43283972598266 -53.582639132125166
6000 100.00001785427867 -80.90174986167018 -58.778486201805165
6100 100.00001747426198 -77.05137791160872 -63.74236155630852
6200 100.00001706423315 -72.89691952372205 -68.45467505436653
6300 100.00001662581037 -68.4547704356921 -72.89682935326505
6400 100.00001616072397 -63.7424617686046 -77.05129330402596
6500 100.00001567080935 -58.77859083975651 -80.9016711388251
6600 100.00001515800002 -53.582747767675556 -84.4327671775946
6700 100.00001462431985 -48.175438159004266 -87.63064579844279
6800 100.0000140718749 -42.578002182367385 -90.4826864352153
6900 100.00001350284549 -36.81253034859931 -92.97763338514665
7000 100.00001291947744 -30.90177632970622 -95.10564023003296
7100 100.00001232407287 -24.86906716062583 -96.8583086956148
7200 100.00001171898157 -18.738211178175494 -98.22872179580955
7300 100.0000111065917 -12.533404060510405 -99.21147113098588
7400 100.00001048931996 -6.279133337910164 -99.80267823254468
7500 100.00000986960247 -8.175174521653228e-05 -100.00000986956906
7600 100.00000924988495 6.278970156979581 -99.8026872571317
7700 100.00000863261317 12.53324184506617 -99.21148912991777
7800 100.0000080202232 18.73805056458767 -98.228748669031
7900 100.00000741513185 24.868908774414006 -96.85834429410973
8000 100.00000681972718 30.90162078135972 -95.10568435709045
8100 100.00000623635901 36.81237822955507 -92.97768579802373
8200 100.00000566732959 42.577854061186905 -90.48274684717644
8300 100.0000051148845 48.17529457778811 -87.63071388115166
8400 100.00000458120421 53.582609238750955 -84.43284256382812
8500 100.00000406839483 58.7784578426692 -80.90175342547614
8600 100.0000035784801 63.74233474744264 -77.05138205535967
8700 100.00000311339362 68.45464979676468 -72.89692410438892
8800 100.00000267497069 72.89680563361073 -68.45477531487283
8900 100.00000226494173 77.05127109560578 -63.74246681404648
9000 100.00000188492496 80.90165040261006 -58.77859592673901
9100 100.00000153642003 84.4327478634704 -53.582752780233704
9200 100.00000122080235 87.63062784656721 -48.17544299097739
9300 100.00000093931756 90.48266977747897 -42.578006738258736
9400 100.00000069307654 92.97761794671167 -36.812534544236804
9500 100.000000483051 95.10562593092368 -30.901780092697216
9600 100.00000031006994 96.85829545233462 -24.869070430601987
9700 100.0000001748159 98.22870952295547 -18.73821390682732
9800 100.00000007782285 99.21145974283534 -12.533406211414873
9900 100.00000001947349 99.80266764458823 -6.279134886157624
10000 99.99999999999805 99.99999999996388 -8.268337527927994e-05

Hier ist eine Live-Version , die auf dem SageMathCell-Server läuft.

Wie Wikipedia erwähnt, gibt es Versionen höherer Ordnung von Leapfrog, die unabhängig voneinander in den 1980er und 1990er Jahren von mehreren Personen (Forest & Ruth, Yoshida und Candy & Rozmus) entdeckt wurden. Die Version 4. Ordnung ist dem einfachen Leapfrog oder Verlet sicherlich überlegen, aber die Verbesserungen gegenüber der erhöhten Anzahl von Berechnungen nehmen mit den höheren Ordnungen ab, sodass es nicht viel Sinn macht, höher als die 4. oder 8. Ordnung zu gehen.

Hier ist eine Version des obigen Codes mit Yoshida-Koeffizienten 4. Ordnung.

''' Simple orbit sim
    4th order Leapfrog with Yoshida coefficients
    Written by PM 2Ring 2017.05.20
    Added Yoshida 2022.02.12
'''

from math import pi, radians, cos, sin

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

# Yoshida 4th order coefficients
def yosh4():
    q = 2 ** (1 / 3)
    w1 = 1 / (2 - q)
    w0 = -q * w1
    c1 = w1 / 2
    c2 = (w0 + w1) / 2
    return (c1, c2, c2, c1), (w1, w0, w1)

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # The Yoshida coefficients, multiplied by the time step
        yc, yd = yosh4()
        self.yc = [k * delta_time for k in yc]
        self.yd = [k * delta_time for k in yd]

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the 4th order synchronised leapfrog algorithm
        '''
        x, y = self.x, self.y
        vx, vy = self.vx, self.vy

        # 4th order Leapfrog integration
        for c, d in zip(self.yc, self.yd):
            x += vx * c
            y += vy * c
            ax, ay = acc(x, y)
            vx += ax * d
            vy += ay * d
        c = self.yc[-1]
        x += vx * c
        y += vy * c

        self.vx, self.vy = vx, vy
        self.x, self.y = x, y

        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 1000
delta_time = period / steps
print('Step:', delta_time)
istep = steps // 100

# Initial direction. Use 0 for a circular orbit
th = radians(0)
vx, vy = v * sin(th), v * cos(th)

body = Body(radius,0, vx,vy, delta_time)
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

Bei nur 1000 Schritten pro Umlaufbahn ist der maximale Fehler im Radius auf 1,1 Teile pro Milliarde gesunken. Hier ist der SageMathCell-Link .

Hier ist ein Diagramm von 10 Perioden einer exzentrischen Umlaufbahn, die mit diesem Code erstellt wurde. Es verwendet denselben Zeitschritt von 1,0 mit einem Anfangswinkel von 60° zur Vertikalen.

Der Radiusfehler ist pro Umlaufbahn auf 5 Teile in zehn Millionen gestiegen, obwohl das in dieser Größenordnung nicht sichtbar ist.

Wenn wir den Zeitschritt auf 2,0 erhöhen, geht der Fehler auf 6 Teile pro Million und wird im Diagramm offensichtlich.

Es gibt also eine Reihe von Dingen, die hier unter der sprichwörtlichen Haube passieren können, die dazu führen können, dass diese Art von Fehler auftritt.

Zunächst einmal wird der Rundungsfehler das Herzstück all Ihrer Probleme sein. Nehmen wir an, Ihre Maschinenpräzision liegt bei etwa$10^{-16}$, was es ist, wenn Sie eine 64-Bit-Gleitkommadarstellung verwenden (dies ist auf der ganzen Linie ziemlich Standard, Sie müssen sich alle Mühe geben, dies zu ändern). Nehmen wir an, Sie lösen die ODE$\frac{d^2 r}{dt^2}=\frac{GM}{r^2}$und Sie zerlegen es in einzelne Komponenten usw. Ihre Entfernung beträgt 10 000 000 oder$10^7$, die dann quadriert wird, damit Sie arbeiten$10^{14}$am Ende dieser Gleichung,$1.9 \times 10^{27}$kg auf die Oberseite, und das alles multipliziert mit$6.67 \times 10^{-11}$. Wie Sie sich vorstellen können, können diese stark unterschiedlichen Zahlen Probleme mit Rundungsfehlern verursachen (wenn dies nicht klar ist, lassen Sie es mich wissen und ich kann ein konkreteres Beispiel hinzufügen). Sie werden beginnen, effektive Sig-Feigen zu verlieren, und dieser Fehler kann zu Dingen führen, die Sie sehen.

Außerdem sind nicht alle Löser gleich. Symplektische Löser sparen Energie, was großartig ist, aber in der Computerphysik gibt es kein kostenloses Mittagessen. sie haben ihren Preis. Die genaue Art dieser Kosten ist mir nicht klar, aber mir wurde immer geraten, sie zu vermeiden, es sei denn, Energieeinsparung auf granularer Ebene ist wichtig (meine aktuelle Forschung befasst sich mit N-Körper-relativistischen Objekten und mein Berater hat uns von symplektischen Integratoren abgeraten , er hat uns einen Grund gegeben, aber ich kann mich nicht genau erinnern, was es war, ich werde dem nachgehen)

Darüber hinaus werden nicht alle Löser gleich erstellt. Wenn Sie einen impliziten Löser finden, ist er langsamer, aber genauer, und je höher die Ordnung, desto besser. Einige Löser haben Fähigkeiten, um mit Steifheit und anderen Dingen umzugehen, während andere dies nicht tun; ab einem bestimmten Level wird es eher zu einer dunklen Kunst als zu einer Wissenschaft, weil die Gründe, warum manche mehr arbeiten als andere, immer düsterer werden.

tl; dr, versuchen Sie eine Reihe verschiedener Löser und skalieren Sie Ihre Einheiten näher an die Einheit (z. B. verwenden Sie die Masse als 1, um 1 Jupitermasse neu zu drucken, verwenden Sie etwas nahe der Einheit für die Entfernung und skalieren Sie G entsprechend neu. Spielen Sie damit herum, bis G innerhalb von a liegt auch einige Größenordnungen bis zur Einheit)

Was folgt, ist zu lang, um ein Kommentar zu sein. Da Kommentare gelöscht oder in den Chat verschoben werden können, mache ich dies außerdem zu einer Antwort.

Jeder numerische Integrator, der auf dem Fortschreiten des kartesischen Zustands von einem Zeitschritt zum nächsten basiert, hat eine optimale Schrittgröße, die davon abhängt

• Die Art des vorliegenden Problems,
• Die Einheiten, die verwendet werden, um das vorliegende Problem darzustellen,
• Die Mechanismen zur Darstellung der reellen Zahlen und
• Die Mechanismen, die verwendet werden, um den Status von einem Zeitschritt zum nächsten zu verschieben.

Ich werde mir die letzten beiden Punkte ansehen, wobei ich mich auf Integrationstechniken konzentriere, die die kartesische Position und Geschwindigkeit basierend auf der Beschleunigung verbessern.

Das numerische Integrieren der Position eines Fahrzeugs im Raum ist inhärent ein Problem einer gewöhnlichen Differentialgleichung zweiter Ordnung (ODE). Einige numerische Integrationstechniken (z. B. kanonisches RK4) können nur ODEs erster Ordnung lösen. Eine ODE zweiter Ordnung kann mit einem ODE-Löser erster Ordnung gelöst werden, indem Hilfsgleichungen hinzugefügt werden, aber dadurch wird die Geometrie aus dem Fenster geworfen. (Dass dies eine ODE zweiter Ordnung ist, ist Geometrie.) Zusätzlich zum Wegwerfen der Geometrie werfen Techniken wie RK4 auch die Erhaltungssätze aus dem Fenster.

Jede numerische Integrationstechnik, die die kartesische Geschwindigkeit und Position basierend auf der Beschleunigung vorantreibt, steht vor zwei grundlegenden Herausforderungen. Wenn der Zeitschritt extrem klein ist, werden Geschwindigkeit und/oder Position eingefroren. Das Problem mit extrem kleinen Zeitschritten ist das$1+10^{-16}$ist genau gleich eins in Arithmetik mit doppelter Genauigkeit. Auf der anderen Seite der Zeitschrittgröße treten Fehler auf, die der Technik selbst innewohnen, wenn der Zeitschritt extrem groß gemacht wird. Kein numerischer Integrator wird mit einer Schrittgröße von zehn Orbits gut funktionieren.

Für ein gegebenes Problem, einen gegebenen Satz von Einheiten, eine gegebene Darstellungstechnik der reellen Zahlen und eine gegebene Integrationstechnik gibt es eine optimale Schrittweite, die irgendwo zwischen zu klein und zu groß liegt. Diese optimale Schrittweite ist möglicherweise nicht konstant. Für eine elliptische Umlaufbahn ist es beispielsweise am besten, kleinere Schritte in der Nähe der Periapsis und größere Schritte in der Nähe der Apoapsis zu machen.

Diese optimale Schrittweite ist für Kreisbahnen etwas leicht zu finden, da dies die optimale Schrittweite konstant macht. Was man sehen wird, ist eine Badewannenfehlerkurve. Idealerweise sollten Integratoren nahe am tiefen Ende dieser Badewanne schwimmen.