# Creating a Gravitational 2-body simulation

I am trying to create a simulation for a gravitational 2 body problem. But I am kind of having trouble to define the equations that can be solve numerically. From an inertial frame I defined the position of the two objects as the $$\vec{r}_1$$ and $$\vec{r}_2$$ with masses $$m_1$$ and $$m_2$$.

Let the $$\vec{R}_{CM}$$ be the position of the CM of the objects. Now from the perspective of the CM, we can write position vectors of the objects in terms of $$\vec{r}'_1$$ and $$\vec{r}'_2$$.

In this case

$$\vec{r}'_1 = \frac{-m_2}{m_1 + m_2} \vec{r} \tag{1}$$

and $$\vec{r}'_2 = \frac{m_1}{m_1 + m_2} \vec{r} \tag{2}$$where $$\vec{r}= \vec{r}'_2 - \vec{r}'_1$$

Now in this case we can use the reduced mass and define the force on this mass. So we have,

$$\vec{F} = \mu \ddot{\vec{r}} = -Gm_1\mu / r^2 \vec{r}$$

Now I need to solve this equation and put back into the (1) and (2) right ?

• Your final equation has the wrong sign, the wrong mass, and the wrong exponent, so you should go back and check your derivation of it. Sep 28, 2019 at 20:47
• For a numerical simulation, you don’t necessarily have to use coordinates relative to the center of mass. Sep 28, 2019 at 20:53
• @G.Smith I need to find positions of the each particle. How can I define the force on the object of the reduced mass than ? Sep 28, 2019 at 21:20
• Your final equation is still wrong. Please see my answer. The fact that you have an $m_1$ without an $m_2$ is a clue that something is wrong, since $m_1$ and $m_2$ enter in a symmetrical way (except for sign) in $\vec{r}$. Also, you don't seem to understand that a radial inverse-square force can be written as $\hat{r}/r^2$ or $\vec{r}/r^3$ but not as $\vec{r}/r^2$. Sep 29, 2019 at 4:45

## 2 Answers

Using $$m\mathbf{a}=\mathbf{F}$$ and Newton's Law of Gravitation for $$\mathbf{F}$$, the two equations of motion are

$$m_1\ddot{\mathbf{r}}_1=-\frac{Gm_1m_2(\mathbf{r}_1-\mathbf{r}_2)}{|\mathbf{r}_1-\mathbf{r}_2|^3}\tag{1}$$

and

$$m_2\ddot{\mathbf{r}}_2=-\frac{Gm_1m_2(\mathbf{r}_2-\mathbf{r}_1)}{|\mathbf{r}_2-\mathbf{r}_1|^3}\tag{2}.$$

You can solve them numerically without worrying about using center-of-mass coordinates. If you choose initial conditions so that the initial momentum is zero, then the center of mass won't move.

If, for some reason, you want to use the center-of-mass vector

$$\mathbf{R}=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2}\tag{3}$$

and the separation vector

$$\mathbf{r}=\mathbf{r}_2-\mathbf{r}_1\tag{4}$$

instead of $$\mathbf{r}_1$$ and $$\mathbf{r}_2$$, then $$\mathbf{R}$$ satifies

$$\ddot{\mathbf{R}}=0\tag{5}$$

as follows by adding (1) and (2), and $$\mathbf{r}$$ satisifies

$$\ddot{\mathbf{r}}=-\frac{G(m_1+m_2)\mathbf{r}}{|\mathbf{r}|^3},\tag{6}$$

as follows by subtracting (1) divided by $$m_1$$ from (2) divided by $$m_2$$.

In terms of the reduced mass

$$\mu=\frac{m_1m_2}{m_1+m_2}\tag{7}$$

this can be written

$$\mu\ddot{\mathbf{r}}=-\frac{Gm_1m_2\mathbf{r}}{|\mathbf{r}|^3}\tag{8}$$

but when working numerically this is not any more useful than (6).

If you have numerically solved for $$\mathbf{r}$$, and have chosen a particular $$\mathbf{R}$$ for the position of the center of mass, then you can get back to $$\mathbf{r}_1$$ and $$\mathbf{r}_2$$ using

$$\mathbf{r}_1=\mathbf{R}-\frac{m_2}{m_1+m_2}\mathbf{r}\tag{9}$$

and

$$\mathbf{r}_2=\mathbf{R}+\frac{m_1}{m_1+m_2}\mathbf{r}\tag{10}$$

which follow from (3) and (4).

• thats for the correction in the unit vector rotatition. In Eqn(8) I thought that since we are defining a new force on the reduced mass the gravitational mass should change according to that. I am new at this topic and in the (1) and (2) the left side has $r_1$ but right has $r_2$ and $r_1$, so It seemed harders for me to solve it in that way. I said numerically but for now I have no idea how to do it. Sep 29, 2019 at 5:47
• @Reign You can solve 2 body systems exactly, see Kepler orbit & the various linked articles, including Kepler's equation. Sep 29, 2019 at 8:24
• @Reign But if you want to integrate the equations of motion numerically you should use a symplectic integrator, otherwise your simulation won't conserve energy. I like synchronized Leapfrog, especially the higher-order versions developed by Yoshida, but Verlet is also popular. Sep 29, 2019 at 8:25

This answer is supplementary to the earlier answer by G.Smith

About numerical analysis: as an introduction, just to get your feet wet, you can try the method described by Micheal Fowler, doing numerical analysis with a spreadsheet program.

You create columns with the values that you want to track (position, velocity). You fill the cells with computation so that the values in each row are used to calculate the values of the row below it.

Incidentally, in the case of gravitational 2-body motion you can use the symmetry in the problem to make the calculation more efficient.

The symmetry is most obvious when the two objects have equal mass. Then the motion of $$m_2$$ will always be symmetric to the motion of $$m_1$$. So you don't have to calculate the position of $$m_2$$ separately; you take the position of $$m_1$$, and you mirror that position with respect to the coordinates of the common center of mass.

When $$m_1$$ and $$m_2$$ have unequal mass you can still use that type of symmetry. For instance, if $$m_1$$ has twice as much mass as $$m_2$$ then the distance of $$m_2$$ to the common center of mass will always be twice as large as the distance of $$m_1$$ to the center of mass.

This mirroring approach has an additional advantage. If you would calculate the motions of $$m_1$$ and $$m_2$$ separately there is the possibility that accumulation of error makes the coordinates of the common center of mass drift away from the initial position, which would be unphysical.

The mirroring approach eliminates possibility of drift of the coordinates of the common center of mass.