Solving the two body problem numerically

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!

3 Answers

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.

• Suppose I set $T=1$ and $r$ is then 5, then the velocity should be 10$\pi$ in the $y$ direction. This inital condition won't provide a circular movement. What am I missing? – Talco Aug 4 '15 at 15:43
• Unstable integration scheme would be my first guess. What integration method are you using? – Kyle Kanos Aug 4 '15 at 15:49
• Euler's Method. But then if I set any velocity in the $y$ direction it should give me circular movement, with different periods which is not the case, nor should it be. I think there's a specific velocity, which yields the desired result. – Talco Aug 4 '15 at 15:58
• Well Euler's method isn't really that great for orbital dynamics (not a symplectic integrator). You'll want to use something like the Verlet integration. – Kyle Kanos Aug 4 '15 at 16:14
• Easiest way to check if your integrator is your (likely) problem is to turn the timestep way down (e.g. 10x smaller). If this meaningfully changes your orbits, you've got numerical rather than physics issues. – Kyle Oman Aug 4 '15 at 16:38

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.

• Covering all the points I was going to make when I read the question - right down to the need for "clever" integration! – Floris Aug 4 '15 at 15:33
• If I followed your explaination correctly, then using my initial placing, the center of mass is at the origin. Hence, $R=5$ for both of the masses. Solving for $v$ yields $v=1/\sqrt(5)$ which is the speed in the $y$ direction. Why can't I set this velocity as an initial condition? Thanks. – Talco Aug 4 '15 at 15:53
• If you want to simplify the setup, you can start one particle at rest at the origin and set up the other one for a circular orbit around it. But note that your "stationary" particle will start to move, and your system will have a total momentum, so the CoM will wander along a straight path. You can always subtract this motion later if you want (or from the ICs!). Verlet is a good integrator, but if you want to try another leapfrog is also good for this problem, and is the usual choice for beginners (and often for experts tbh), so you may find tutorials using it. – Kyle Oman Aug 4 '15 at 16:36
• @Talco No, you made a mistake. The radius of orbit is not the same as the separation distance. Yes, R=5, but the distance between the particles is 10. You can set this value as the initial speed of each particle, one in the +y direction and the other in the -y direction, initially. – Bill N Aug 5 '15 at 0:26

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.