$5\mathrm{.}2\mathrm{.}3$ Gravitational systems

The motion of $n$ planets $(P_0,P_1,$dots,$P_{n-1})$ is governed by the set of $6n$ first order ODEs

$r_i^{\text{'}}=v_i,$ and $,v_i^{\text{'}}={\displaystyle \underset{j=1ji}{\overset{n}{}}}G{m}_j\frac{r_j-r_i}{d_{ij}^3}$

where $m_i,r_i(t)$ and $v_i(t)$ are the mass, position vector and velocity of the planet $P_i,$ respectively. The scalar $d_{ij}=d_{ji}=||r_j-r_i||$ is the distance between the two planets $P_i$ and $P_j.$ The gravitational constant, $G,$ can be taken equal to $1$ without loss of generality.

$5\mathrm{.}2\mathrm{.}3\mathrm{.}1$ The two$-$body problem

Two point particles $($or two planetary bodies$)$ interacting only with each other under Newton's laws of gravitational attraction are well understood and the equation of motion can be integrated analytically. However, here, the trajectories of the two planets will be computed numerically.

$(1)$ Write down, in vector form, the four ODEs governing the motion of the two planets interacting gravitationally. $($If written in component form, these four equations would give $12$ equations.$)$

$(2)$ Write a Python function two$\_$body$($x$)$ which returns the vector $f(x),$ as defined in $(5\mathrm{.}1),$ for the two$-$body problem. The NumPy array $x$ must be of size $12,$ and its shape must be $2,2,3.$ For the entry $x[i,j,k],i=0,1$ corresponds to the first or second planet; $j=0,1$ corresponds to the position vector or the velocity of the planet considered; and $k=0,1,2$ corresponds to the $x,y$ or $z$ component of the vector $($position or velocity$).$ You can assume that the masses of the planets, $m_0$ and $m_1,$ are defined outside the function two$\_$body.

$(3)$ Consider two planets, $P_0$ and $P_1,$ of mass $m_0=1$ and $m_1=2,$ respectively. At time $t=0,P_0$ is at $r_0=(1,0,0)$ with velocity $v_0=(0,\frac{1}{2},0)$ and $P_1$ is at $r_1=(-1,0,0)$ with velocity $v_1=(0,-\frac{1}{4},0).$ The trajectories of the two planets are periodic and remain in the $(x,y)-$plane. Solve the equations of motion numerically, using your function rk$4,$ then plot the complete orbits of the two planets.