I have this code to emulate the movement of N points in a closed perimeter, taking into account some laws of physics. It works, but the animation is delayed, frame-by-frame and slow. How can I fix it?
import functools
import math
import random
import warnings
import heapq
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):
        self._rx += self._vx * dt
        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

        # Force factor
        magnitude = 2 * self.mass * that.mass * dvdr / ((self.mass + that.mass) * dist)

        # Normalization factor
        fx = magnitude * dx / dist
        fy = magnitude * dy / dist

        # Update velocities
        self._vx += fx / self.mass
        self._vy += fy / self.mass
        that._vx -= fx / that.mass
        that._vy -= fy / 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

    def time_to_hit_vertical_wall(self):
        if self._vx > 0:
            return (x_lim[1] - self._radius - self._rx) / self._vx
        elif self._vx < 0:
            return (x_lim[0] + self._radius - self._rx) / self._vx
        else:
            return math.inf

    def time_to_hit_horizontal_wall(self):
        if self._vy > 0:
            return (y_lim[1] - self._radius - self._ry) / self._vy
        elif self._vy < 0:
            return (y_lim[0] + self._radius - self._ry) / self._vy
        else:
            return math.inf

    def bounce_off_vertical_wall(self):
        self._vx = -self._vx
        self.collisions_count += 1

    def bounce_off_horizontal_wall(self):
        self._vy = -self._vy
        self.collisions_count += 1

class Event:
    def __init__(self, time, a, b):
        self.time = time  # Time when the event occurs
        self.a = a        # First particle (can be None for wall collisions)
        self.b = b        # Second particle (can be None for wall collisions)
        self.count_a = a.collisions_count if a else -1
        self.count_b = b.collisions_count if b else -1

    def __lt__(self, other):  # For the priority queue
        return self.time < other.time

    def is_valid(self):
        if self.a and self.a.collisions_count != self.count_a:
            return False
        if self.b and self.b.collisions_count != self.count_b:
            return False
        return True

# Initialize particles
delta_t = 100  # milliseconds
n_balls = 5
balls = []

for i in range(n_balls):
    x = random.uniform(x_lim[0] + 1, x_lim[1] - 1)
    y = random.uniform(y_lim[0] + 1, y_lim[1] - 1)
    vx = random.uniform(-2, 2)
    vy = random.uniform(-2, 2)
    mass = random.uniform(0.5, 1.5)
    radius = random.uniform(2.2, 2.5)
    print(f"Creating a new Particle({x}, {y}, {vx}, {vy}, {radius}, {mass})")
    balls.append(Particle(x, y, vx, vy, radius, mass))

# Event-driven simulation
event_queue = []
heapq.heapify(event_queue)
current_time = 0.0
simulation_limit = 100.0  # seconds

def predict_events(particle, current_time, limit):
    # Predict particle's collisions with other particles
    for other in balls:
        if particle != other:
            dt = particle.time_to_hit(other)
            event_time = current_time + dt
            if 0 <= dt <= limit:
                event = Event(event_time, particle, other)
                heapq.heappush(event_queue, event)
    # Predict particle's collisions with vertical walls
    dt = particle.time_to_hit_vertical_wall()
    event_time = current_time + dt
    if 0 <= dt <= limit:
        event = Event(event_time, particle, None)
        heapq.heappush(event_queue, event)
    # Predict particle's collisions with horizontal walls
    dt = particle.time_to_hit_horizontal_wall()
    event_time = current_time + dt
    if 0 <= dt <= limit:
        event = Event(event_time, None, particle)
        heapq.heappush(event_queue, event)

# Schedule initial events
for particle in balls:
    predict_events(particle, current_time, simulation_limit)

def simulate():
    global current_time
    if not event_queue:
        return
    event = heapq.heappop(event_queue)
    if not event.is_valid():
        simulate()
        return  # Skip invalid events

    # Move all particles to the time of the event
    dt = event.time - current_time
    for particle in balls:
        particle.move(dt)
    current_time = event.time

    # Process the event
    a, b = event.a, event.b
    if a and b:
        # Particle-particle collision
        a.bounce_off(b)
    elif a and not b:
        # Particle-vertical wall collision
        a.bounce_off_vertical_wall()
    elif b and not a:
        # Particle-horizontal wall collision
        b.bounce_off_horizontal_wall()

    # After processing, predict new events involving the affected particles
    if a:
        predict_events(a, current_time, simulation_limit)
    if b:
        predict_events(b, current_time, simulation_limit)

def animate(_):
    simulate()
    # Return updated scatter plots for redrawing
    return [obj.scatter for obj in balls]

ani = animation.FuncAnimation(
    fig, animate, interval=delta_t, blit=True)

plt.show()