2 body equations of motion in 3D

Question

How do you write the 2 body equations of motion in 3D as a system of differential equations

Answer 1

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

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

$a=\dfrac{\mu}{r^2}$ (Newton's universal law of gravitation / two body acceleration in scalar form)

Division be zero errors are avoided by using the vector equation for acceleration, since any one of the position components can be equal to zero, but the norm will be a positive number (assuming the norm itself isn't 0)

$\vec{a}=\dfrac{-\mu}{|\vec{r}|^3}\vec{r}$ (Newton's universal law of gravitation / two body acceleration in vector form. Note that the negative sign was added since gravitational acceleration is in the opposite direction of the position vector.

Here is the Python version of this ODE:

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

And a plot for fun:

Answer 2

This answer to How to solve the two-body problem in the ECI frame through numerical integration? says:

Instead of a single second order differential equation

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

We can solve the following pair of first order differential equations in parallel

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

using various simple methods including standard libraries or homegrown implementations of Runge-Kutta including my favorite simple RK4/5 with variable step size.

[...]

From this answer you can see a 2D implementation using not only the monopole $1/r^2$ gravity term but the additional quadrupole $J_2$ term for Earth's oblate shape and field. For more on that see this answer to Trouble deriving rectangular components of acceleration of satellite in orbit around Earth with J2 consideration.

Here $\mu$ with the standard gravitational parameter which would be the product $GM$ if we knew each parameter accurately. It turns out that for solar system bodies we can determine their product much more accurately than either one alone. If you like to work unitless and your orbiting body has negligible mass, set $\mu = 1$ and your period will be $2 \pi$.

So the question asks:

How do you write the 2 body equations of motion in 3D as a system of differential equations in the form Xdot = f(X).

Written out, $\dot{\mathbf{r}} = \mathbf{v}$ and $\dot{\mathbf{v}} = -\frac{\mu}{r^3}\mathbf{r}$ become:

$$\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}$$

If:

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

then

$$\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} $$