How to Solve N-body Problem with Python using Euler's Method?

[epr2008]

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 what's 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()

 

[gsal]

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.

 

[D H]

You have an error in Vector crossproduct(self,other). Look at the computation of the z component.

 

[D H]

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)$$

 

[epr2008]

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?

 

[epr2008]

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.

 

[D H]

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.

 

[epr2008]

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?

 

[epr2008]

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!

 

[D H]

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.

 

[D H]

My last post was about the software development process. This one is about numerical integration. You are, I think, running into typical numerical integration problems.

I see four basic classes of errors that arise in numerically integrating an ordinary differential equation. I'll look at what happens as one increases the step size from a very small number to a very large number.

Regime of frozen integration.
Suppose you are using double precision arithmetic (e.g. double in C and C++, double precision in Fortran). If you calculate 1+10^-100 or even 1+10^-16 on your computer, you are going to get one as a result. Exactly one. That additional 10^-100 or 10^-16 just vanishes. Those double precision numbers are but a poor man's substitute for the real numbers. When the step size is too small, x+v*Δt won't change the position one iota. The numerical integration remains frozen at the initial value in this regime of step sizes that are way, way too small.

Regime of machine errors.
Eventually, increasing the step size will result in values that do change. There's still a problem here. Suppose that x has a value of 1 and v*Δt has a value of 1.01×10^-15. When you compute x_new=x_old+v*Δt, that 0.01×10^-15 is still going to vanish. The number you are adding effectively has one decimal digit of precision. These truncation errors add up, and in this extreme case they add up rather quickly. Increase the step size a bit and the problem with truncation error decreases a bit. In this regime, increasing the step size tends to improve accuracy.

Regime of technique errors.
Increase the step size even more and you'll find that accuracy now starts dropping as the step size increases. You are now in the regime where the errors inherent in the integration technique itself dominate over the errors that result from using finite precision arithmetic. This transition from machine errors to technique errors varies dramatically with different integration techniques. This transition also depends one the specific problem that one is trying to solve. For a simple circular orbit, this transition occurs at the ridiculously small step size of about 10^-12 of an orbit for symplectic Euler (Euler is even worse). For some (very) advanced integration techniques, the transition between numerical errors and technique errors doesn't occur unless the step size becomes multiple orbits long.

Regime of catastrophic errors.
Many numerical integration techniques have a stability boundary. Step outside this boundary by making the step size way too large and the technique becomes numerically unstable. When this happens, things go bad, very very bad, very very quickly. You don't want to cross this stability boundary. Numbers tend to go to infinity when you do.The very best one can do is to try to find that sweet spot where machine error and technique error are more or less the same size. Some adaptive integration techniques try to find this sweet spot automagically. Sometimes the magic works, sometimes it doesn't. Sans using such techniques, you need to be aware of these four different ways in which errors arise using numerical integration. The burden of hitting that sweet spot is all upon you.