Our task is to emulate the motion of N points in a closed perimeter, taking into account some laws of physics.

Laws of interaction:
Moving points interact with each other and with the walls.
Each point has a known position, velocity, mass, and radius.
There are no other influencing forces.

We need a model that allows us to work with time intervals of length dt
We will update the positions of all points after dt time and check for overlap
If we have an overlap (the event should have happened a little earlier than this time period), we rewind the time, calculate the collision, update the velocities and directions of the points, and continue the simulation.

We change the state only if something has happened
between collisions, the points move along their own trajectory
we focus only on the time points where collisions occur
create a priority queue containing collision events ordered by time
remove min and get the next collision
We have a method that can predict collisions in time:
it receives the position, velocity, and radius of a point and calculates when it will collide with another point or wall
If a collision occurs, we have to update the colliding objects.

Collision prediction:
we have a point of radius s at position (rx,ry)
the point is moving with velocity (vx,vy)
question: will it collide with the vertical wall and if so, when?

Initialization:
Fill the queue with all potential collisions with walls
Fill the queue with all potential collisions between points
potential because they may not happen, as the points may collide earlier and change direction
The main cycle:
we remove the event with the lowest priority
if the event is no longer valid, ignore it (for this we have count)
execute the event
calculate new collisions and add them to the queue
update the positions of all particles by time t


Code examples:
import functools
import random
import warnings

import matplotlib.pyplot as plt
import matplotlib.animation as animation

# disable warnings
warnings.filterwarnings("ignore")
# matplotlib configuration
# bounds of the room
x_lim = (0, 30)
y_lim = (0, 20)
# scaler for bigger points in plots
PLT_RADIUS_SCALER = 25

# set up plot to draw animations
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=x_lim, ylim=y_lim)
# disable x and y axes labels
ax.set_yticks([])
ax.set_xticks([])


class Ball:
    def __init__(self, rx, ry, vx, vy, radius):
        self._rx = rx
        self._ry = ry
        self._vx = vx
        self._vy = vy
        self._radius = radius

        # scatter plot for drawing point
        self.scatter, = ax.plot([], [], 'o', markersize=self._radius * PLT_RADIUS_SCALER)

    def move(self, dt):
        if (self._rx + self._vx * dt < x_lim[0] + self._radius) or \
                (self._rx + self._vx * dt > x_lim[1] - self._radius):
            self._vx = -self._vx

        if (self._ry + self._vy * dt < y_lim[0] + self._radius) or \
                (self._ry + self._vy * dt > y_lim[1] - self._radius):
            self._vy = -self._vy

        self._rx = self._rx + self._vx * dt
        self._ry = self._ry + self._vy * dt

        # updates point position in plot
        self.scatter.set_data([self._rx], [self._ry])


# this function moves all objects passing delta time, returns
# scatter plot of points that should be drawn
# the last parameter is time in seconds, it goes from matplotlib animator
def animate(objects, dt, _t):
    for obj in objects:
        obj.move(dt)

    return [obj.scatter for obj in objects]


if __name__ == "__main__":
    # 10 milliseconds delta t
    delta_t = 10
    n_balls = 5
    balls = []

    for i in range(n_balls):
        x = random.randint(1, 25)
        y = random.randint(1, 15)
        vx = random.randint(1, 20)
        vy = random.randint(1, 20)
        radius = random.random()
        print(f"Creating a new Ball({x}, {y}, {vx}, {vy}, {radius})")
        balls.append(Ball(x, y, vx, vy, radius))

    # interval in milliseconds
    ani = animation.FuncAnimation(
        fig, functools.partial(animate, balls, delta_t / 1_000),
        interval=delta_t, blit=True, save_count=100,
    )

    plt.show()

import functools
import math
import random
import warnings

import matplotlib.pyplot as plt
import matplotlib.animation as animation

# disable warnings
warnings.filterwarnings("ignore")
# matplotlib configuration
# bounds of the room
x_lim = (0, 30)
y_lim = (0, 20)
# scaler for bigger points in plots
PLT_RADIUS_SCALER = 25

# set up plot to draw animations
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=x_lim, ylim=y_lim)
# disable x and y axes labels
ax.set_yticks([])
ax.set_xticks([])


class Particle:
    def __init__(self, rx, ry, vx, vy, radius, mass):
        self._rx = rx
        self._ry = ry
        self._vx = vx
        self._vy = vy
        self._radius = radius
        self._mass = mass
        self._collisions_count = 0

        # scatter plot for drawing point
        self._scatter, = ax.plot([], [], 'o', markersize=self._radius * PLT_RADIUS_SCALER)

    @property
    def rx(self):
        return self._rx

    @property
    def ry(self):
        return self._ry

    @property
    def vx(self):
        return self._vx

    @property
    def vy(self):
        return self._vy

    @vx.setter
    def vx(self, value):
        self._vx = value

    @vy.setter
    def vy(self, value):
        self._vy = value

    @property
    def radius(self):
        return self._radius

    @property
    def mass(self):
        return self._mass

    @property
    def collisions_count(self):
        return self._collisions_count

    @collisions_count.setter
    def collisions_count(self, value):
        self._collisions_count = value

    @property
    def scatter(self):
        return self._scatter

    def move(self, dt):
        # Ball implementation below
        # if (self._rx + self._vx * dt < x_lim[0] + self._radius) or \
        #         (self._rx + self._vx * dt > x_lim[1] - self._radius):
        #     self._vx = -self._vx
        #
        # if (self._ry + self._vy * dt < y_lim[0] + self._radius) or \
        #         (self._ry + self._vy * dt > y_lim[1] - self._radius):
        #     self._vy = -self._vy
        #
        # self._rx = self._rx + self._vx * dt
        # self._ry = self._ry + self._vy * dt

        # updates point position in plot
        self._scatter.set_data([self._rx], [self._ry])

    def bounce_off(self, that):
        dx = that.rx - self.rx
        dy = that.ry - self.ry

        dvx = that.vx - self.vx
        dvy = that.vy - self.vy

        dvdr = dx * dvx + dy * dvy

        dist = self.radius + that.radius

        j = 2 * self.mass * that.mass * dvdr / ((self.mass + that.mass) * dist)

        j_x = j * dx / dist

        j_y = j * dy / dist

        self.vx += j_x / self.mass
        self.vy += j_y / self.mass

        that.vx -= j_x / that.mass
        that.vy -= j_y / that.mass

        self.collisions_count += 1
        that.collisions_count += 1

    def time_to_hit(self, that):
        if self == that:
            return math.inf

        dx = that.rx - self.rx
        dy = that.ry - self.ry

        dvx = that.vx - self.vx
        dvy = that.vy - self.vy

        dvdr = dx * dvx + dy * dvy

        if dvdr > 0:
            return math.inf

        dvdv = dvx * dvx + dvy * dvy
        drdr = dx * dx + dy * dy

        sigma = self.radius + that.radius

        d = (dvdr * dvdr) - dvdv * (drdr - sigma * sigma)

        if d < 0:
            return math.inf

        return -(dvdr + math.sqrt(d)) / dvdv


# this function moves all objects passing delta time, returns
# scatter plot of points that should be drawn
# the last parameter is time in seconds, it goes from matplotlib animator
def animate(objects, dt, _t):
    for obj in objects:
        obj.move(dt)

    return [obj.scatter for obj in objects]


if __name__ == "__main__":
    # 10 milliseconds delta t
    delta_t = 10
    n_balls = 5
    balls = []

    for i in range(n_balls):
        x = random.randint(1, 25)
        y = random.randint(1, 15)
        vx = random.randint(1, 20)
        vy = random.randint(1, 20)
        mass = random.random()
        radius = random.random()
        print(f"Creating a new Particle({x}, {y}, {vx}, {vy}, {radius}, {mass})")
        balls.append(Particle(x, y, vx, vy, radius, mass))

    # interval in milliseconds
    ani = animation.FuncAnimation(
        fig, functools.partial(animate, balls, delta_t / 1_000), interval=delta_t, blit=True, save_count=100,
    )

    plt.show()