Berechnung der Planeten und Monde nach Newtons Gravitationskraft

Abgesehen von numerischen Problemen kann "Mit Sonne als Zentrum" Teil Ihres Problems sein. Holen Sie sich alle Daten von Horizons relativ zum Baryzentrum des Sonnensystems, nicht zur Sonne, die sich relativ zum Baryzentrum bewegt. Dieses Baryzentrum ist ein Trägheitsbezugssystem, das Zentrum der Sonne hingegen nicht. Stellen Sie außerdem sicher, dass Sie die anfängliche Position und Geschwindigkeit der Sonne eingeben und sie sich bewegen lassen, wenn Sie dies nicht bereits tun.

Ich werde die numerischen Methoden in dieser Antwort und die vollständige Berechnung des Sonnensystems (einschließlich Relativitätstheorie und möglicher Auswirkungen der abgeflachten Form der Sonne) als zweite Antwort veröffentlichen. Es ist einfach zu viel, um alles in einer Antwort zusammenzufassen.

Die Methode, die Sie beschreiben, sieht aus wie die Euler-Methode oder die so genannte Euler-Forward-Methode . Bei jedem Zeitschritt$dt$Sie berechnen alle Kräfte und ihre resultierenden Nettobeschleunigungen$\mathbf{a}$, und erhöhen Sie dann einfach die Geschwindigkeiten$\mathbf{v}$von$\mathbf{a} dt$und alle Positionen$\mathbf{x}$von$\mathbf{v} dt$. Sie brauchen absurd kleine Zeitschritte, um dem auch nur nahe zu kommen. Sie haben 2 Sekunden und 250 Millisekunden erwähnt, ein bisschen kurz für die Zeitskalen des Sonnensystems.

Im folgenden Skript habe ich die Euler-Forward-Methode und zwei weitere Runge-Kutta- Methoden niedriger Ordnung geschrieben, die normalerweise RK2 und RK4 genannt werden. Für jede Methode wird eine vereinfachte (falsche) Merkurumlaufbahn um eine feste Sonne für 100 Tage berechnet, wobei die Anzahl der Iterationen zwischen 10 und 10.000 variiert. Außerdem verwende ich für jeden den SciPy - Bibliothekslöser odeint()mit einer relativen Genauigkeitstoleranz von 1E-12 pro Schritt. odint ist ein Python-Wrapper für lsodadie FORTRAN-Bibliothek odepack.

Sie können sehen, dass sogar RK4 odeint nach 100 Tagen mit dem Pegel von Metern übereinstimmt, wenn Sie einen Zeitschritt von etwa 15 Minuten verwenden, und die Euler-Forward-Methode (die Sie verwenden) eine absurde Anzahl von Schritten benötigt, um sich diesem überhaupt anzunähern.

Numerische Techniken sind hier nicht das einzige Problem, ich werde in ein paar Tagen eine zweite Antwort mit dem Rest veröffentlichen, was Sie brauchen. Ich kann dies unter einer separaten Frage tun, anstatt zwei Antworten auf dieselbe Frage zu posten .

Aber das sollte Sie dazu bringen, entweder RK4 in Java zu programmieren oder nach einer numerischen Java-Bibliothek wie der in den Kommentaren erwähnten Apache zu suchen.

Der einfachste Standardweg zum Lösen eines Orbitalproblems mit einem ODE-Löser besteht darin, alle kartesischen Koordinaten und Geschwindigkeiten in einen langen Vektor zu packen, nennen Sie ihn$y$, und schreiben Sie eine einzelne Funktion für die Zeitableitung:

Wenn Sie also zwei Körper und drei Dimensionen hätten, der Vektor$y$wäre:

mit sechs Elementen pro Körper. Die Ableitung von$x_1$ist$v_{x1}$, und die Ableitung von$v_{x1}$ist die Beschleunigung$a_{x1}$wegen all der anderen Körper.

Angenommen, Sie haben einen Körper in einer zentralen Kraft mit einem Standard-Gravitationsparameter$GM$. Die Änderungsrate der Position ist nur die Geschwindigkeit,

und die Geschwindigkeitsänderungsrate ist die Beschleunigung aufgrund der Kraft;

Wenn Sie mehrere Körper hatten, die$\mathbf{r}$wäre der Abstand zwischen Körperpaaren und für jeden Körper würden Sie alle anderen summieren, wie Sie bereits beschrieben haben.

Also, wenn Sie ausgeschrieben haben$f$, es wäre:

Die Euler-Forward-Methode , von der ich glaube, dass Sie sie verwenden, wird nur iteriert

mit Zeitschritten von$h$. Die verbesserte RK2-Methode (unten gezeigt) würde geschrieben werden:

und die allgegenwärtige RK4-Methode wird geschrieben als

Hier ist der Abschnitt, der die Ähnlichkeiten und Unterschiede zwischen der Euler-Forward-Methode (der einfachsten Runge-Kutta-Methode) und zwei RK-Methoden höherer Ordnung zeigt. Dies wurde aus Gründen der Klarheit geschrieben, offensichtlich nicht der Geschwindigkeit.

def deriv(t, X):

    x, v = X.reshape(2, -1)
    acc  = -GMs * x * ((x**2).sum())**-1.5

    return np.hstack((v, acc))

def derivv(X, t):
    return deriv(t, X) # historical reasons

def RK_all(F, t0, X0, n, h, method=None):

    hov2, hov6   = h/2.0, h/6.0
    t, X         = t0, X0
    answer, time = [], []

    answer.append(X)
    time.append(t)

    for i in range(n):
        k1 = F(t,        X            )
        k2 = F(t + hov2, X + hov2 * k1)
        k3 = F(t + hov2, X + hov2 * k2)
        k4 = F(t + h,    X + h    * k3)

        if   method == 'EulerFwd':
            X = X + h*k1        # X + h*F(t, X)

        elif method == 'RK2':
            X = X + h*k2

        elif method == 'RK4':
            X = X + hov6*(k1 + 2.*(k2+k3) + k4)

        else:
            pass

        t += h

        answer.append(X)
        time.append(t)

    return np.array(time), np.array(answer)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

GMs = 1.327E+20 # approx

X0 = np.array([-2.094E+10,  4.303E+10,  5.412E+09,
               -5.328E+04, -2.011E+04,  3.243E+03]) # approx

methodnames = ['RK4', 'RK2', 'EulerFwd']
niterations = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]

Time = 100*24*3600.  # total time

methdict = dict()

for methodname in methodnames:

    times, answers, ODEanswers, posdevs = [], [], [], []

    for n in niterations:

        h  = Time/float(n)
        t0 = 0.0

        time, answer    = RK_all(deriv, t0, X0, n, h, method=methodname)
        # recalculate using library ODE solver for same times, to compare
        ODEanswer, info = ODEint(derivv, X0, time,
                                 rtol=1E-12, full_output=True)
        posdev = np.sqrt((((answer - ODEanswer)[:,:3])**2).sum(axis=1))


        times.append(time)
        answers.append(answer)
        ODEanswers.append(ODEanswer)
        posdevs.append(posdev)

    methdict[methodname] = (times, answers, ODEanswers, posdevs)

if 1 == 1:
    plt.figure()
    for i, meth in enumerate(methodnames):
        plt.subplot(1, 3, i+1)
        for time, answer, ODEanswer, posdev in zip(*methdict[meth]):
            x, y, z = answer.T[:3]
            plt.plot(x, y)
        plt.ylim(-2.8E+11, 0.8E+11)
        plt.xlim(-1.2E+11, 0.8E+11)
        plt.title(meth, fontsize=16)
        plt.plot([0],[0], 'ok')
    plt.show()

if 1 == 1:
    plt.figure()
    for i, meth in enumerate(methodnames):
        plt.subplot(1, 3, i+1)
        for time, answer, ODEanswer, posdev in zip(*methdict[meth]):
            plt.plot(time/(24*3600.), posdev)
        plt.yscale('log')
        plt.ylim(1E-01, 1E+12)
        plt.title(meth+' vs odeint', fontsize=16)
    plt.suptitle('RKmethod - odeint (meters) vs time (days)', fontsize=18)
    plt.xticks([0, 20, 40, 60, 80, 100])
    plt.show()