- #1

- 44

- 0

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:

This is my class for Body objects:

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:

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:

Code:

```
'''
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)
```

Code:

```
'''
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:

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

Last edited: