# System of vector differential equations

How to solve numerically this system of vector differential equations (along with its initial conditions) without using any physical trick (like the reduced mass $\mu=\frac{m_1m_2}{m_1+m_2}$):

$$G\frac{\require{cancel}\cancel{m_1}m_2}{||\vec{r}_1-\vec{r}_2||^2}\frac{\vec{r}_2-\vec{r}_1}{||\vec{r}_1-\vec{r}_2||}=\require{cancel}\cancel{m_1}\frac{d^2\vec{r}_1}{dt^2} \tag{1}$$

$$G\frac{m_1\require{cancel}\cancel{m_2}}{||\vec{r}_1-\vec{r}_2||^2}\frac{\vec{r}_1-\vec{r}_2}{||\vec{r}_1-\vec{r}_2||}=\require{cancel}\cancel{m_2}\frac{d^2\vec{r}_2}{dt^2} \tag{2}$$

where

$$\vec{r}_1=\vec{r}_1(t)=x_1(t)\hat{i}+y_1\hat{j}+z_1(t)\hat{k}$$ $$\vec{r}_2=\vec{r}_2(t)=x_2(t)\hat{i}+y_2\hat{j}+z_2(t)\hat{k}$$

This is my attempt to solve it:

1. Separate into scalar equations:

$$\left\{\begin{matrix} \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(x_2-x_1)=\ddot{x}_1\\ \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(y_2-y_1)=\ddot{y}_1 \\ \frac{Gm_2}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(z_2-z_1)=\ddot{z}_1 \end{matrix}\right.$$

$$\left\{\begin{matrix} \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(x_1-x_2)=\ddot{x}_2 \\ \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(y_1-y_2)=\ddot{y}_2 \\ \frac{Gm_1}{[(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2)]^{3/2}}(z_1-z_2)=\ddot{z}_2 \end{matrix}\right.$$

1. Reduce order:
• I define the new variables $$\begin{matrix} u_1 \triangleq x_1 & u_2 \triangleq \dot{x}_{1}=\dot{u}_{1} & \dot{u}_{2}=\ddot{x}_1\\ u_3 \triangleq y_1 & u_4 \triangleq \dot{y}_{1}=\dot{u}_{3} & \dot{u}_{4}=\ddot{y}_1\\ u_5 \triangleq z_1 & u_6 \triangleq \dot{x}_{1}=\dot{u}_{5} & \dot{u}_{6}=\ddot{x}_1\\ u_7 \triangleq x_2 & u_8 \triangleq \dot{x}_{2}=\dot{u}_{7} & \dot{u}_{8}=\ddot{x}_2\\ u_9 \triangleq y_2 & u_{10} \triangleq \dot{y}_{2}=\dot{u}_{9} & \dot{u}_{10}=\ddot{y}_2\\ u_{11} \triangleq z_2 & u_{12} \triangleq \dot{z}_{2}=\dot{u}_{11} & \dot{u}_{12}=\ddot{z}_2\\ \end{matrix}$$
• I write the system of diferential equations in the form $\vec{u}'(t)=\vec{F}(\vec{u},t)$ with

$$\vec{u}'(t)=\begin{bmatrix} \dot{u}_{1}\\ \dot{u}_{2}\\ \dot{u}_{3}\\ \dot{u}_{4}\\ \dot{u}_{5}\\ \dot{u}_{6}\\ \dot{u}_{7}\\ \dot{u}_{8}\\ \dot{u}_{9}\\ \dot{u}_{10}\\ \dot{u}_{11}\\ \dot{u}_{12}\\ \end{bmatrix}$$

$$\vec{F}(\vec{u},t)=\begin{bmatrix} {u}_{2}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_7-u_1)\\ {u}_{4}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_9-u_3)\\ u_{6}\\ \frac{Gm_2}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_{11}-u_5)\\ u_{8}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_1-u_7)\\ u_{10}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_3-u_9)\\ u_{12}\\ \frac{Gm_1}{[(u_1-u_7)^2+(u_3-u_9)^2+(u_5-u_{11})^2)]^{3/2}}(u_5-u_{11})\\ \end{bmatrix}$$

1. Solve the system using RK4(5), for example, in Matlab, using the function ode45, with the initial condition vector $\vec{u}_0$ and the time span $tspan$ from $t0=0$ to $tmax=31536000$ (1 year).

G=6.67408e-11; %gravitational constant m^3 kg^-1 s^-2

m1=1.98855e30; %mass sun kg

m2=5.97237e24; %mass earth kg

r0=149.6e9; %distance sun-earth m

v0=29.78e3; %earth-to-sun speedm/s

r = @(U) ((U(1)-U(7))^2+(U(3)-U(9))^2+(U(5)-U(11))^2)^(3/2);

dUUdt = @(t,U,G,m1,m2) [U(2);...

G * m2 / r(U) * (U(7)-U(1));...

U(4);...

G * m2 / r(U) * (U(9)-U(3));...

U(6);...

G * m2 / r(U) * (U(11)-U(5));...

U(8);...

G * m1 / r(U) * (U(1)-U(7));...

U(10);...

G * m1 / r(U) * (U(3)-U(9));...

U(12);...

G * m1 / r(U) * (U(5)-U(11))];

dUdt = @(t,U) dUUdt(t,U,G,m1,m2);

U0 = zeros(12,1);

U0(7) = r0; %x0

U0(10) = v0; %vy0

tspan=0:86400:31536000; % 1 year in seconds with a step of 1 day

[T,Y]=ode45(dUdt,tspan,U0);

plot(Y(:,7),Y(:,9)) % corroboration that orbit is aprox circular

(Y(:,7).^2+Y(:,9).^2).^(1/2) %corroboration that earth-sun distance is aprox constant

(Y(:,8).^2+Y(:,10).^2).^(1/2) %corroboration that earth-to-sun speed is aprox constant

• The two body problem in Newtonian gravitation has been extensively studied and there is little purpose to solving it numerically, given the existence of exact, closed form solutions. See Wikipedia on the Central Force Problem. So why numerical ? The numerical is less accurate. – StephenG Dec 24 '15 at 20:04
• Look for the (or a) "moving stars around" web page. It has multiple explorations to the numerical treatment of multiple star systems. – LutzL Dec 24 '15 at 22:56

$\vec {\ddot x} = \dfrac{1}{m}\bigg(G\dfrac{m_1m_2}{(||\vec{x_1}-\vec{x_2}||)^2}-\dfrac{\vec{x_2}-\vec{x_1}}{||\vec{x_1}-\vec{x_2}||} \bigg)$
Let $\dot x= v_x$ and $\ddot x = \dot v_x$
$\dot x= v_x$
$\dot v_x = \dfrac{1}{m}\bigg(G\dfrac{m_1m_2}{(||\vec{x_1}-\vec{x_2}||)^2}-\dfrac{\vec{x_2}-\vec{x_1}}{||\vec{x_1}-\vec{x_2}||} \bigg)$
$x(0) = x_0\\ \dot v_x(0) = \dot x_0$
same thing can be done for $y$ and $z$ coordinates for each equation, and then you can use any forward integration scheme such as Euler's or Runge-Kutta to integrate numerically till desired time.