I wrote python code to solve the N-body problem using the Euler method. The code works without problems and seems to give a reasonable answer (for example, if there are two particles, then the beginning of the movement towards each other). However, when I run this simulation on a large number of iterations, I see that the particles (say, I run them with two particles) pass with each other (I do not consider collisions) and continue to move in their directions for an indefinite period. This violates energy conservation, and therefore there should be a flaw in my code, but I cannot find it. Can anyone find it and explain my mistake.
Thank.
Thanks to @samgak for updating the particles twice. I have now fixed this, but the problem is still ongoing. I also reproduced the result when I run this simulation with two stationary particles (0,0) and (1,0) in 1 second increments and 100,000 iterations:
Particle with mass: 1 and position: [234.8268420043934, 0.0] and speed: [0.011249111128594091, 0.0]
Particle with mass: 1 and position: [-233.82684200439311, 0.0] and speed: [-0.011249111128594091, 0.0]
Also thanks to @ PM2Ring for pointing out some optimizations I could do and the dangers of using the Euler method.
the code:
import math
class Particle:
"""
Class to represent a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11
self.time_interval = 10**0
self.mass = mass
self.position = position
self.velocity = velocity
self.updated_position = position
self.updated_velocity = velocity
def __str__(self):
"""
String representation of particle
"""
return "Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity: " + str(self.velocity)
def get_mass(self):
"""
Returns the mass of the particle
"""
return self.mass
def get_position(self):
"""
returns the position of the particle
"""
return self.position
def get_velocity(self):
"""
returns the velocity of the particle
"""
return self.velocity
def get_updated_position(self):
"""
calculates the future position of the particle
"""
for i in range(len(self.position)):
self.updated_position[i] = self.updated_position[i] + self.time_interval*self.velocity[i]
def update_position(self):
"""
updates the position of the particle
"""
self.position = self.updated_position.copy()
def get_distance(self,other_particle):
"""
returns the distance between the particle and another given particle
"""
tot = 0
other = other_particle.get_position()
for i in range(len(self.position)):
tot += (self.position[i]-other[i])**2
return math.sqrt(tot)
def get_updated_velocity(self,other_particle):
"""
updates the future velocity of the particle due to the acceleration
by another particle
"""
distance_vector = []
other = other_particle.get_position()
for i in range(len(self.position)):
distance_vector.append(self.position[i]-other[i])
distance_squared = 0
for item in distance_vector:
distance_squared += item**2
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/(distance_squared)
for i in range(len(self.velocity)):
self.updated_velocity[i] = self.updated_velocity[i]+self.time_interval*force*(distance_vector[i])/(self.mass*(distance))
def update_velocity(self):
"""
updates the velocity of the particle
"""
self.velocity = self.updated_velocity.copy()
def update_particles(particle_list):
"""
updates the position of all the particles
"""
for i in range(len(particle_list)):
for j in range(i+1,len(particle_list)):
particle_list[i].get_updated_velocity(particle_list[j])
particle_list[j].get_updated_velocity(particle_list[i])
for i in range(len(particle_list)):
particle_list[i].update_velocity()
particle_list[i].get_updated_position()
for i in range(len(particle_list)):
particle_list[i].update_position()
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
for i in range(100000):
update_particles(partList)
for item in partList:
print(item)
----------------------------------------------- -------------------------------------------------- -------------------------------------------------- -------------------------------------------------- -------------------------------------------------- ----------------- :
Leapfrog, , , , ( , ). , , . , , , . , , . , , . , .
- , .
:
import math
import matplotlib.pyplot as plt
class Particle:
"""
Represents a single particle
"""
def __init__(self,mass,position,velocity):
"""
Initialize the particle
"""
self.G = 6.67408*10**-11
self.time_step = 10**2
self.mass = mass
self.dimensions = len(position)
self.position = position
self.velocity = velocity
self.acceleration = [0 for i in range(len(position))]
self.next_position = position
self.next_velocity = velocity
self.next_acceleration = [0 for i in range(len(position))]
def __str__(self):
"""
A string representation of the particle
"""
return "A Particle with mass: " + str(self.mass) + " and position: " + str(self.position) + " and velocity:" + str(self.velocity)
def get_mass(self):
return self.mass
def get_position(self):
return self.position
def get_velocity(self):
return self.velocity
def get_acceleration(self):
return self.acceleration
def get_next_position(self):
return self.next_position
def put_next_position(self):
for i in range(self.dimensions):
self.next_position[i] = self.position[i] + self.time_step*self.velocity[i]+0.5*self.time_step**2*self.acceleration[i]
def put_next_velocity(self):
for i in range(self.dimensions):
self.next_velocity[i] = self.velocity[i] + 0.5*self.time_step*(self.acceleration[i]+self.next_acceleration[i])
def update_position(self):
self.position = self.next_position.copy()
def update_velocity(self):
self.velocity = self.next_velocity.copy()
def update_acceleration(self):
self.acceleration = self.next_acceleration.copy()
def reset_acceleration(self):
self.acceleration = [0 for i in range(self.dimensions)]
def reset_future_acceleration(self):
self.next_acceleration = [0 for i in range(self.dimensions)]
def calculate_acceleration(self,other_particle):
"""
Increments the acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.position[i]-other[i])**2
distances.append(self.position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.acceleration[i] += acc[i]
def calculate_future_acceleration(self,other_particle):
"""
Increments the future acceleration of the particle due to the force from
a single other particle
"""
distances = []
other = other_particle.get_next_position()
distance_squared = 0
for i in range(self.dimensions):
distance_squared += (self.next_position[i]-other[i])**2
distances.append(self.next_position[i]-other[i])
distance = math.sqrt(distance_squared)
force = -self.G*self.mass*other_particle.get_mass()/distance_squared
acc = []
for i in range(self.dimensions):
acc.append(force*distances[i]/(distance*self.mass))
for i in range(self.dimensions):
self.next_acceleration[i] += acc[i]
def update_all(particleList):
for i in range(len(particleList)):
particleList[i].reset_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_position()
for i in range(len(particleList)):
particleList[i].reset_future_acceleration()
for j in range(len(particleList)):
if i != j:
particleList[i].calculate_future_acceleration(particleList[j])
for i in range(len(particleList)):
particleList[i].put_next_velocity()
for i in range(len(particleList)):
particleList[i].update_position()
particleList[i].update_velocity()
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
Alist = [[],[]]
Blist = [[],[]]
for i in range(10000):
Alist[0].append(partList[0].get_position()[0])
Alist[1].append(partList[0].get_position()[1])
Blist[0].append(partList[1].get_position()[0])
Blist[1].append(partList[1].get_position()[1])
update_all(partList)
plt.scatter(Alist[0],Alist[1],color="r")
plt.scatter(Blist[0],Blist[1],color="b")
plt.grid()
plt.show()
for item in partList:
print(item)
- , , .