2 уравнения движения тела в 3D

$state=[rx,ry,rz,vx,vy,vz]$

$\dot{state}=[vx,vy,vz,ax,ay,az]$

$a=\dfrac{\mu}{r^2}$(Универсальный закон тяготения Ньютона / ускорение двух тел в скалярной форме)

Ошибок деления на ноль можно избежать, используя векторное уравнение для ускорения, поскольку любой из компонентов положения может быть равен нулю, но норма будет положительным числом (при условии, что сама норма не равна 0)

$\vec{a}=\dfrac{-\mu}{|\vec{r}|^3}\vec{r}$(Универсальный закон тяготения Ньютона / ускорение двух тел в векторной форме. Обратите внимание, что был добавлен отрицательный знак, поскольку ускорение свободного падения направлено в направлении, противоположном вектору положения.

Вот версия этого ODE для Python:

def two_body_ode( t, state, mu = pd.earth[ 'mu' ] ):
    r = state[ :3 ]
    a = -mu * r / np.linalg.norm( r ) ** 3

    return np.array( [
        state[ 3 ], state[ 4 ], state[ 5 ],
        a[ 0 ], a[ 1 ], a[ 2 ] ] )

И сюжет для развлечения:

Этот ответ на Как решить задачу двух тел в кадре ECI с помощью численного интегрирования? говорит:

Вместо одного дифференциального уравнения второго порядка

$$\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}$$

Мы можем решить следующую пару дифференциальных уравнений первого порядка параллельно

$$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$$ $$\dot{\mathbf{r}} = \mathbf{v}$$

используя различные простые методы, включая стандартные библиотеки или доморощенные реализации Рунге-Кутта, включая мой любимый простой RK4/5 с переменным размером шага .

[...]

Из этого ответа вы можете увидеть 2D-реализацию, использующую не только монополь$1/r^2$гравитационный член, но дополнительный квадруполь$J_2$термин для обозначения сплюснутой формы и поля Земли. Подробнее об этом см. в этом ответе на проблему с получением прямоугольных составляющих ускорения спутника на орбите вокруг Земли с учетом J2.

Здесь$\mu$со стандартным гравитационным параметром , который был бы произведением$GM$если бы мы точно знали каждый параметр. Оказывается, для тел Солнечной системы мы можем определить их произведение гораздо точнее, чем любое из них по отдельности. Если вам нравится работать без единиц измерения, а ваше вращающееся тело имеет незначительную массу, установите$\mu = 1$и ваш период будет$2 \pi$.

Итак, вопрос:

Как записать уравнения движения двух тел в 3D в виде системы дифференциальных уравнений в виде Xdot = f(X).

Выписано,$\dot{\mathbf{r}} = \mathbf{v}$и$\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$становиться:

$$\dot{x} = v_x$$ $$\dot{y} = v_y$$ $$\dot{z} = v_z$$ $$\dot{v_x} = -\frac{\mu}{r^3}x$$ $$\dot{v_y} = -\frac{\mu}{r^3}y$$ $$\dot{v_z} = -\frac{\mu}{r^3}z$$ $$r = \left(x^2 + y^2 + z^2 \right)^{1/2}$$

Если:

$$\mathbf{X} = \begin{bmatrix} x \\ y \\ z \\ v_x \\ v_y \\ v_z \\ \end{bmatrix}$$

затем

$$\mathbf{\dot{X}} = \begin{bmatrix} v_x \\ v_y \\ v_z \\ -\frac{x \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ -\frac{y \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ -\frac{z \ \mu}{(x^2 + y^2 + z^2)^{3/2}} \\ \end{bmatrix}$$