Creating a Gravitational 2-body simulation

Question

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 ?

Answer 1

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

Answer 2

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.