Почему шарики такие неустойчивые?

Это физическая симуляция, ограничивающая шары в круглой области.

Я сделал исходный код в Scratch и преобразовал его в Python в Pygame. Когда я запускал симуляцию, все шары тряслись по сравнению с исходным кодом.

Я ограничил максимальную скорость до 20, но это не помогло.

Я создал подшаги для каждого кадра, но это тоже не помогло.

import pygame
import math
import random

screen = pygame.display.set_mode((16*120,9*120))
pygame.display.set_caption('')
clock = pygame.time.Clock()

mx = 16*60
my = 9*60

global x
global y
x = []
y = []
xv = []
yv = []
prevx = []
prevy = []

for i in range(20):
    x.append(random.randint(-200,200))
    y.append(random.randint(-200,200))
    xv.append(0)
    yv.append(0)
    prevx.append(0)
    prevy.append(0)

r = 25
size = 300

sub = 20
maxvel = 20

global dist
global dx
global dy
global d

#points at something
def pointat(px,py):
    global d
    dx = px-x[i]
    dy = py-y[i]

    if dy == 0:
        if dx < 0:
            d = -90
        else:
            d = 90
    else:
        if dy < 0:
            d = 180+math.atan(dx/dy)
        else:
            d = math.atan(dx/dy)

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            exit()
    screen.fill((10,10,10))
    pygame.draw.circle(screen,(100,100,100),(mx,my),size)
    for i in range(len(x)):
        pygame.draw.circle(screen,(255,255,255),(mx+x[i],my+y[i]*-1),r)

    for i in range(len(x)):
        prevx[i] = x[i]
        prevy[i] = y[i]
        for j in range(sub):

            x[i] += xv[i]/sub
            y[i] += yv[i]/sub

            y[i] -= 1/sub
            
            a = 0
            for k in range(len(x)-1):
                if a == i:
                    a += 1
                dist = math.sqrt(((x[i]-x[a])*(x[i]-x[a]))+((y[i]-y[a])*(y[i]-y[a])))
                if dist < r*2:
                    pointat(x[a],y[a])
                    x[i] += (math.sin(d)*(dist-(r*2)))
                    y[i] += (math.cos(d)*(dist-(r*2)))
                    x[a] -= (math.sin(d)*(dist-(r*2)))
                    y[a] -= (math.cos(d)*(dist-(r*2)))
                
                dist = math.sqrt((x[i]*x[i])+(y[i]*y[i]))
                if dist > size-r:
                    pointat(0,0)
                    x[i] += (math.sin(d)*(dist-(size-r)))
                    y[i] += (math.cos(d)*(dist-(size-r)))
                a += 1

        xv[i] = x[i]-prevx[i]
        yv[i] = y[i]-prevy[i]
        if xv[i] > maxvel:
            xv[i] = maxvel
        if xv[i] < -maxvel:
            xv[i] = -maxvel
        if yv[i] > maxvel:
            yv[i] = maxvel
        if yv[i] < -maxvel:
            yv[i] = -maxvel

    pygame.display.update()
    clock.tick(60)