I've been adding some fun features to my N-body code. Now you can grab the masses with your mouse to adjust their initial positions. Additionally, you can change the mass of the particles, which is reflected in their size. Here's a little sample of some of the new stuff. I've got it set up for four bodies right now. Adding more bodies is not really a big deal now, which makes me super happy. I'm planning on doing a little write up of the Physics of this code in the next couple of days. I also want to talk about how I should have been using Hamiltonian mechanics the whole time, as it results in a system of first order ODEs instead of a single second order ODE. For now, you're stuck with this .gif
L = \frac{1}{2}\sum^{N}_{i=0} m_i(\dot{x_i}^2 + \dot{y_i}^2) - \frac{1}{2}\sum^{N}_{i=0}\sum^{N}_{j=0,j\ne i}\frac{Gm_im_j}{\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}} $$
Notice that we divide by a half in the potential term because we're double counting everything in the sums. Solving for the equations of motion is actually quite easy, as we just plug everything into the Euler-Lagrange equation. For the $ith$ particle, the $x$ equation is as follows:
$$ m_i\ddot{x_i} = \sum^{N}_{j=0,j\ne i}\frac{G m_i m_j (x_i - x_j)}{[(x_i - x_j)^2 + (y_i - y_j)^2]^{3/2}} $$
Here, I think it becomes clear why this becomes computationally difficult for large N. As N increases, we slap on another term to each of our equations of motion, in addition to adding another two equations of motion (for each of the coordinates of the new particle). My java implementation starts to noticeably slow down after about seven particles.
For those of you who have been wondering about the physics of this simulation, I'll do a little write up. This is a relatively simple extension from the two or three body problem that you might find at the beginning of an introductory text on E&M or mechanics. Using the Lagrangian is definitely overkill here, but I really prefer using scalar quantities with generalized coordinates instead of messing around with vectors. Anyways, we start with our Lagrangian, $L$, in this case defined as follows.
$$
L = T-V \\L = \frac{1}{2}\sum^{N}_{i=0} m_i(\dot{x_i}^2 + \dot{y_i}^2) - \frac{1}{2}\sum^{N}_{i=0}\sum^{N}_{j=0,j\ne i}\frac{Gm_im_j}{\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}} $$
Notice that we divide by a half in the potential term because we're double counting everything in the sums. Solving for the equations of motion is actually quite easy, as we just plug everything into the Euler-Lagrange equation. For the $ith$ particle, the $x$ equation is as follows:
$$ m_i\ddot{x_i} = \sum^{N}_{j=0,j\ne i}\frac{G m_i m_j (x_i - x_j)}{[(x_i - x_j)^2 + (y_i - y_j)^2]^{3/2}} $$
Here, I think it becomes clear why this becomes computationally difficult for large N. As N increases, we slap on another term to each of our equations of motion, in addition to adding another two equations of motion (for each of the coordinates of the new particle). My java implementation starts to noticeably slow down after about seven particles.