Solving the two body problem numerically

Question

I'm trying to solve the two body problem numerically, setting up $G$, $m1$ and $m2$ to be equal to 1. then I located each mass on positions -5 and 5 respectively along the $x$ axis and gave them both 0 on the $y$ axis. I'm having trouble finding the initial conditions fitting for circular movement. Please help!

Answer 1

For common center (barycenter) orbits, the velocities will be \begin{align} v_1&=\sqrt{\frac{Gm_2r_2}{\left(r_1+r_2\right)^2}}\\ v_2&=\sqrt{\frac{Gm_1r_1}{\left(r_1+r_2\right)^2}} \end{align} which, since $m_1=m_2$ and $r_1=r_2$, will be the same value, $v\approx0.22$ for your values of $G,\,m,\,r$. Since you've placed the two objects along the $x$ axis, then you need to give them this velocity in the $\pm y$ direction and 0 in the $x$ direction (the positive/negative value being associated with the direction of motion of the particle).

You can see why with the below diagram from the Wikipedia entry on circular motion,

At any one of the four cardinal points, the velocity is purely in the tangential direction. For example, consider the top point; it is at $(x,y)=(0,1)$ with a velocity of $(v_x,\,v_y)=(-v,0)$. You can do the same analysis with the other 3 points. So rather than trying to determine a random $(x,y)$ point, you can simply use the cardinal points to set the positions and velocities.

Answer 2

I assume that you want both particles to be free to move, each in a circle. That circle should be about the center of mass of the system, which will be at $$\frac{m_1\vec{r}_1+m_2\vec{r}_2}{m_1+m_2}.$$

The radius, $R$ of each circle will be the distance from the mass to the center-of-mass point. Calculate the Newtonian gravity force on a particle, set it equal to $mv^2/R$, where $R$ is the circle radius but the Newtonian gravity is calculated using the distance to the other particle. Solve for $v$. This is the speed perpendicular to the line separating the particles. The direction of this velocity needs to be updated every time step, so you'll have to do good vector algebra.

Setting the initial conditions is only the beginning of the problem. Consider using the Verlet velocity algorithm. A simple Euler's method will deteriorate quickly.

Answer 3

Your question is about finding the initial conditions to a circular orbit, which has already been answered very well. This answer instead tackles determining the initial conditions for more general orbits. This will allow you to simulate different orbit shapes.

Below are a few shapes of different orbits:

There are two parameters that govern the shape: the eccentricity, $\epsilon$, and the semi-latus rectum, $l$.

$\epsilon$ is how much the shape varies from a circular orbit. $\epsilon = 0$ is a circular orbit. For low values, the orbit is elliptical. Higher values yield open "orbits", where the body only passes by once.

$l$ is a particular length dimension perpendicular to the line of symmetry, governing the general size of the orbit. Doubling $l$ doubles the orbit size.

So you can define the shape of an orbit. For the general orbital set-up below:

You have two masses, $m_1$ and $m_2$, each with a position vector from the barycentre ($\vec r_1$, $\vec r_2$), a distance from the barycentre ($r_1$, $r_2$), a radial velocity ($v_{r1}$, $v_{r2}$) which is parallel to the line connecting the mass and the barycentre, and a tangential velocity ($v_{t1}$, $v_{t2}$) which is perpendicular to the line connecting the mass and the barycentre. Also note the angle between the apex of the orbit (sharpest bend) and the mass, $\theta$.

The two masses have orbit shapes with different semi-latus recti ($l_1$, $l_2$) but the same eccentricity.

You can determine the initial distance, and velocities from the shape of the desired orbit you wish to simulate, and the starting angle, $\theta$. Here are the relationships for mass 1:

$$r_1 = \frac{l_1}{1 + \epsilon{\cos \theta}}$$

$$v_{r1} = \sqrt{\frac{G \mu_2}{l_1}} \epsilon \sin{\theta}$$

$$v_{t1} = \sqrt{\frac{G \mu_2}{l_1}} (1+\epsilon \cos{\theta})$$

where $\mu_2$ is a mass factor equal to $\frac{m_2^3}{(m_1 + m_2)^2}$ which is approximately equal to $m_2$ for large $m_2$.

Then the initial condition for mass 2 can be given by:

$$r_2 = r_1 \frac{m_1}{m_2}$$

$$v_{r2} = v_{r1} \frac{m_1}{m_2}$$

$$v_{t2} = v_{t1} \frac{m_1}{m_2}$$

And that should do it! Note that putting $\epsilon = 0$ will get you the same results for the usual circular orbit.