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


$$\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));...


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


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


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


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


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


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

  • $\begingroup$ 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. $\endgroup$ – StephenG Dec 24 '15 at 20:04
  • $\begingroup$ Look for the (or a) "moving stars around" web page. It has multiple explorations to the numerical treatment of multiple star systems. $\endgroup$ – LutzL Dec 24 '15 at 22:56

It's quite simple actually. You strat by wtiting the second order dfferential equation into first order equations. Let's start by considering the x ccordinate here, equation can be re-written as
$ \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$
which will give you set of two equations written as
$\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)$
and your initial conditions become
$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.
Hope this is helpful enough


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.