Die numerische Analyse wird verwendet, um Annäherungen an Dinge zu berechnen: den Wert einer Funktion an einem bestimmten Punkt, an dem sich eine Wurzel einer Gleichung befindet, oder die Lösungen für einen Satz von Differentialgleichungen. Es ist ein großes und wichtiges Thema, da in der Praxis die meisten realen Probleme in Mathematik, Naturwissenschaften und Technik keine explizite geschlossene Lösung haben werden (und selbst wenn dies der Fall ist, ist es möglicherweise nicht möglich, damit zu rechnen unendliche Präzision - Computer repräsentieren schließlich Zahlen mit endlicher Präzision). Im Allgemeinen gibt es Kompromisse zwischen Genauigkeit und Rechengeschwindigkeit.

Für das Drei-Körper-Problem haben wir drei Punktmassen an Startpositionen $ \ mathbf {x} _i (0) $ span> mit Geschwindigkeiten $ \ mathbf {v} _i (0) $ span>, den wir für spätere Zeiten berechnen möchten $ t $ span>. Mathematisch wollen wir die Lösung für das System $$ \ mathbf {x} '_ i (t) = \ mathbf {v} _i (t), $$ span> $$ \ mathbf {v} '_ i (t) = \ mathbf {f} _i (t) / m_i, $$ span> $$ \ mathbf {f} _i (t) = Gm_i \ sum_ {j \ neq i} \ frac {m_j (\ mathbf {x} _j- \ mathbf {x} _i)} {|| \ mathbf {x} _i - \ mathbf {x} _j || ^ 3}. $$ span>

Die naheliegende Methode besteht darin, zu denken, "wenn wir einen winzigen Schritt $ h $ span> in der Zeit vorwärts gehen, können wir alles als linear approximieren", also machen wir a Formel, in der wir den Zustand zum Zeitpunkt $ t + h $ span> aus dem Zustand zum Zeitpunkt $ t $ span> berechnen ( und so weiter für $ t + 2h $ span> und weiter): $$ \ mathbf {x} _i (t + h) = \ mathbf {x} _i (t) + h \ mathbf {v} _i (t), $$ span> $$ \ mathbf {v} _i (t + h ) = \ mathbf {v} _i (t) + h \ mathbf {f} _i (t). $$ span> Dies wird als Eulers Methode bezeichnet. Es ist einfach, neigt aber dazu, ungenau zu sein; Der Fehler pro Schritt ist $ \ approx O (h ^ 2) $ span> und sie neigen dazu, sich aufzubauen. Wenn Sie es bei einem Zweikörperproblem versuchen, führen die umlaufenden Massen aufgrund des Fehleraufbaus eine vorlaufende Rosettenumlaufbahn durch, insbesondere wenn sie sich nahe kommen.

Es gibt eine Reihe von Methoden zum numerischen Lösen von ODEs. Man kann Methoden höherer Ordnung verwenden, die die Funktionen an mehr Punkten abtasten und sie daher besser approximieren. Es gibt implizite Methoden, die, anstatt zu einem späteren Zeitpunkt zu versuchen, einen Zustand nur auf der Grundlage des aktuellen Zustands zu finden, nach einem selbstkonsistenten späten und mittleren Zustand suchen. Die schwerwiegendsten Methoden zum Lösen von ODEs reduzieren auch die Schrittgröße $ h $ span>, wenn die Kräfte bei engen Begegnungen groß werden, um sicherzustellen, dass die Genauigkeit akzeptabel bleibt. Wie gesagt, das ist ein großes Thema.

Für mechanische Simulationen sollten Sie sich jedoch insbesondere Methoden ansehen, mit denen Energie und andere konservierte Mengen eingespart werden können ( symplektische Methoden - diese werden von Fachleuten für langfristige Umlaufbahnberechnungen verwendet). . Am einfachsten ist vielleicht die semi-implizite Euler-Methode. Es gibt auch die Verlet-Methode und die Sprungintegration. Ich mag die semi-implizite Euler-Methode, weil sie sehr einfach ist (aber als Methode erster Ordnung immer noch nicht besonders genau ist): $$ \ mathbf {v} _i (t + h) = \ mathbf {v} _i (t) + h \ mathbf {f} _i (t), $$ span> $$ \ mathbf {x} _i ( t + h) = \ mathbf {x} _i (t) + h \ mathbf {v} _i (t + h). $$ span> Sehen Sie den Unterschied? Sie berechnen zuerst die aktualisierte Geschwindigkeit und verwenden sie dann, um die Positionen zu aktualisieren - ein winziger Trick, aber plötzlich verhalten sich 2-Körper-Bahnen gut.

Das Drei-Körper-Problem ist im wahrsten Sinne des Wortes chaotisch. Wir wissen, dass es Situationen gibt, in denen winzige Unterschiede in den Anfangsbedingungen auf beliebig große Unterschiede in späteren Positionen skaliert werden (selbst wenn wir sehr enge Durchgänge zwischen Massen ausschließen). Selbst mit einer willkürlich feinen numerischen Genauigkeit wird es einen Zeitpunkt geben, an dem unsere berechneten Umlaufbahnen völlig falsch sind. Der allgemeine "Stil" der Flugbahn mag immer noch korrekt sein, weshalb es in Ordnung ist, mit semi-implizitem Euler herumzuspielen, solange man keine Weltraummission basierend auf den Ergebnissen plant.

Die "Lösung" des Dreikörperproblems kann als Paar von Differentialgleichungen geschrieben werden. \ begin {align} \ vec {v} & = \ frac {\ mathrm d \ vec {x}} {\ mathrm dt} \\ \ vec {a} & = \ frac {\ mathrm d \ vec {v}} {\ mathrm dt} \ end {align} span> wo letzteres normalerweise in Form der Kraft geschrieben wird, $ m \ vec {a} = \ vec {F} $ span>. Verwenden Sie dann die Definition der Ableitung, $$ \ frac {\ mathrm df} {\ mathrm dx} = \ lim_ {h \ to0} \ frac {f (x + h) -f (x)} {h}, $$ span> Die Differentialgleichungen, mit denen wir begonnen haben, können als finite Differenz -Gleichungen geschrieben werden. \ begin {align} \ vec x (t + \ mathrm dt) & \ simeq \ vec x (t) + \ vec v \, \ mathrm dt \\ \ vec v (t + \ mathrm dt) & \ simeq \ vec v (t) + \ vec a \, \ mathrm dt \ end {align} span> Dies geschieht unter der Annahme, dass $ \ mathrm dt $ span> eher "nur klein" als infinitesimal ist.

Es sind diese Gleichungen, die ein Computer verwendet, um das Drei-Körper-Problem zu "lösen": Wenn eine Anfangsbedingung für jeden der Körper gegeben ist, tritt der Computer unter Verwendung der obigen Finite-Differenzen-Gleichungen iterativ zeitlich vorwärts Die Kraft ist die $ n $ span> -Körperkraft.

Wie in den Kommentaren erwähnt, heißt das vereinfachende Modell oben Euler-Integration und ist für dieses Problem überhaupt nicht gut geeignet, da es keine Energie spart. Eine bessere Option heißt Velocity Verlet , was in anderen Posts auf Physics.SE nicht erwähnt wird (ich würde diesen Post von mir empfehlen da es einige anständige, aber kurze Details der Implementierung enthält).

Ein einigermaßen zugänglicher Artikel über das Drei-Körper-Problem, einschließlich einer Diskussion über die numerischen Algorithmen, ist Musielak & Quarles (2014).Für eine umfassendere Diskussion nützlicherer Techniken höherer Ordnung, wie von homocomputeris vorgeschlagen, siehe Farrés et al (2012)

Vor einigen Jahren habe ich ein Modell des inneren Sonnensystems, der Sonne, des Erdmondes und der Planeten von Merkur bis Mars zusammengestellt.

Ich habe eine Leapfrog-Variante von Rk4 verwendet und anständige Ergebnisse erzielt. Die Umlaufzeiten stimmen mit den ursprünglichen Konfigurationen überein, die über die Zeit stabil blieben. Sogar einige der schlechten Dinge, die ich bekam, deuteten auf gute numerische Praktiken hin.

Ich würde meine Planeten mit anfänglichem Impuls initialisieren, wobei die Sonne keinen anfänglichen Impuls hat. Der Schwung der Planeten ließ nicht nach, so dass es diese seltsame Netzdrift in ungefähr gerader Richtung gab, die eine Weile dauerte, um herauszufinden.

Zusätzlich zum Aussortieren war eine wichtige Funktion die Reihenfolge, in der Änderungen an Position und Geschwindigkeit vorgenommen wurden.

Sie möchten nicht durch Ihren Körper iterieren, die Nettokraft in einem Durchgang Körper für Körper berechnen und dann Ihre dynamischen Änderungen Körper für Körper anwenden. Andernfalls wird beispielsweise die Position der Venus aktualisiert, während die Position des Mars unverändert bleibt, während Sie die Position der Erde aktualisieren. Ihre Zeitschritte sind nicht mehr synchron.

Speichern Sie irgendwo im Speicher die neuen Verschiebungen und Geschwindigkeiten, während Sie durch die Körper iterieren, und führen Sie Ihre Änderungen erst durch, nachdem Sie Ihre Dynamik beim aktuellen Tick gefunden haben. Entschuldigung und Glückwunsch, wenn Sie das beim ersten Versuch bemerkt haben. Es hat mich ein bisschen gestolpert.

Eine Möglichkeit, drei Körper zu simulieren, besteht in der Verwendung von Zeitschrittinkrementen, in denen Sie so tun, als würde jeder Körper für einen kurzen Zeitraum mit konstanter Geschwindigkeit beschleunigen. Sie beginnen mit den Anfangspositionsvektoren, Geschwindigkeitsvektoren und den Massen der Körper und verwenden sie dann, um den Beschleunigungsvektor jedes Körpers zu berechnen. In der Newtonschen Physik ist der Gravitationsbeschleunigungsvektor von Körper A aus Körper B durch die Gleichung $$ \ vec {a_ {GA}} = - \ frac {G \ left (\ vec {) gegeben A} - \ vec {B} \ right) M_B} {|| \ vec {A} - \ vec {B} || ^ 3} $$ span> mit $ G $ span> ist die Gravitationskonstante, $ \ vec {A} $ span> ist der Positionsvektor von Körper A, $ \ vec {B} $ span> ist der Positionsvektor von Körper B, $ M_B $ span> ist die Masse von Körper B und $ \ vec {a_ {GA}} $ span> ist der Gravitationsbeschleunigungsvektor von Körper A. Um den gesamten Gravitationsbeschleunigungsvektor von Körper A zu ermitteln, addieren Sie jeden der Gravitationsbeschleunigungsvektoren von Körper A. A von jedem der anderen Körper.

Sie finden die Positionsvektoren und Geschwindigkeitsvektoren am Ende des Zeitschritts $ n $ span> unter Verwendung der allgemeinen Formeln $ $ \ vec {x_n} = \ frac {1} {2} \ vec {a_ {n-1}} s ^ 2 + \ vec {v_ {n-1}} s + \ vec {x_ {n-1}} $$ span> $$ \ vec {v_n} = \ vec {a_ {n-1}} s + \ vec {v_ {n-1}} $$ span > wobei $ \ vec {a_ {n-1}} $ span> der Beschleunigungsvektor für den Zeitschritt $ n-1 $ ist span>, $ \ vec {v_ {n-1}} $ span> ist der Geschwindigkeitsvektor für den Zeitschritt $ n -1 $ span>, $ \ vec {x_ {n-1}} $ span> ist der Positionsvektor $ n -1 $ span>, $ s $ span> ist die Länge der Zeitschrittinkremente, $ \ vec {x_n} $ span> ist der Positionsvektor für den Zeitschritt $ n $ span> und $ \ vec {v_n} $ Spanne > ist der Geschwindigkeitsvektor für den Zeitschritt $ n $ span>. Bei jedem Zeitschritt berechnen Sie auch die Beschleunigungsvektoren neu, da sie nicht konstant sind.

Diese Methode zum Simulieren von $ n $ span> -Körpern funktioniert für jedes Kraftgesetz, dh Sie können sie zum Simulieren von $ verwenden n $ span> Körper für ein Kraftgesetz $ F = f (r) $ span> mit $ f (r) $ span> ist eine Funktion der Entfernung, obwohl im Fall des Kraftgesetzes $ F = f (r) $ span> die Gleichung ist, die ich zuvor zur Berechnung der Beschleunigung angegeben habe von Körper A aus Körper B wird $$ \ vec {a_ {GA}} = - \ frac {\ vec {A} - \ vec {B}} {|| \ vec { A} - \ vec {B} ||} GM_Bf \ left (|| \ vec {A} - \ vec {B} || \ right) $$ span>, aber die anderen Formeln, die ich erwähnt habe, einschließlich des Teils über Das Hinzufügen der Beschleunigungsvektoren bleibt gleich.

Es stellt sich heraus, dass man durch Kenntnis der paarweisen Trajektorien die Bewegungsgleichung von Körpern bestimmen kann. Die Bewegungsgleichungen sind geschrieben als $$ \ frac {d ^ 2 \ vec r ^ k} {ds ^ 2} = - G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n (\ vec r ^ k- \ vec r ^ n)} {| \ vec r ^ k- \ vec r ^ n | ^ 3} \ tag {1} $$ span> Wir lösen das Hilfsproblem der Paarwechselwirkung von Körpern mit reduzierter Trägheitsmasse $$ \ frac {m_n m_k} {\ sum_ {k = 1} ^ N m_k} \ frac {d ^ 2 \ vec R ^ {kn}} {ds ^ 2} = -G \ frac {m_n m_k \ vec R ^ {kn}} {| \ vec R ^ {kn} | ^ 3} \ tag {2} $$ span> Diese Gleichungen unterscheiden sich in verschiedenen Anfangsbedingungen. Subtrahieren Sie von Gleichung (2) Gleichung (1), wir erhalten $$ \ frac {d ^ 2 R_0 ^ k} {ds ^ 2} - \ frac {d ^ 2 r ^ k} {ds ^ 2} = G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n (\ vec r ^ k- \ vec r ^ n)} {| \ vec r ^ k- \ vec r ^ n | ^ 3} - G \ sum_ {n = 1, n \ ne k} ^ N \ frac {m_n \ vec R ^ {kn}} {| \ vec R ^ {kn} | ^ 3} = 0, \ vec R ^ {kn} = \ vec r ^ k- \ vec r ^ n $$ span> Wobei $ \ vec R_0 ^ k = \ sum_ {n = 1 n \ ne k} ^ N \ frac {m_n \ vec R ^ {kn}} {\ sum_ {p = 1 } ^ N m_p} $ span> Holen Sie sich die Formel $$ \ frac {d ^ 2 \ vec R_0 ^ k} {ds ^ 2} = \ frac {d ^ 2 \ vec r ^ k} {ds ^ 2} $$ Durch Integration dieser Gleichheit erhalten wir die Bewegungsgleichung für jeden der N Körper. $$ \ vec r_0 ^ k = \ frac {d \ vec r ^ k (0)} {ds} s + \ vec r ^ k (0) + \ vec R_0 ^ k ( s) - \ frac {d \ vec R ^ k (0)} {ds} s- \ vec R ^ k (0) = \ vec R_0 ^ k (s) $$ span> Die Bewegung jedes Körpers wird durch die Bewegung des Trägheitszentrums des Paarkörpersystems bestimmt. Die Anfangsbedingungen werden aus der Bedingung $$ \ frac {d \ vecr ^ k (0)} {ds} s + \ vecr ^ k (0) = \ frac {d \ bestimmt vec R ^ k (0)} {ds} s + \ vec R ^ k (0) \ tag {3} $$ span> Wir wählen das Koordinatensystem $ \ sum_ {n = 1} ^ N m_n \ vec r ^ n (0) = 0, \ sum_ {n = 1} ^ N m_n \ frac { d \ vec r ^ n (0)} {ds} = 0 $ span> und dann ist Bedingung (3) identisch erfüllt. Erhält $ \ vec R ^ {kn} (0) = \ vec r ^ k (0) - \ vec r ^ n (0), \ frac {d \ vec R ^ {kn}} {ds} (0) = \ frac {d \ vec r ^ k} {ds} (0) - \ frac {d \ vec r ^ n} {ds} (0) $ span> Die Die Energie der Paartrajektorie des Systems wird nach der Formel $$ E_ {kn} = \ frac {m_k m_n} {\ sum_ {q = 1} ^ N m_q} (\) berechnet frac {d \ vec R ^ {kn}} {ds} (0)) ^ 2 / 2- \ frac {G m_k m_n} {| \ vec R ^ {kn} (0) |} + \ frac {M_ { kn} ^ 2 \ sum_ {q = 1} ^ N m_q} {2 m_k m_n | R ^ {kn} | ^ 2} = \ frac {m_k m_n} {\ sum_ {q = 1} ^ N m_q} (\ frac {d \ vec R ^ {kn}} {ds} (0)) ^ 2 / 2- \ frac {G m_k m_n} {| \ vec R ^ {kn} (0) |} + \ frac {m_k m_n } {2 \ sum_ {q = 1} ^ N m_q | R ^ {kn} (0) | ^ 2} (\ frac {d \ vec R ^ {kn}} {ds} (0) \ times \ vec R. ^ {kn} (0)) ^ 2 $$ span> Der Moment des Impulses wird durch das Formelvektorprodukt $$ M_ {kn} = \ frac {m_k m_n bestimmt } {\ sum_ {q = 1} ^ N m_q} \ frac {d \ vec R ^ {kn}} {ds} (0) \ times \ vec R ^ {kn} (0) $$ span> In diesem Fall wird die konstruierte Lösung in die berechneten Trajektorien $ \ vec r_0 ^ k- \ vec r_0 ^ p $ span> der Planeten eingesetzt und $ \ vec R ^ {kp} = r ^ kr ^ p $ span>, dann erhalten wir $ \ vec R ^ {kp} = \ vec r_0 ^ k- \ vec r_0 ^ p $ span>. Tatsächlich $$ r_0 ^ k-r_0 ^ p = \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {kn}} {\ sum_ {q = 1} ^ N m_q} - \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {pn}} {\ sum_ {q = 1} ^ N m_q} = \ sum_ {n = 1} ^ N \ frac { m_n (\ vec r ^ k- \ vec r ^ n)} {\ sum_ {q = 1} ^ N m_q} - \ sum_ {n = 1} ^ N \ frac {m_n (\ vec r ^ p- \ vec r ^ n)} {\ sum_ {q = 1} ^ N m_q} = \ vec r ^ k- \ vec r ^ p = R_ {kp} $$ span> Vereinfachte Formel $$ \ vec r_0 ^ u = \ sum_ {n = 1} ^ N \ frac {m_n \ vec R ^ {un} (s)} {\ sum_ {q = 1} ^ N m_q | R ^ {un} (s) |} \ frac {p_ {un}} {1 + e_ {un} cos (\ varphi_ {un} - \ varphi_ {un0})} $$ span> Der Winkel wird durch die Formel $ \ varphi_ {un} - \ varphi_ {un0} = \ int_0 ^ s \ frac {M_ {un}} {| \ vec R ^ {pn berechnet } (u) | ^ 2} du $ span> Die Größe $$ \ frac {R ^ {un} (\ varphi)} {| R ^ {un} (\ varphi) |} = \ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} (0) |} cos \ varphi + \ frac {\ vec R ^ {un} (0) \ times \ vec M_ {un}} {| \ vec R ^ {un} (0) \ times \ vec M_ {un} |} sin \ varphi $$ span> des Winkels $ \ varphi $ span>, das allen Paartrajektorien gemeinsam ist, entspricht der Richtung $ \ frac {R ^ {un} (\ varphi)} {| R ^ {un} (\ varphi) |} $ span> mit einzelnen orths $$ \ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} ( 0) |}, \ frac {\ vec R ^ {un} (0) \ times \ vec M_ {un}} {| \ vec R ^ {un} (0) \ times \ vec M_ {un} |} $ $ span>. Wir haben die endgültige Formel für die Flugbahn von Körpern $$ \ vec r_0 ^ u (\ varphi) = \ sum_ {n = 1} ^ N \ frac {m_n} {\ sum_ { q = 1} ^ N m_q} [\ frac {\ vec R ^ {un} (0)} {| \ vec R ^ {un} (0) |} cos \ varphi + \ frac {\ vec R ^ {un} (0) \ times \ vec M_ {un}} {| \ vec R ^ {un} (0) \ times \ vec M_ {un} |} sin \ varphi] \ frac {p_ {un}} {1 + e_ {un} cos \ varphi} $$ span> Da es jedoch auch für kurze Zeit keine unendliche Lösung gibt, muss $$ | \ vec R ^ {angenommen werden un} (\ varphi) | = \ frac {p_ {un}} {1 + e_ {un} cos \ varphi} >0, e_ {un} <1 $$ span> Dies führt zu einer Einschränkung des Werts der Parameter ;; Die Energie der Paarwechselwirkung muss negativ sein. Somit werden Körper mit geringer Masse ausgeschlossen, deren Beitrag zur Energie des Systems nicht signifikant ist. Nach der Formel der Paarwechselwirkung variiert der Radius nach dem Gesetz, siehe Landau Lifshits Band 1 $$ r = a (ecosh \ xi-1); t = \ sqrt {\ frac {a ^ 3} {G \ sum_ {q = 1} ^ N m_q}} ( esinh \ xi- \ xi), r \ approx c \ sqrt {\ frac {r_g} {a}} t $$ span> Nach einfachen Transformationen wird diese Formel auf $$ r = \ sqrt {V ^ 2 (0) + \ frac {(\ vec V (0) \ times \ vec R (0) reduziert)) ^ 2} {R (0) ^ 2} - \ frac {2G \ sum_ {q = 1} ^ N m_q} {R (0)}} t $$ span> Diese lineare Zunahme des Radius verschiebt dieSchwerpunkt des Systems.Ein Teilchen, das sich ins Unendliche bewegt, führt dazu, dass das gesamte System mitgeführt wird und sich der Schwerpunkt verschiebt.Die zweite kosmische Geschwindigkeit entspricht einem positiven Wert des Radikanden oder einer positiven Energie des Systems