1
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – G. Smith
    Feb 28, 2020 at 1:08
  • $\begingroup$ 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/… $\endgroup$ Feb 28, 2020 at 1:34
  • $\begingroup$ 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? $\endgroup$ Feb 28, 2020 at 10:31

2 Answers 2

4
$\begingroup$

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

enter image description here

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.

$\endgroup$
7
  • 1
    $\begingroup$ Thank you soo much G.Smith! This makes complete sense! $\endgroup$ Feb 29, 2020 at 23:20
  • $\begingroup$ 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 $\endgroup$ Mar 1, 2020 at 4:28
  • $\begingroup$ Hi, Can you run me through your calculation? $\endgroup$ Mar 2, 2020 at 0:14
  • $\begingroup$ @Warrenmovic If you want to get Mitul’s attention, address him with @ like I here addressed you. $\endgroup$
    – G. Smith
    Mar 2, 2020 at 0:17
  • $\begingroup$ Thanks @G.Smith . $\endgroup$ Mar 2, 2020 at 10:36
1
$\begingroup$

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.