Recently I decided to write a discrete collision simulation from scratch to help learn python and for personal enjoyment. I wrote this program to create a video by saving each graph after a discrete time update as an image in a folder, then compiling all those images together (using ffmpeg) resulting in a video. At first I wanted to see what would happen if I simulated a ball in a box under the influence of gravity. I know what the answer should have come out to look like, I've played with bouncy balls as a kid. However whenever the ball gets near the end of its bouncing, it will start to clip through the floor at a stuck velocity rather than coming to a rest on the floor and I can not for the life of me figure out why.
I also tried adding in another particle to the simulation and took out gravity so that these particles bounce around in free space, however the same problem will occur sometimes when the particles collide. They stick together rather than bouncing off, much like how a single ball will stick into the floor.
This is the code that I used. I tried to comment it as best as I could. I'm not a very proficient coder so I apologize if there are conventional errors.
I think the problem is in the "handleBoxCollision()" function. My best guess is that at small numbers the velocity after collision isn't enough to get the particle out on the next frame, so it stays stuck in the wall getting velocity flipped again and again and again and it can never get out. I have tried changing the velocity calculation to an absolute value (since when the ball hits the ground it has negative velocity, then switches to positive) but the ball will still continue to go through the floor when it is supposed to be at rest. The code I used is included below, as well as a discord link that goes to the video. where Does anyone have any insight?
video: https://cdn.discordapp.com/attachments/368098721994506240/1104520524488527883/video.webm
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import math
import os
#Class defining a particle used in the simulation
class Particle:
def __init__(self, mass, position, velocity, acceleration):
self.mass = mass
self.radius = math.sqrt(self.mass/(math.pi*1.5))
self.position = position
self.velocity = velocity
self.acceleration = acceleration
self.KE = (1/2)*self.mass*np.dot(self.velocity,self.velocity)
self.left = position[0]-self.radius
self.right = position[0]+self.radius
self.top = position[1]+self.radius
self.bottom = position[1]-self.radius
#just defining the box that the simulation takes place in
class Box:
def __init__(self):
self.left = -10
self.right = 10
self.bottom = -10
self.top = 10
#Function that detects if there is a colision between the particle and the box
#This is where i think theres a problem but cant for the life of me figure out what it is
def handleBoxCollision():
#cor is the coefficient of restitution
cor = 0.8
if p.left <= box.left or p.right >= box.right:
p.velocity[0]=-cor*p.velocity[0]
if p.bottom <= box.bottom or p.top >= box.top:
p.velocity[1]=-cor*p.velocity[1]
#Since this is a discreet colission simulation this function is for updating the state of the simulation
def update(dt):
p.velocity = p.velocity+(p.acceleration*dt)
p.position = p.position+(p.velocity*dt)
p.left = p.position[0]-p.radius
p.right = p.position[0]+p.radius
p.top = p.position[1]+p.radius
p.bottom = p.position[1]-p.radius
handleBoxCollision()
#Because I run this simulation many times I first delete all the contents in my directory
#Then as simulation runs it saves each updated graph as a frame. Compile all frames into a video
#I took out the actual directory on my computer for privacy
dir = 'C:/PathToImageFile'
for f in os.listdir(dir):
os.remove(os.path.join(dir, f))
#Initial mass, position, velocity and acceleration
mass = 10
position = np.array([0,0])
velocity = np.array([5,0])
acceleration = np.array([0,-9.8])
#time step = 1/framerate
dt = 1/60
p = Particle(mass, position, velocity, acceleration)
box = Box()
#Run this loop for however many frames I want. In this case 600
for i in range(600):
figure, axes = plt.subplots()
update(dt)
#left of box
plt.plot([-10,-10],[-10,10],color='black')
#right of box
plt.plot([10,10],[-10,10],color='black')
#top of box
plt.plot([-10,10],[10,10],color='black')
#bottom of box
plt.plot([-10,10],[-10,-10],color='black')
cc = plt.Circle((p.position[0] ,p.position[1]), p.radius)
plt.scatter(p.position[0],p.position[1])
plt.xlim(-11,11)
plt.ylim(-11,11)
plt.text(-10,-10.7,"time="+str(i*dt))
plt.text(-10,10.3,"velocity="+str(p.velocity[1]))
axes=plt.gca()
axes.add_artist(cc)
axes.set_aspect(1)
figure.savefig('/PathToImageFile'+str(i)+'.png')
plt.close('all')
The problem is that you're allowing particle.position[1] to move below box.bottom + particle.radius. Then after handleBoxCollision() reverses AND reduces (due to cor < 1) the velocity in that direction, the next tick doesn't move the particle back inside of the box bounds, so it just keeps reversing directions while stuck along the edge. If you set cor = 1 so that velocity is never reduced, it will keep bouncing as expected.
I made two changes below:
handleBoxCollision() I added the lines: if p.bottom <= box.bottom:
p.position[1] = box.bottom + p.radius
This will always force the particle back inside of the box. I only
handled the bottom edge as it is the most relevant for your gravity simulation. Adding similar resets for the other edges is up to you.
I didn't make a video, but based on checking the last 50-60 images, it seems to have resolved the issue.
An alternative solution would be to handle this within update and prevent the particle position from getting too close to the edges in the first place, by setting box.side +/- particle.radius as a hard limit of particle.position.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import math
import os
#Class defining a particle used in the simulation
class Particle:
def __init__(self, mass, position, velocity, acceleration):
self.mass = mass
self.radius = math.sqrt(self.mass/(math.pi*1.5))
self.position = position
self.velocity = velocity
self.acceleration = acceleration
self.KE = (1/2)*self.mass*np.dot(self.velocity,self.velocity)
# Use properties for the particle boundaries
@property
def left(self):
return self.position[0] - self.radius
@property
def right(self):
return self.position[0] + self.radius
@property
def top(self):
return self.position[1] + self.radius
@property
def bottom(self):
return self.position[1] - self.radius
#just defining the box that the simulation takes place in
class Box:
def __init__(self):
self.left = -10
self.right = 10
self.bottom = -10
self.top = 10
#Function that detects if there is a colision between the particle and the box
#This is where i think theres a problem but cant for the life of me figure out what it is
def handleBoxCollision():
#cor is the coefficient of restitution
cor = 0.8
if p.left <= box.left or p.right >= box.right:
p.velocity[0]=-cor*p.velocity[0]
if p.bottom <= box.bottom or p.top >= box.top:
p.velocity[1]=-cor*p.velocity[1]
if p.bottom <= box.bottom:
p.position[1] = box.bottom + p.radius
#Since this is a discreet colission simulation this function is for updating the state of the simulation
def update(dt):
p.velocity = p.velocity+(p.acceleration*dt)
p.position = p.position+(p.velocity*dt)
handleBoxCollision()
#Because I run this simulation many times I first delete all the contents in my directory
#Then as simulation runs it saves each updated graph as a frame. Compile all frames into a video
#I took out the actual directory on my computer for privacy
dir = './particle_image'
for f in os.listdir(dir):
os.remove(os.path.join(dir, f))
#Initial mass, position, velocity and acceleration
mass = 10
position = np.array([0,0])
velocity = np.array([5,0])
acceleration = np.array([0,-9.8])
#time step = 1/framerate
dt = 1/60
p = Particle(mass, position, velocity, acceleration)
box = Box()
#Run this loop for however many frames I want. In this case 600
for i in range(600):
figure, axes = plt.subplots()
update(dt)
#left of box
plt.plot([-10,-10],[-10,10],color='black')
#right of box
plt.plot([10,10],[-10,10],color='black')
#top of box
plt.plot([-10,10],[10,10],color='black')
#bottom of box
plt.plot([-10,10],[-10,-10],color='black')
cc = plt.Circle((p.position[0] ,p.position[1]), p.radius)
plt.scatter(p.position[0],p.position[1])
plt.xlim(-11,11)
plt.ylim(-11,11)
plt.text(-10,-10.7,"time="+str(i*dt))
plt.text(-10,10.3,"velocity="+str(p.velocity[1]))
axes=plt.gca()
axes.add_artist(cc)
axes.set_aspect(1)
figure.savefig('./particle_image/'+str(i)+'.png')
plt.close('all')
Suggestions for other possible improvements:
Update handleBoxCollision() to take the particle and box as arguments. This will make it easier to support multiple particles and won't rely on global variables. It might work best as a method of the Box class, then you just pass it the particle to check collisions for. Or make it a method of the particle class and just pass the box in; but that feels backwards for some reason.
Likewise, update() should probably be a method of the Particle class.