Python N-body simulation code giving incorrect answer
<p>I have written some python code to solve the N-body problem using the Euler method. The code runs without problems and seems to give a reasonable answer (e.g. if there are two particles then the start moving towards each other). However when I run this simulation over a large number of iterations I see that the particles (say I run it with two particles) pass by each other (I do not consider collisions) and keep going in their directions indefinitely. This violates the conservation of energy and so there must be a flaw in my code but I am unable to find it. Could anyone please find it and explain my mistake.</p>

<p>Thank you.</p>

<p>Thanks to @samgak for pointing out that I was updating the particles twice. I have now fixed this but the problem still keeps coming. I have also replicated  the output I get when I run this simulation with two stationary particles at (0,0) and (1,0) with a time step of 1 second and 100000 iterations:</p>

<p>Particle with mass: 1 and position: [234.8268420043934, 0.0] and velocity: [0.011249111128594091, 0.0]</p>

<p>Particle with mass: 1 and position: [-233.82684200439311, 0.0] and velocity: [-0.011249111128594091, 0.0]</p>

<p>Also thanks to @PM2Ring for pointing out some optimizations I could make and the perils of using the Euler method.</p>

<p>Code:</p>

<pre><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 #fixed throughout the simulation
        self.time_interval = 10**0 #fixed throughout the simulation, gives the interval between updates
        self.mass = mass
        self.position = position #should be a list
        self.velocity = velocity #should be a list
        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()      
#the list of particles
partList = [Particle(1,[0,0],[0,0]),Particle(1,[1,0],[0,0])]
#how many iterations I perform
for i in range(100000):
    update_particles(partList)
#prints out the final position of all the particles
for item in partList:
    print(item)
</code></pre>

<p>------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------Further Edit:</p>

<p>I decided to implement the Leapfrog method and I have developed some code that once again runs and seems to work well (at least in the command line). However when I added plotting functionality and analysed it there seemed to be another problem. Again the system seemed to go too far and the energy again increased without bound. I have attached a picture of the output I get to showcase the problem. If I again had just two particles of equal mass they again pass each other and continue away from each other without stopping. Thus there must be a bug in my code I am not finding.</p>

<p>If anyone can help it will be much appreciated.</p>

<p>My Code:</p>

<pre><code>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)
</code></pre>

<p><a href="https://i.stack.imgur.com/0PHtE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/0PHtE.png" alt="A zoomed in plot"></a>  </p>

<p>Could someone please tell me where is the error I am making in my code.</p>