Sunday, November 8, 2015

N body update

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

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.

Tuesday, November 3, 2015

N body in Processing

This was bound to happen. I coded up a super basic n-body simulation in Processing. It works by just iteratively adding terms to the equations of motion for each particle. Below is a .gif for 7 (7!) particles. The central one is more massive than the others, but there is (at least in principle) still interaction between the other particles. Notice how slow it goes! That means my computer is really taking some time to crunch those numbers! There's probably a good reason why we don't use this naive method for real N-body simulations.

Sunday, November 1, 2015

Three body problem in Processing

I've been messing around with doing some very basic physics simulations in Processing. Last week I did the two body problem, which describes the motion of two particles under mutual gravitation. This week I extended the code a little to do the three body problem -- three particles under mutual gravitation. Coding this up was pretty easy, but finding initial conditions that resulted in something viewable was quite hard! I want to build in some interactivity such that a user can dynamically adjust initial conditions, but I've been too lazy to get around to doing that. Maybe next time I'll have done that! For now, however, you're stuck with a .gif I made:



From this little animation, you can see how chaotic this system is! You can set up a situation in which you get stable orbits. In the simulation below I set the mass of the central particle to be about 200 times that of the outer one. The one in the middle orbit is about 100 times smaller than that of the outer one.



Here you can see that the central mass is moving a little. This is because it is orbiting around the center of mass of the outer/central mass sub system. For those who are interested, I changed from float to double arrays in this code and my previous two body code. The stability of the solution increased in the long term as a result. As always, using an implicit ODE method would probably be the best way to ensure long term stability of my simulation.

UPDATE: As promised, I added some sliders to my code so that users can adjust the initial parameters of the system before letting it loose. Check out the gif below.