solar system - Three-body problem - sun velocity around barycenter

I am really, really stuck! I am trying to write Java code to model the solar system, or at least currently just earth, the sun and Mercury. I just wanted to check that my actual physics equations are right, and my initial conditions!

I've set the sun to start at the origin, and the earth's position is $$x = 1.496 \times 10^{11}m$$, and Mercury is at $$x = 5.971 \times 10^{10}m$$

I want everything to be done in the centre-of-mass frame, so I then calculate the position of the centre of mass and subtract it from the initial positions, e.g. so the sun is at x - centreOfMass. The velocity of each body is given by $$v = \frac{2\pi r}{T}$$. This is probably where I go wrong. I know the orbital period of Mercury and earth around the sun, so I can calculate their velocities around the centre of mass. But the sun will also move, even if only a tiny, tiny bit, and I don't have its orbital period around the centre of mass! So how can I work out its velocity?

Assuming I have the sun's velocity around the centre of mass, I then work out the velocity of the centre of mass. Because I want to stay in the centre-of-mass frame, I then subtract this velocity from the velocities of the planets.

My code is then written so that the position is recalculated, given a time step, the centre-of-mass position is recalculated and subtracted from the new position. A new velocity for each planet is calculated, then the centre of mass velocity is recalculated and again subtracted from each velocity.

Is this the correct procedure? Do I not need to recalculate the centre-of-mass position/velocity after each time step? I've made a big mistake somewhere; the planets are very unstable and either move in a line or spiral inwards!

Answer

Is this the correct procedure? Do I not need to recalculate the centre-of-mass position/velocity after each time step? I've made a big mistake somewhere; the planets are very unstable and either move in a line or spiral inwards!

This is probably down to the stability of your numerical system and not the center of mass of the Sun. (Which by the way is a non-issue: the position of the Sun can be calculated uniquely from the positions of the Earth and Mercury, simply by placing the center of mass at $(0, 0)$ in your coordinates and then figuring out where the Sun now has to be for this point to be the center of mass.) This is a notoriously difficult problem to get good numerical solutions for.

If you want a good reference on how to transform this so that key parameters like total angular momentum and total energy are conserved, see this PDF from UAlberta for instance, or search for "N-body conserves energy angular momentum" on Google or some such; it's one of those fields where there's a ton of literature on it and it hasn't conclusively recommended "this is the best way!" for any one approach.

Three things that might help: so, your direct simulation should look like: $$\begin{array}{cl} \vec r_i^+ ~=& \vec r_i + \vec v_i ~\Delta t\\ \vec v_i^+ ~=& \vec v_i - G \Delta t \sum_j m_j (\vec r_i - \vec r_j) / |\vec r_i - \vec r_j|^3\end{array},$$ where the superscript $+$ signifies "for the next time step" and unadorned numbers are for the current time step.

1. 
   

   One huge problem newcomers to simulations have is that they start updating some array of positions incrementally. This always causes bugs. You should not be updating the forces on the Earth using the position of Mercury from the next time step, just because you happened to process Mercury first. Do not overwrite their present positions, instead, have two copies: their present positions and their future positions at the next time step. First re-calculate all of these these future positions; then copy those all at once to the current positions, so that for each future-calculation you know it's based on the current time.

   
   



2. 
   

   It helps in these circumstances to actually interpret $\vec v_i$ as being $\vec v_i(t - \Delta t/2)$ rather than to $\vec v_i(t)$, because that "balances" the derivatives a little bit more. Similarly $\vec v_i^+ = \vec v_i(t + \Delta t/2).$ This has no effect on the second calculation above, but it changes the first to use $\vec v_i^+$ instead of $\vec v_i$.

   
   


3. 
   

   Consider cheating. Consider after each timestep nudging your objects a little bit so as to preserve angular momentum and total energy. There's probably a formal way to do this with Lagrange multipliers and such, but if you pretend that the center of the Sun is the center of the solar system (which it roughly will be, until you get to Jupiter), you should be able to just do it ad-hoc on the assumption that these planets are not stealing each other's energy too much.