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 '20 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$ – QuantumApple Feb 28 '20 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$ – Warrenmovic Feb 28 '20 at 10:31
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$ – Warrenmovic Feb 29 '20 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$ – Mitul Agrawal Mar 1 '20 at 4:28
  • $\begingroup$ Hi, Can you run me through your calculation? $\endgroup$ – Warrenmovic Mar 2 '20 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 '20 at 0:17
  • $\begingroup$ Thanks @G.Smith . $\endgroup$ – Warrenmovic Mar 2 '20 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, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.