Решение задачи асболютно упругого соударения двух шаров?
Я использую такое решение для эмуляции движения и соударения пятнадцати шаров (Object Pascal):

for i:=1 to 15 do   // просчитываем соударения шаров
for j:=i+1 to 16 do begin
 dist:=sqrt(sqr(balls[i].x-balls[j].x)+sqr(balls[i].y-balls[j].y)); //расстояние между центрами шаров
 if dist<diametr then begin //если расстояние меньше диаметра, значит есть факт соударения
  a:=balls[i].x-balls[j].x; //вспомогательные переменные типа extended
  b:=balls[i].y-balls[j].y;
  p1:=a*b/sqr(dist);
  p2:=sqr(a/dist);
  p3:=sqr(b/dist);
  d1:=balls[i].dy*p1+balls[i].dx*p2-balls[j].dy*p1-balls[j].dx*p2;
  d2:=balls[i].dx*p1+balls[i].dy*p3-balls[j].dx*p1-balls[j].dy*p3;
  balls[i].dx:=balls[i].dx-d1; //меняем значение приращения координаты шаров при движении
  balls[i].dy:=balls[i].dy-d2;
  balls[j].dx:=balls[j].dx+d1;
  balls[j].dy:=balls[j].dy+d2;

  p3:=(diametr-dist)/2; //при соударении шары всегда "проникают" друг в друга, поэтому раздвигаем их
  p1:=p3*(a/dist);
  p2:=p3*(b/dist);
  balls[i].x:=balls[i].x+p1;
  balls[i].y:=balls[i].y+p2;
  balls[j].x:=balls[j].x-p1;
  balls[j].y:=balls[j].y-p2;
 end;
end;

for i:=1 to 15 do
 balls[i].x:=balls[i].x+balls[i].dx;  //эмулируем движение
 balls[i].y:=balls[i].y+balls[i].dy;
end;

Пока такое решение, но его нужно переделывать…

Кодclc; clear all; close all;

% Start points
x1 = 0; y1 = 0; r1 = 2;
x2 = 6; y2 = 1; r2 = 2;

% Absolute Velocity
v1 = 6;
v2 = 3;

% Velocity vectors
L1 = rand(1,2);
L2 = rand(1,2);

% Time
t = 0:0.02:10;

% Vectors with coordinates by steps (time)
vec1(1,:) = x1 + v1*L1(1)*t; vec1(2,:) = y1 + v1*L1(2)*t;
vec2(1,:) = x2 + v2*L2(1)*t; vec2(2,:) = y2 + v2*L2(2)*t;

% Draw circles
ang = 0:0.01:2*pi;
xp1 = r1*cos(ang); yp1 = r1*sin(ang);
xp2 = r2*cos(ang); yp2 = r2*sin(ang);

for i=1:length(vec1)
    plot(vec1(1,i)+xp1,vec1(2,i)+yp1);
    hold on
    plot(vec2(1,i)+xp2,vec2(2,i)+yp2);
    hold off
    axis([min(min(vec1(1,:), vec2(1,:))), max(max(vec1(1,:), vec2(1,:))), min(min(vec1(2,:), vec2(2,:))), max(max(vec1(2,:), vec2(2,:)))])
    pause(0.02) % delay
    
    if sqrt((vec1(1,i)-vec2(1,i)).^2 + (vec1(2,i)-vec2(2,i)).^2) < r1+r2
        
        disp('COLLAPSE!')
        
        L1 = rot90(L1,1);
        L2 = rot90(L2,1);
        tt = 0:0.02:10;
        
        newvec1(1,:) = vec1(1,i) + v1*L1(1)*tt;
        newvec1(2,:) = vec1(2,i) + v1*L1(2)*tt;
        newvec2(1,:) = vec2(1,i) + v2*L2(1)*tt;
        newvec2(2,:) = vec2(2,i) + v2*L2(2)*tt;
        
        for j=i:length(vec1)
            plot(newvec1(1,j-i+1)+xp1,newvec1(2,j-i+1)+yp1);
            hold on;
            
            plot(newvec2(1,j-i+1)+xp2,newvec2(2,j-i+1)+yp2);
            hold off;
            
            axis([min(min(vec1(1,:), vec2(1,:))), max(max(vec1(1,:), vec2(1,:))), min(min(vec1(2,:), vec2(2,:))), max(max(vec1(2,:), vec2(2,:)))])
            pause(0.1) % delay
        end
        
        return
        
    end
end

