# Solving Keplerian two-body problem numerically?

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?

• You can try different initial velocities and see what eccentricity they produce. Adjust them until you get 0.9. If you are doing a numerical simulation isn’t that more satisfying than having us work out a formula? – G. Smith Feb 28 '20 at 1:08
• If you want to go through the pain of finding the right initial velocity that will yield a given eccentricity, this answer might help you get started: physics.stackexchange.com/questions/522208/… – QuantumApple Feb 28 '20 at 1:34
• I have used the vis-viva equation v = sqrt( GM * (2/r -1/a)) at apocentre but that yields a incorrect velocity. Is there no simple equation to calculate the initial velocitites? – Warrenmovic Feb 28 '20 at 10:31

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.

• Thank you soo much G.Smith! This makes complete sense! – Warrenmovic Feb 29 '20 at 23:20
• Hi, I have a doubt......if we calculate the velocity required for circular orbit by gravitational force = centrifugal we get the same velocity as that of an elliptical orbit with e = 0.5 so in that case what will happen – Mitul Agrawal Mar 1 '20 at 4:28
• Hi, Can you run me through your calculation? – Warrenmovic Mar 2 '20 at 0:14
• @Warrenmovic If you want to get Mitul’s attention, address him with @ like I here addressed you. – G. Smith Mar 2 '20 at 0:17
• Thanks @G.Smith . – Warrenmovic Mar 2 '20 at 10:36

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