Solving Keplerian two-body problem numerically?

Question

I'm trying to solve the two-body problem numerically, setting up G, m1 and m2 =1.0. The masses are placed at the positions -10 and 10 respectively along the x-axis and gave them both 0 on the y-axis. I am having some real issue with the initial conditions fitting for a keplerian orbit with a given eccentricity (e=0.9)! Also how would the initial condition change for different mass ratios?

Answer 1

If we give the two masses initial velocities in the $\hat{y}$ and $-\hat{y}$ directions, perpendicular to their separation, we'll be starting at either periapsis or apoapsis.

The vis-viva equation for a binary system says that

$$v^2=G(m_1+m_2)\left(\frac{2}{r}-\frac{1}{a}\right)\tag1.$$

Here $r$ is the distance between the two masses, which is the magnitude of the separation vector

$$\mathbf{r}=\mathbf{r}_1-\mathbf{r}_2\tag2;$$

note that this is not the distance of either from their barycenter.

The speed $v$ is the magnitude of the relative velocity vector

$$\mathbf{v}=\frac{d\mathbf{r}}{dt}=\mathbf{v}_1-\mathbf{v}_2\tag3.$$

The semimajor axis $a$ is for the ellipse formed by $\mathbf{r}$, not the smaller ellipses formed by $\mathbf{r}_1$ or $\mathbf{r}_2$.

The ellipse formed by $\mathbf{r}$ has the form

$$r=\frac{a(1-e^2)}{1+e\cos\theta}\tag4$$

where $\theta=0$ at periapsis and $e$ is the eccentricity.

The value of $r$ at apoapsis is

$$r_a=a(1+e)\tag5.$$

So the relative speed at apoapsis is

$$\begin{align} v_a&=\sqrt{G(m_1+m_2)\left(\frac{2}{r_a}-\frac{1}{a}\right)}\\ &=\sqrt{G(m_1+m_2)\left(\frac{2}{r_a}-\frac{1+e}{r_a}\right)}\\ &=\sqrt{\frac{G(m_1+m_2)}{r_a}(1-e)}\\ \end{align}\tag6.$$

Thus, given your initial data, we can find the initial relative speed.

To determine the initial speed of each mass, we use the fact that

$$\mathbf{r_1}=\frac{m_2}{m_1+m_2}\mathbf{r}\tag7$$

and

$$\mathbf{r_2}=-\frac{m_1}{m_1+m_2}\mathbf{r}\tag8$$

which follow from (2) plus the condition that the center of mass is at the origin:

$$\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2}=0\tag9.$$

We thus have

$$\mathbf{v_1}=\frac{m_2}{m_1+m_2}\mathbf{v}\tag{10}$$

and

$$\mathbf{v_2}=-\frac{m_1}{m_1+m_2}\mathbf{v}\tag{11}$$

Putting in your numbers $G=m_1=m_2=1$, $r_a=(10)-(-10)=20$, and $e=9/10$, we find the initial speeds should be

$$v_1=v_2=\frac{1}{1+1}\sqrt{\frac{(1)(1+1)}{20}\left(1-\frac{9}{10}\right)}=\frac{1}{20}\tag{12}.$$

Trying this in Mathematica using

data = NBodySimulation[ "InverseSquare", {<|"Mass" -> 1, "Position" -> {10, 0}, "Velocity" -> {0, 1/20}|>, <|"Mass" -> 1, "Position" -> {-10, 0}, "Velocity" -> {0, -1/20}|>}, 400]

ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 400}]

one gets

As a check, the separation at periapsis should be

$$r_p=a(1-e)=r_a\frac{1-e}{1+e}=20\frac{1-\frac{9}{10}}{1+\frac{9}{10}}=\frac{20}{19}\approx 1.05\tag{13}$$

and it seems to be.

Answer 2

You can use Mathematica to solve the numerical solution of this problem:

data  = NBodySimulation[
  "InverseSquare", {<|"Mass" -> 1, "Position" -> {0, 0}, 
    "Velocity" -> {0, .5}|>,
   <|"Mass" -> 1, "Position" -> {1, 1}, "Velocity" -> {0, -.5}|>}, 4]
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 4}]