# N-body problem in Python

by epr2008
Tags: nbody, python
 P: 44 I was wondering if anyone could by any chance give me some help. I used Euler's method and the program runs, however, it doesn't give the correct solution. I honestly have no idea whats happening in the program, the problem has to be in the way I implemented Euler's Method. But I can't seem to find any errors. Then again, I really just started learning how to program on the fly about a month ago, and I am by no means any good at it, nor have I had the time to read any literature on programming, so my code may be a little bit cumbersome. I also, have no clue how to properly debug, I just kind of use techniques that are intuitive... Anyway here it is This is the class I created for vectors: ''' Created on Nov 19, 2012 @author: Penn ''' import math ''' Created on Nov 19, 2012 @author: Penn ''' import math class Vector: ''' classdocs ''' def __init__(self, data = []): ''' Constructor ''' self.data = data def __repr__(self): return repr(self.data) def __add__(self, other): data = [] for i in range(len(self.data)): data.append(self.data[i] + other.data[i]) return Vector(data) def __sub__(self, other): data = [] for i in range(len(self.data)): data.append(self.data[i] - other.data[i]) return Vector(data) def __mul__(self, a): data = [] for i in range(len(self.data)): data.append(a*self.data[i]) return Vector(data) def dot(self, other): dot = 0 for i in range(len(self.data)): dot += self.data[i]*other.data[i] return dot def mag(self): dot = Vector.dot(self, self) mag = math.pow(dot,1.0/2.0) return mag def unit(self): mag = Vector.mag(self) unit = Vector.__mul__(self, 1.0/mag) return unit def angle(self, other): dot = Vector.dot(self, other) A = Vector.mag(self) B = Vector.mag(other) theta = math.acos(dot / (A*B)) return theta def crossproduct(self,other): data = [] if (len(self.data) == 3): cross_x = self.data[1]*other.data[2] - self.data[2]*other.data[1] data.append(cross_x) cross_y = self.data[2]*other.data[0] - self.data[0]*other.data[2] data.append(cross_y) cross_z = self.data[0]*other.data[1] - self.data[1]*other.data[0] data.append(cross_z) return Vector(data) This is my class for Body objects: ''' Created on Nov 20, 2012 @author: Penn ''' from Vector import Vector class Body(Vector): ''' classdocs ''' def __init__(self, m, r=Vector(), v=Vector(), a=Vector()): ''' Constructor ''' self.mass = m self.position = r self.velocity = v self.acceleration = a self.positions = [] self.velocities = [] self.accelerations = [] def current_state(self): print 'Current Position: ', self.position print 'Current Velocity: ', self.velocity print 'Current Acceleration: ', self.acceleration print 'Mass: ', self.mass def update_acceleration(self, other): G = 6.6742 * pow(10,-11) position = Vector() self.acceleration = ((other.position - self.position) * (G * other.mass)) * (1.0/pow(Vector.mag(other.position - self.position),3.0)) return self.acceleration def update_velocity(self, h): self.velocity = self.velocity + (self.acceleration * h) return self.velocity def update_position(self, h): self.position = self.position + (self.velocity * h) return self.position def update_positions(self, h): self.positions.append(Body.update_position(self, h)) return self.positions def update_velocities(self, h): self.velocities.append(Body.update_velocity(self, h)) return self.velocities def update_accelerations(self, other): self.accelerations.append(Body.update_acceleration(self, other)) return self.accelerations And finally here is my main program. I set the number of bodies equal to 2 and just plugged in simple initial values to see if it worked: ''' Created on Dec 3, 2012 @author: Penn ''' #import tkSimpleDialog as tksd from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt import math from Vector import Vector from Body import Body if __name__ == '__main__': pass SOLAR_MASS = 4 * math.pi * math.pi DAYS_PER_YEAR = 365.24 Bodies = [] initialpositions = [Vector([0.0,0.0,0.0]),Vector([4.84143144246472090e+00, -1.16032004402742839e+00, -1.03622044471123109e-01])] initialvelocities = [Vector([0.0,0.0,0.0]),Vector([1.66007664274403694e-03 * DAYS_PER_YEAR, 7.69901118419740425e-03 * DAYS_PER_YEAR, -6.90460016972063023e-05 * DAYS_PER_YEAR])] initialaccelerations = [Vector([]),Vector([])] masses = [SOLAR_MASS, 9.54791938424326609e-04 * SOLAR_MASS] G = 6.6742 * pow(10,-11) #Get Number of bodies and time step from user '''root = tksd.Tk() N = tksd.askinteger("N-body Problem", "Enter the preferred number of bodies") t_ini = tksd.askfloat("N-Body Problem", "Enter the initial time:") t_fin = tksd.askfloat("N-Body Problem", "Enter the final time:") steps = tksd.askinteger("N-body Problem", "Enter the number of steps:") root.mainloop()''' N = 2 t_ini = 0.0 t_fin = 1000.0 steps = 100000 #Initialize each body and create an array of bodies for organizational purposes for i in range(N): '''Bodies.append(Body(masses[i], initialpositions[i], initialvelocities[i], initialaccelerations[i])) Body.__current_state__(Bodies[i])''' initialaccelerations[i] = Vector([0.0,0.0,0.0]) for k in range(N): if (k != i): initialaccelerations[i] = initialaccelerations[i] + (((initialpositions[k] - initialpositions[i]) * (G * masses[k])) * (1.0/pow(Vector.mag(initialpositions[k] - initialpositions[i]),3.0))) Bodies.append(Body(masses[i], initialpositions[i], initialvelocities[i], initialaccelerations[i])) #Body.current_state(Bodies[i]) #Calculate time step h = abs(t_fin - t_ini)/float(steps) for i in range(N): Bodies[i].positions.append(Bodies[i].position) Bodies[i].velocities.append(Bodies[i].velocity) Bodies[i].accelerations.append(Bodies[i].acceleration) for n in range(steps): for i in range(N): Bodies[i].velocity = Body.update_velocity(Bodies[i], h) Bodies[i].velocities.append(Bodies[i].velocity) Bodies[i].position = Body.update_position(Bodies[i], h) Bodies[i].positions.append(Bodies[i].position) Bodies[i].acceleration = Vector([0.0,0.0,0.0]) for k in range(N): if(k != i): Bodies[i].acceleration = Bodies[i].acceleration + Body.update_acceleration(Bodies[i], Bodies[k]) Bodies[i].accelerations.append(Bodies[i].acceleration) #Body.current_state(Bodies[i]) fig = plt.figure() ax = Axes3D(fig) for i in range(N): xlist = [] ylist = [] zlist = [] for n in range(steps+1): x = Vector.dot(Bodies[i].positions[n], Vector([1,0,0])) xlist.append(x) y = Vector.dot(Bodies[i].positions[n], Vector([0,1,0])) ylist.append(y) z = Vector.dot(Bodies[i].positions[n], Vector([0,0,1])) zlist.append(z) ax.plot(xlist,ylist,zlist) plt.show()
 P: 824 don't have time for much...just browsing for easy ones so, just one thing...are you sure you need to implement your own Vector class...look into numpy, I am sure you can simply use arrays for what you are doing and handle the operations that way, etc.
 Mentor P: 13,654 You have an error in Vector crossproduct(self,other). Look at the computation of the z component.
Mentor
P: 13,654

## N-body problem in Python

You also have (at least) a couple of errors in your calculation of the accelerations. One problem is that range(N-1); this should berange(N).

The other problem that I saw is that is your calculation of the acceleration of body i toward body k is wrong. You are using
$$a_{i,k} = \frac{Gm_k}{||\vec x_k-\vec x_i||}(\vec x_k - \vec x_i)$$
You should be using
$$a_{i,k} = \frac{Gm_k}{||\vec x_k-\vec x_i||^3}(\vec x_k - \vec x_i)$$
 P: 44 Thanks DH, caught this earlier but I didn't have time to post the updated code. It still does not change the result though, the planar components of the satelite body are still moving away from the origin, linearly. I tried it with actual values for jupiter and the sun, and still no change. I'll post the updated code now. And thanks also for catching my error in the cross product. No I probably did not need to write a class for Vectors. But I like doing things myself. I want to create a plotting tool of my own, but I am not well versed enough in python yet to even know where to start. And btw, I was wondering if anyone knew of a way to customize the syntax used to call methods of a class. I know that for "special methods" such as __add__, __append__, etc, are associated with special keys to call them, namely +, append() respectively. Is there away you can define syntax like this associated with user defined methods?
 P: 44 Also, I was wondering if my coding technique was convoluted or relatively clean for someone who just began learning? I usually go for organization as opposed to efficiency, then again I have realized that sometimes they go hand in hand.
 Mentor P: 13,654 A few last points. Please don't take these comments too negatively. As the saying goes in code review, you need to check your ego in at the door.You are using what is called symplectic Euler integration. About the only thing that is worse than this is Euler's method itself. This is the starting point rather than the ending point of numerical integration. As such, you should have your code structured in a way such that you can easily swap in a better integration technique. Your code is not structured in this way. Everything is too intertwined. Using G and mass is a sign of an amateurish simulation. It's far better to use the product μ=G*m, the standard gravitational parameter, for the body. The reason is precision. Physicists know G to only four decimal places, but astronomers can observe μ to much greater precision. The value for the Sun's standard gravitational parameter, for example, is known to over ten decimal places. Test, test, and then test some more. Python, more than any other language, is all about test driven development. I should not have found those basic errors that I found. Your testing should have ferreted them out.
 P: 44 No, not at all, I definitely appreciate the criticism. That was what I was hoping for. Like I said I just started learning how to program about a month ago. I have found that it's a lot more difficult to catch errors than with using a pencil and paper. Do you have any recommendations for learning to write a less intertwined code?
 P: 44 I'm guessing you are referring all of the for loops in my main module btw and the fact that I explicitly wrote the method for Euler's method in my class module. I was planning on cleaning it up once I got it working. I have been working on thinking more logically when programming. Thanks for the input!
Mentor
P: 13,654
 Quote by epr2008 Do you have any recommendations for learning to write a less intertwined code?
Practice. Lots and lots of practice. I've been practicing for 40+ years and I am still quite adept at writing convoluted, intertwined code.

Also its important to realize that there is no such thing as perfection, and even if there was, you wouldn't recognize that perfection even if it smacked you in the face. Write old-style and you'll have functions that are hundreds of lines long. These functions are so fragile and so internally intertwined that no one dare touch them. Write too new-style and you'll have hundreds of one liners that are so externally intertwined, forming an incredibly convoluted call tree that no one will touch this, either.

One of the most important skills you can learn is the ability to recognize your own garbage and be willing to throw that garbage out. (Aside: Managers probably need to learn this even more than programmers do.) There's nothing per se wrong with the way you wrote your code. You'll find out what's wrong when you try to plug in a more advanced integration technique. I know what's wrong with your code because I've dealt with numerical integration for far too long.

You don't necessarily want to anticipate ways you can do things better. There's a pejorative name for doing that: Premature optimization. AKA "if it ain't broke, don't fix it." The time to untangle your code is when you need/want to go beyond the simple scheme that does work. General rule: When you do do that, you'll find that your code is essentially broken; you need to refactor. This happens all the time. People generally aren't that good at anticipating what enhancements will need to be made to the code. The need for the enhancements they do anticipate never arises, and the things they did not anticipate oftentimes break the architecture.

The software community is finally seeing the light in this regard. The modern view is that if it ain't broke, don't fix it, and if it is broke, be willing to throw the old stuff out. Don't just patch over it. Refactor instead.