newtonian-gravitynewtonian-mechanicsorbital-motion
Why is the center-of-mass of 2 bodies (which interact only via Newtonian gravity) located at a focus of each of the elliptical orbits?
I know that when there are no external forces, the center of mass moves at a constant speed, but that doesn't explain it.
The second (empty) focus is relevant in the theory of tides. In an elliptical orbit, the line joining the planet and the empty focus rotates at the same frequency as the mean motion of the planet; therefore, if spin rotation period is equal to the orbital period (the planet is locked in a synchronous rotation), the planet rotates with one face pointing to the empty focus.
Importantly, a tidal bulge will try to point to the massive object (the occupied focus) while the planet itself will be pointing to the empty focus, causing a "librational tide".
As Guillermo Angeris correctly pointed out, this is essentially a numerical roundoff problem, not a physical situation.
As a physical example, there are sungrazing comets that get very close to to Sun, yet they maintain their original elliptical (or hyperbolic) orbit, without the orbit precessing a full third of a circle as you seem to be seeing.
Computationally, there are a few interesting issues. As Kyle pointed out in a comment, many integration schemes are indeed unreliable in that roundoff error (which is always present in floating-point computations) can accumulate in a runaway feedback. Indeed I often advise using leapfrog methods over Euler (used by Box2D) or even Runge-Kutta (see for instance What is the correct way of integrating in astronomy simulations? over at the Computational Science Stackexchange).
However, I suspect your problem is even simpler, in the sense that even an unstable numerical scheme should work for one or two orbits. Given that everything is going wrong in just one pass, it seems that your timesteps are simply too large. A brief glimpse at the Box2D documentation suggests you don't change the timestep mid-simulation, so I presume you are just using a good value to simulate the whole process in reasonable time. The problem is that when gravitating bodies get close in their orbit, they move quickly, sometimes very quickly. The way the code works is it updates each object's position and velocity at each timestep, where the new velocity is determined by the force. As far as I can tell, this is done in line 206 of b2Island.cpp (v. 2.2.1):
v += h * (b->m_gravityScale * gravity + b->m_invMass * b->m_force);
Without looking at your code, I am guessing you simply calculate the gravitational force the body should feel at that moment, and have the simulation chug away. The problem is this moves the orbiting object in a straight line for the next timestep, and that straight line takes it too far away from the gravitating mass for that mass to properly curve its orbit into a closed ellipse. The quick schematic below shows the blue object moving to the tip of the red arrow, rather than staying on the path.
Physically, your timestep should be smaller than any timescale you encounter in the problem. Now for an orbit conserving angular momentum, the product of the orbiting body's mass, tangential velocity $v$, and distance from the other object $r$ should be constant: $v \sim 1/r$. At the same time, the acceleration $a$ it feels is given by Newton's law of gravity: $a \sim 1/r^2$. So one natural timescale in this problem is $$ t \sim \frac{v}{a} \sim \frac{1/r}{1/r^2} \sim r $$ (omitting dimensional constants), which goes to show that if your timescale is just barely small enough and then you tweak the orbit so as to half the periapsis distance (distance of closest approach), then you would expect to need timesteps at least twice as small in order to preserve the integrity of the simulation.
Best Answer
I'm going to assume that Omer is specifically asking why the centre of mass is at the focus (well, one of the foci) of the orbits. Omer, if this isn't what you meant please ignore what follows because it's completely irrelevant!
If you have a body moving in a central field (i.e. the force is always pointing towards the centre), and the field is inversely proportional to the square of the distance from body to the centre, then the orbit is an ellipse with the centre at one of the foci. For now let's just assume this and we can come back to prove it later.
So if we can show that both of the bodies feel a central inverse square force, with the COM at the centre, this guarantees the orbits will be ellipses with a focus at the COM. Given that the force is due to the two bodies attracting each other, and that both bodies are orbiting around, it may seem a bit odd that each body just feels a central inverse square force, but actually this is easy to show.
The picture shows the two masses and the COM. I haven't shown the velocities because it doesn't matter what they are. For now let's just consider $m_1$ and calculate the force on it. By Newton's law this is simply:
$$ F_1 = \frac{Gm_1m_2}{(r_1 + r_2)^2}$$
First is this force central? We know the centre of mass doesn't move. For two bodies this seems obvious to me, but in any case dmckee proved it in his answer. If the COM doesn't move it must lie on the line joining the two mases, otherwise there'd be a net force on it. So the force $F_1$ must always point towards the COM i.e. the force is central.
Second is this an inverse square law force i.e. is $F_1 \propto 1/r_1^2$? Well the definition of the centre of mass is that:
$$m_1r_1 = m_2r_2 $$
or
$$ r_2 = r_1 \frac{m_1}{m_2} $$
If we substitute for $r_2$ in the expression for $F_1$ we get:
$$ F_1 = \frac{Gm_1m_2}{(r_1 + r_1(m_1/m_2))^2}$$
or with a quick rearrangement:
$$ F_1 = \frac{1}{(1 + m_1/m_2)^2} \frac{Gm_1m_2}{r_1^2}$$
and this shows that $F_1$ is inversely proportional to $r_1^2$. I won't work through it, but it should be fairly obvious that exactly the same reasoning applies to $F_2$ so:
$$ F_2 = \frac{1}{(1 + m_2/m_1)^2} \frac{Gm_1m_2}{r_2^2}$$
This is the key result. Even though the two bodies are whizzing around each other, each body just behaves as if it were in a static gravity field, but the strength of the field is reduced by a factor of $(1 + m_1/m_2)^2$ for $m_1$ or $(1 + m_2/m_1)^2$ for $m_2$. This applies to all two body systems, even such unequal ones as the Sun and the Earth (ignoring perturbations from Jupiter etc).
I did start by assuming that a body in a central gravity field orbits in an ellipse with the foci at the centre, but I'm going to wimp out of proving this since it would double the length of this answer and you'd all go to sleep. The proof is easily Googled.
NB this only applies to two body systems. For three or more body systems the orbits are generally not ellipses with the centre of mass at the focus.