Lorenz Attractor in Processing

I coded up the Lorenz system instead of doing the rest of my work. By moving the mouse, the user can rotate the system, and then by using the arrow keys he or she can move the point about which the plot rotates. What's cool about this to me is that its solving the system of differential equations in real time, while you're rotating it. The color reflects the depth (z-axis) of the plot. I'd like to set it up such that one drags the mouse to rotate. Maybe once the semester is over.... For now, you're stuck with this GIF.

I had some trouble getting ArrayLists to work in Processing, as they feel very foreign to me (they're so Python) but I eventually figured it out. If you're interested, you can see the code for this project here.

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.

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.

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.

Two body problem in processing

Learning Processing has been a fun (and frustrating) experience for me. Coding in Java has me realizing how much Python spoils me. I keep getting NullPointerExceptions in Java because of how I deal with objects that I haven't initialized yet. I like that I have to declare the type for every variable. This organizes my thinking, and forces me to be more deliberate with my code. Processing is a powerful tool for visualization, but it's not as easy to use for doing math as Numpy. In Numpy adding two arrays together is trivial. So is performing operations. In Java, at least with regular arrays, you have to iterate through the array. Ugh. 

#Python
x = np.random.random((10,))
y = np.sin(x)
h = np.random.random((10,)) + x

//Java/Processing
float[] x = new float[10];
float[] y = new float[10]; 
float[] h = new float[10];
for (int i=0; i < 10; i++){
x[i] = random(10);
y[i] = sin(x[i]);  
h[i] = random(10) + x[i];
}

Anyways, I set about solving the two body problem in Processing. Normally when solving the two body problem, you do a coordinate transform, which ultimately ends up turning the problem into a 1 body problem. I want to turn my code into a three or four body code, so I stuck with solving things in cartesian coordinates. You might be thinking that I'm crazy, but I don't mind the extra lines of code for the simplicity of my approach. Let's dive in. The first thing to do is to set up the Lagrangian for our system. The idea of the two body problem is solve for the equations of motion of two particles who are mutually gravitationally attracted to each other. The Lagrangian for the system is as follows. $$ L = T-V\\ L = \frac{m_1}{2}(\dot{x_1}^2 + \dot{y_1}^2) + \frac{m_2}{2}(\dot{x_2}^2 + \dot{y_2}^2) - \frac{Gm_1m_2}{\sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2}} $$ By applying the Euler-Lagrange equations we can arrive at the following equation of motion for the four coordinates in our system. $$ \ddot{x_1} = -Gm_2(x_1-x_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{y_1} = -Gm_2(y_1-y_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{x_2} = Gm_1(x_1-x_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2}\\ \ddot{y_2} = Gm_1(y_1-y_2)[(x_1 - x_2)^2 + (y_1 - y_2)^2]^{-3/2} $$ Now we have to solve a system of nonlinear second order differential equations. We can simplify this problem a lot by turning each of the second order equations into a first order equation. I've talked about this in previous posts, but I'll go over it again briefly. Take, for instance, the first coordinate $x_1$. Let me create new variables $x_1^a$ and $x_1^b$ defined as follows. $$ x_1^a = x_1\\ x_1^b = \dot{x_1}\\ \dot{x_1}^a = \dot{x_1} = x_1^b\\ \dot{x_1}^b = \ddot{x_1} $$ Now I can rewrite my second order equation for $x_1$ as a system of two first order differential equations in terms of $x_1^a$ and $x_1^b$. Note that I've gone ahead and introduced the 'a' and 'b' notation for the other coordinates. $$ \dot{x_1}^a = x_1^b\\ \dot{x_1}^b = -Gm_2(x_1^a-x_2^a)[(x_1^a - x_2^a)^2 + (y_1^a - y_2^a)^2]^{-3/2} $$ With these equations in hand I can use any first order ODE method to step my equations through time. I used Runge Kutta in my code. I coded all this up in Processing, and I made a .gif for your enjoyment! Check it out below.

You might notice some decaying of the orbit. This is not physical, rather it represents a loss of stability in my ODE solver. This is even more pronounced when using a simple method like forward Euler. Also notice that I make no assumptions about the position or mass of the particle at the center of the frame -- I just defined it as having about 100x the mass of the second particle.

Two point charges in Processing

I've been messing around with Processing a little bit lately. In my last post, I animated the double pendulum, and it became pretty clear that I had messed up the math somewhere in deriving the equations of motion for the system. This time, I wanted to show the vector field around a point charge, or a pair of point charges. Here's a gif of the processing sketch I coded up:

Double Pendulum in Processing

I'm taking a class called Introduction to Interactive Media this semester. The second half of the semester concentrates mostly on making cool stuff in Processing. As much as I like Python, it has a few shortcomings. Making animations is one. I thought it would be cool to take some of code I used for solving differential equations (in Python) and plug it into Processing to see what happens. My first experiment involves animating the double pendulum. From this animation, I think it's pretty clear that I messed up the math somewhere, as it has some rather weird behavior. Check it out:

I'll have to re solve the equations of motion of this system to be sure that I have it right. Otherwise, I think part of the issue might be the method I'm using to solve the system of differential equations. I haven't (yet) incorporated a linear algebra package into my Processing code, so I'm stuck using forward Euler method. Alternatively, I think I could use another explicit method that allows for better stability, like Runge-Kutta. Another thing I'd like to do is to make the speed at which the animation plays independent of the step sized used in my solver.

I see Processing as a really cool tool for Physics visualizations. Stay tuned; I plan on doing a few more of these simple animations over my Fall break.