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.(adsbygoogle = window.adsbygoogle || []).push({});

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:

This is my class for Body objects:Code (Text):

'''

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)

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:Code (Text):

'''

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

Code (Text):

'''

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

**Physics Forums - The Fusion of Science and Community**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# N-body problem in Python

Loading...

Similar Threads - body problem Python | Date |
---|---|

The mind-body problem for computers | May 2, 2017 |

Runge-Kutta methond for n body C++ | Mar 15, 2016 |

3D N-Body Problem (Solar System) C# | Jan 10, 2014 |

Two Body problem with RK4 in C++ | Dec 14, 2012 |

Two Body problem in python using RK4 | Nov 11, 2012 |

**Physics Forums - The Fusion of Science and Community**