Cálculo de los planetas y lunas en base a la fuerza gravitatoria de Newton

Aparte de los problemas numéricos, "Con el Sol como centro" puede ser parte de su problema. Obtenga todos los datos de Horizons en relación con el Baricentro del Sistema Solar, no el Sol, que se mueve en relación con el baricentro. Ese baricentro es un marco de referencia inercial, mientras que el centro del Sol no lo es. También asegúrese de colocar la posición inicial y la velocidad del Sol y dejar que se mueva, si aún no lo ha hecho.

Publicaré los métodos numéricos en esta respuesta y el cálculo completo del sistema solar (incluida la relatividad y los posibles efectos de la forma achatada del Sol) como una segunda respuesta. Hay demasiado para ponerlo todo en una respuesta.

El método que describe se parece al método de Euler o lo que se puede llamar el método Euler Forward . En cada paso de tiempo$dt$usted calcula todas las fuerzas y sus aceleraciones netas resultantes$\mathbf{a}$, y luego simplemente incrementar las velocidades$\mathbf{v}$por$\mathbf{a} dt$y todas las posiciones$\mathbf{x}$por$\mathbf{v} dt$. Necesitas pasos de tiempo absurdamente pequeños para acercarte a esto. Usted mencionó 2 segundos y 250 milisegundos, un poco corto para las escalas de tiempo del sistema solar.

En el siguiente script, he escrito el método Euler Forward y otros dos métodos de orden bajo de Runge-Kutta , generalmente llamados RK2 y RK4. Para cada método, se calcula una órbita simplificada (falsa) de Mercurio alrededor de un Sol fijo durante 100 días, con un número de iteraciones que varía de 10 a 10 000. Además, para cada uno, uso el solucionador de la biblioteca SciPyodeint() con una tolerancia de precisión relativa de 1E-12 por paso. odeint es un contenedor de python lsodade la biblioteca FORTRAN odepack.

Puede ver que incluso RK4 está de acuerdo con odeint el nivel de metros después de 100 días si usa un paso de tiempo de aproximadamente 15 minutos, y el método Euler Forward (lo que usa) necesitará una cantidad absurda de pasos para siquiera acercarse a esto.

Las técnicas numéricas no son el único problema aquí, publicaré una segunda respuesta en unos días con el resto de lo que necesita. Puedo hacerlo bajo una pregunta separada en lugar de publicar dos respuestas a la misma pregunta .

Pero esto debería ayudarlo a comenzar a codificar RK4 en Java o buscar una biblioteca numérica de Java como la de Apache mencionada en los comentarios.

La forma estándar más sencilla de resolver un problema orbital usando un solucionador ODE es poner todas las coordenadas y velocidades cartesianas en un vector largo, llámelo$y$, y escribe una sola función para la derivada del tiempo:

Entonces, si tuvieras dos cuerpos y tres dimensiones, el vector$y$sería:

teniendo seis elementos por cuerpo. el derivado de$x_1$es$v_{x1}$, y la derivada de$v_{x1}$es la aceleración$a_{x1}$debido a todos los otros cuerpos.

Digamos que tienes un cuerpo en una fuerza central con un parámetro gravitacional estándar$GM$. La tasa de cambio de posición es solo la velocidad,

y la razón de cambio de la velocidad es la aceleración, debida a la fuerza;

Si tuviera varios cuerpos, el$\mathbf{r}$sería la distancia entre pares de cuerpos y para cada cuerpo sumarías todos los demás como ya lo describiste.

Así que si escribiste$f$, podría ser:

El método Euler Forward que creo que está utilizando simplemente itera

con pasos de tiempo de$h$. El método RK2 mejorado (que se muestra a continuación) se escribiría:

y el omnipresente método RK4 se escribe como

Aquí está la sección que muestra la similitud y las diferencias entre el método Euler Forward (el método Runge-Kutta más simple) y dos métodos RK de orden superior. Esto está escrito para mayor claridad, no para acelerar, obviamente.

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