Расчет планет и лун на основе гравитационной силы Ньютона

Помимо числовых проблем, «С Солнцем в центре» может быть частью вашей проблемы. Получите все данные от Horizons относительно барицентра Солнечной системы, а не Солнца, которое движется относительно барицентра. Этот барицентр является инерциальной системой отсчета, а центр Солнца — нет. Также убедитесь, что вы вводите начальное положение и скорость Солнца и позволяете ему двигаться, если вы еще этого не сделали.

Я опубликую численные методы в этом ответе и полномасштабный расчет солнечной системы (включая теорию относительности и возможные эффекты сплюснутой формы Солнца) в качестве второго ответа. Слишком много, чтобы поместить все это в один ответ.

Метод, который вы описываете, выглядит как метод Эйлера или то, что можно назвать методом Эйлера вперед . На каждом временном шаге$dt$вы вычисляете все силы и их результирующие чистые ускорения$\mathbf{a}$, а затем просто увеличивайте скорость$\mathbf{v}$от$\mathbf{a} dt$и все позиции$\mathbf{x}$от$\mathbf{v} dt$. Вам нужны абсурдно крошечные временные шаги, чтобы приблизиться к этому. Вы упомянули 2 секунды и 250 миллисекунд, что немного меньше для шкалы времени Солнечной системы.

В приведенном ниже сценарии я написал метод Эйлера вперед и еще два метода низкого порядка Рунге-Кутты , обычно называемые RK2 и RK4. Для каждого метода рассчитывается упрощенная (фиктивная) орбита Меркурия вокруг фиксированного Солнца за 100 дней с количеством итераций от 10 до 10 000. Кроме того, для каждого я использую решатель библиотеки SciPyodeint() с относительным допуском точности 1E-12 на шаг. odeint — это оболочка Python lsodaиз библиотеки FORTRAN odepack.

Вы можете видеть, что даже RK4 соглашается с odeint уровнем метров через 100 дней. , если вы используете временной шаг около 15 минут, а метод Эйлера вперед (то, что вы используете) потребует абсурдного количества шагов, чтобы хотя бы приблизиться к этому.

Численные методы - не единственная проблема здесь, я опубликую второй ответ через несколько дней с остальным, что вам нужно. Я могу сделать это под отдельным вопросом, а не публиковать два ответа на один и тот же вопрос .

Но это должно заставить вас начать либо кодировать RK4 в java, либо искать числовую библиотеку java, такую ​​​​как Apache, упомянутая в комментариях.

Самый простой стандартный способ решить орбитальную задачу с помощью решателя ОДУ — поместить все декартовы координаты и скорости в один длинный вектор, назовем его$y$, и напишите одну функцию для производной по времени:

Итак, если бы у вас было два тела и три измерения, вектор$y$было бы:

с шестью элементами на тело. производная от$x_1$является$v_{x1}$, и производная от$v_{x1}$это ускорение$a_{x1}$из-за всех остальных тел.

Скажем, у вас есть одно тело в центральной силе со стандартным гравитационным параметром$GM$. Скорость изменения положения - это просто скорость,

а скорость изменения скорости есть ускорение, вызванное силой;

Если у вас было несколько тел,$\mathbf{r}$будет расстояние между парами тел, и для каждого тела вы суммируете все остальные, как вы уже описали.

Итак, если вы написали$f$, это было бы:

Метод Euler Forward, который, как я полагаю, вы используете, просто повторяет

с временными шагами$h$. Улучшенный метод RK2 (показан ниже) будет записан так:

а вездесущий метод RK4 записывается как

Вот раздел, который показывает сходство и различия между методом Эйлера Форвард (самый простой метод Рунге-Кутта) и двумя методами РК более высокого порядка. Это написано для ясности, а не скорости явно.

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