Метод Верле́ — один из самых элегантных и простых численных способов решать уравнения движения. Его можно встретить и в молекулярной физики, и в геймдеве (cloth sim). Недавно я сделал короткое видео с его наглядной демонстрацией (см. YouTube Shorts под катом и зеркало на GitHub).

В этом посте я хочу показать, как идея «чисто позиционного» интегрирования без явного использования скоростей превращается в рабочую анимацию, начиная с простых примеров, заканчивая сетками и игрой в бильярд.

Я уже публиковал пост о методе Верле ✂️, кратко упомянув идеи, лежащие в его основе.

На Хабре есть пара статей на подобную тему, но довольно мало за последние 10 лет

• Симуляция ткани методом Стёрмера — Верле

• Верле: разрешаем коллизии

• Рэгдолл физика своими руками

Основы

Все заслуги принадлежат этим людям:

Луп Верле

французский физик, пионер в Молекулярно-кинетической теории (1967)

Карл Стёрмер

норвежский математик и физик (1907)

Как показано в видео, простая система из дискретизированных уравнений движения материальной точки

приводит нас к следующему численному решению как кинематических, так и динамических задач:

где — номер итерации или шага по времени (). Как видно, здесь нет скорости, а только координаты и силы.

Превращаем уравнения в анимацию

Для демонстрации я использую WLJS Notebook ? в силу удобства для меня реализации динамики и анимаций в блокноте

Это уравнение удобно переводится в код:

With[{uid = CreateUUID[], ball = Unique["ball"]},

  (* начальное положение шарика x3 *)  

  ball = Table[{0, 0}, {3}];

  (* функция обновления координат *) 
     
  EventHandler[uid, Function[Null,
    (* метод Верле *)  
                             
    ball[[3]] = ball[[2]];
    ball[[2]] = ball[[1]];

    ball[[1]] = 2 ball[[2]] - ball[[3]] + {0, -1} 0.01;
    If[ball[[1, 2]] < -1.0, ball[[{1, 2}]] = ball[[{2, 1}]]];

    (* дергаем еще раз, чтобы обновить шарик на экране *)      
                             
    ball = ball;
  ]];

  (* динамическая графика *)
     
  Graphics[{
    {Gray, Line[{{-1, 0}, {1, 0}}]},
    Point[ball[[1]] // Offload],
    AnimationFrameListener[ball // Offload, "Event" -> uid]
  }, PlotRange -> {{-1, 1}, {-1, 1}}]
]

Все три , , одного подпрыгивающего шарика хранятся в символе ball. Для решения этой задачи уравнение Верле преобразуется раскрывается в

где последний член соответствует действию силы тяжести на шар единичной массы помноженной на квадрат постоянной времени, ака .

Граничные условия: если шарик пересекает поверхность земли, мы меняем местами и — это неявно означает «переворот» вектора скорости

Демонстрация

Потери на тепло

Чтобы учесть неупругое столкновение, можно уменьшить отражённый импульс (скорость), скажем в 0.95 - 0.8 раз

  If[ball[[1, 2]] < -1.0, 
    With[{vel = ball[[1]]-ball[[2]]},
      ball[[2]] = ball[[1]];
      ball[[1]] = ball[[1]] - 0.95 vel;
    ];
  ];

Точность численных методов

Даже при малом шаге времени возможны ошибки выхода за границы. Чаще всего после сильного удара:

Да, здесь приходится «жульничать», как это делают почти все разработчики физических движков.

В таких случаях можно просто прижать (clamp) одну из координат к границам области вручную мимо любой физики

В коде это станет такой поправкой

  If[ball[[1, 2]] < -1.0, 
    With[{vel = ball[[1]]-ball[[2]]},
      ball[[2]] = ball[[1]];
      ball[[1]] = Clip[ball[[1]] - 0.95 vel, {-1, Infinity}]; (* <-- *)
    ];
  ];

Если у вас есть более элегантное решение, напишите об этом в комментариях ??‍♂️

Общие граничные условия

Если поверхность наклонена, нужно отражать скорость относительно нормали. Это хорошее упражнение, которое можно попробовать сделать самостоятельно

Здесь же я просто приведу ответ:

Шарик в ободе

Скажем

Circle[] // Graphics

будет нашим граничным условием. Т.е. делаем все то же самое, только отражаем вектор скорости от поверхности заданной внутренней "стороной" обода

With[{
  p = Unique["bbb"], 
  frame = CreateUUID[],
  circ = {0,0}
},
  p = Table[{0,0}, {3}];

  EventHandler[frame, Function[Null,
    Module[{b = p, \[Delta]t = 0.01},
      Do[
        If[Norm[b[[1]]-circ] > 1.0,
          With[{n = -Normalize[b[[1]]-circ], delta = b[[1]]-b[[2]]},
            b[[2]] = b[[1]];
            b[[1]] = b[[1]] -  n (n.delta);
          ]
        ];

        b[[3]] = b[[2]];
        b[[2]] = b[[1]];
        b[[1]] = 2 b[[2]] - b[[3]] + {0,-1} \[Delta]t \[Delta]t;
      , {5}];
    
      p = b;
    ];
  ]];

  Graphics[{
    Circle[circ, 1.0], 
    Lighter[Lighter[Red]], Disk[p[[3]]//Offload, 0.01],
    Lighter[Red], Disk[p[[2]]//Offload, 0.01],
    Red, Disk[p[[1]]//Offload, 0.01],
    AnimationFrameListener[p//Offload, "Event"->frame]
  }, PlotRange->{{-1,1}, {-1,1}}, AspectRatio->1]
]

Правда если обод статичен, мы никак не увидим разницы между ним и плоскостью, когда шарик падает строго вертикально... А давайте сделаем обод подвижным

With[{
  p = Unique["bbb"], 
  frame = CreateUUID[],
  circ = Unique["bbb"]
},
  circ = {0,0};
  p = Table[{0,0}, {3}];

  EventHandler[frame, Function[Null,
    Module[{b = p, \[Delta]t = 0.01},
      Do[
        If[Norm[b[[1]]-circ] > 1.0,
          With[{n = -Normalize[b[[1]]-circ], delta = b[[1]]-b[[2]]},
            b[[2]] = b[[1]];
            b[[1]] = b[[1]] -  n (n.delta);
          ]
        ];

        b[[3]] = b[[2]];
        b[[2]] = b[[1]];
        b[[1]] = 2 b[[2]] - b[[3]] + {0,-1} \[Delta]t \[Delta]t;
      , {5}];
    
      p = b;
    ];
  ]];

 {
  Graphics[{
    Circle[circ//Offload, 1.0], 
    Lighter[Lighter[Red]], Disk[p[[3]]//Offload, 0.01],
    Lighter[Red], Disk[p[[2]]//Offload, 0.01],
    Red, Disk[p[[1]]//Offload, 0.01],
    AnimationFrameListener[p//Offload, "Event"->frame]
  }, PlotRange->{{-1,1}, {-1,1}}, AspectRatio->1],

  (* новый элемент ввода *)            
  EventHandler[InputJoystick[], Function[xy, circ = xy]]
 } // Column
]

Всё работает! Почти... Но чего-то не хватает. Сам обод (или окружность, если угодно) тоже имеет скорость и импульс. Однако сейчас мы не передаём их шарику.

Передача импульса

Если всё сложно — переходите в другую систему отсчёта (СО).

Перенесём вычисления в систему отсчёта обода (движущуюся вместе с ним)

В той системе скорость обода равна нулю, он статичен. Мы решаем задачу с граничными условиями как и раньше, а затем возвращаемся в исходную (стационарную) систему отсчёта

u[[1]] = u[[1]] - rimVel;
                
With[{
  n = -Normalize[u[[1]]-circ[[2]]], 
  delta = u[[1]]-u[[2]]
},
  u[[2]] = u[[1]];
  u[[1]] = u[[1]] -  0.8 n (n.delta);
];
  
u[[1]] = u[[1]] + rimVel;

Полное решение для нескольких шариков

With[{
  p = Unique["bbb"],
  circ = Unique["bbb"],
  target = Unique["bbb"],
  frame = CreateUUID[]
},
  p = Table[RandomReal[{-1,1}, 2]0.001, {3}, {100}];
  circ = {{0,0}, {0,0}};
  target = {0,0};
  
  EventHandler[frame, Function[Null,
    Module[{b = p, \[Delta]t = 0.005},
      Do[
        circ[[1]] = circ[[1]] + \[Delta]t (target - circ[[1]]);

        (* граничные условия для каждого шарика *)
        
        b[[{1,2}]] = MapThread[If[Norm[#1-circ[[2]]] > 1.0, Module[{
          u = {#1, #2},
          vel = (circ[[1]] - circ[[2]])
        },

          (* меняем СО *)
                                                                   
          u[[1]] = u[[1]] - vel;
          
          
          With[{n = -Normalize[u[[1]]-circ[[2]]], delta = u[[1]]-u[[2]]},
            u[[2]] = u[[1]];
            u[[1]] = u[[1]] -  0.8 n (n.delta);
          ];

          (* меняем СО обратно *)    
                                                                   
          u[[1]] = u[[1]] + vel;

          (* клипуем координату если вышла за пределы *)  
                                                                   
          If[Norm[u[[1]]-circ[[2]]] > 1.0, 
            u[[1]] += (1.0 - Norm[u[[1]]-circ[[2]]]) Normalize[u[[1]]-circ[[2]]];
          ];
  
          u
        ], {#1, #2}]&, {b[[1]], b[[2]]}] // Transpose;

        (* обычный метод Верле *)
        
        b[[3]] = b[[2]];
        b[[2]] = b[[1]];
        b[[1]] = 2 b[[2]] - b[[3]] + Table[{0,-1}, {Length[b[[1]]]}] \[Delta]t \[Delta]t;
  
  
        circ[[2]] = circ[[1]];
      , {10}];
  
      circ = circ;
  
      p = b;
    ];
  ]];

  {
    Graphics[{
      Circle[circ[[1]] // Offload, 1.0], 
      Pink, Point[p[[1]] // Offload],
      AnimationFrameListener[p//Offload, "Event"->frame]
    }, PlotRange->{1.4{-1,1}, 1.4{-1,1}}, AspectRatio->1, ImageSize->Medium, TransitionType->None],
  
    EventHandler[InputJoystick[], Function[xy, target = xy;]]
  } // Row
]

Так то лучше!

Шарики и стержни (Constraint Algorithm)

При моделировании системы с шариками, соединенными стержнями, явно рассчитывая натяжения, передачу импульса - это быстро превращается в месиво. Здесь же метод Верле заметно упрощает все расчёты, так как вся возня со скоростями и силами происходит автоматически и за кадром.

Как я уже упоминал в предыдущем посту ✂️

... если хотим сохранить фиксированные расстояния между точками согласно графу для этого потребовалась бы нетривиальная работа со скоростью и координатой, но метод Верле оставляет нам только 

Жесткость таких связей задаётся не через дифференциальные уравнения, а специальным алгоритмом, который аппроксимирует реальные потенциалы этих связей. На этот счёт есть целая наука, как не решать диффуры, а обойтись манипуляциями координат. Здесь же представлен наиболее наивный и простой вариант.

Сначала мы выполняем нашу обычную итерацию по методу Верле, а затем применяем алгоритм - дисциплина со шпицрутенами. Поскольку каждая связь влияет сразу на обе вершины, «информация» (импульс, напряжение и т.д.) распространяется по всей системе узлов. Это приводит к тому, что если, например, наша сетка из связей ударяется о землю одной из вершин, всё тело останавливается и отскакивает обратно.

В дополнение приведу пример, похожий на тот, что вы (возможно) видели в YT Shorts:

Image to Mesh

Преобразование чёрно-белого изображения в сетку взаимосвязанных вершин

Возьмём такой пример изображения

dude = "1:eJztXWeIHVUU3uxmk6hBYwM7idg7iknEmkQJCRp1N7GAhU3ydo0ku7q7EQv+UYP/NFhTFLGggg2JgjVGNBqjoGAEEeIaBXuNdU053uOcj3fezZ15896Wd2fevfDtm51y55zz3XPb3DkzYW5XS/vIhoaGnmbzZ8aito5CeyP/O9b8mbV4UaF7wbwp3d1tN8xtMjvmC/gCIsoaRhg0GYwUNDrQ5IGceUej2D/t+U3CXTX3GSFo9EBvH6HtwpxMN2g3WGqw1mCd4F2D1Qaz1fnVcBIQD9Q9BxksMfiE0qVeg90q4AR+caTBWwZvynbwkx25mGbwk7L1NoMtgq0ObJPzNhkcm9KmqN8eVfd51JKjngH7TTH4U+zzr7J1ufSP/N4t+ZRre3A/ru/A65vWsawBfZumAeqA9nRfgx/ErlvTkKAS+w5zt4LK84H77WnwpcF2yWO9Ol5r29YSsF2v2OXfVAyUJvhRH0XtCGzuuh/qo+PU/ZiTpZY8WQF84XCDRwyepPT1dhIfK8UuW9IQ4EjwqTll7Iq24145H/xPUcdrbeNq+FirbLHeOlYNH8skr2r54OuYz6co3q7g4lTrXl8b7ELJfuUrIO/bog+3pf0GhyXYIQ0fyy0bVZrgHx9Ycrru9SxF3PXLdVdUKbsPgE5zRP9++V1epU6Dzcc6cvOB+1xkyf1KlXL7Avj0KIONFJUz9BnHU3FeqVI+BlpfJfEBeSYYfE9R+w8+zqLsz4PBhteJTn/L752yv7mKvFZIHgPl430q5QPl4wCDz+UcjFeq9WnfgPmG/Q1+FVug33g5ldo5LR93i40G2p4/TqU2Rv4vy3nwi88MdlW61Nqmg+UjnaLfP2IPHtMdbNkkCTjnLMkn7ZjcTuCxRckHP10kx/pFxj4lY1bH40mcPCz6oh74WB1Lwwls8gYV26NKErjgOd9myQ9c8HxYv5wDH55kyZ8XwNd57mGj2ATjK64fDkqpt+0jlfABf+J5r6Mln9HyO9XgLyr2bTndKMcqaeOyBJTtQwx+t3Tn9nMCFTlJqhtw7Hm5Ng0n4IL7E2fL9aOolAtO8NtVVMp/XoHyf6WyFWywSdkKtgA3uh3FfMtYitpabW9XwjzgLwZnWPKcT8V5YsjxGkVj8Ly03+WAMnca7divZNvx84UTHNeNUEAdwn00PO9wpe2C3yhqH8YY7GMwk6LxHRL89HWDnSTvPLXfaTnhOupTZROUc7Yv83KZwRGUXG+8qK6JS5w3zz1xn/sXtR/PrzixX9QjFwDqrnEUzf1q2+n6h7eZs/cc4PnJ7+S87eROfL3NFfpRSM8Z7Ez1ywWgdWdf+FDZSM/jlUtxXJRLXGedo2Soh/aiHPTzEK6XmJd3qDi3Um2Cj/EYh+dY+LnJMgH3ZU9RMtRL210J7DaC25aZYjtuI7B+x1Vv/Si2136C+ujphHtiHV2tdfcVsE+lZXWJxYHmZjNF401ur0ZRcQ1j4KEyYG5bj0FsjJTfMymqn+xxCNqeVskzb/MevgF+xOX+C7G95gRzuZj7CHwMPdAXWCMc6P6YfvYU2uzhAcr8rWJ7VxvyB0XrthoCJ0MO+Md4iuaiMEdic4K589CODy1Q3vkZHtbzuviYGPgYVj72MPg58FFzgI/dAx9eIPDhFwIffiHw4RcCH34h8OEXAh9+IfDhFwIffiHw4RcCH34h8OEXAh9+IfDhFwIffiHw4RcCH34h8OEXAh9+IfDhFwIffiEtH5OodM1vWBc3+DxwWW8WG+9Vho+THHkgXmy1cWDrHbBfXIykJD74fU72oT3ld3RMHvX8vlpauN6b4fVvJ1P0ju1tFL0DiPXsrsTvgfys0EfR+7R3GFxIxXc5wXvgxQ1tF8Tb5XgXeCdqsNJGyfsodb/Q9heBtpe3R4utfrVsiBgmSe+g63M18N6tfS3H8XjG4BgqclLvbYv2CY4N8I6yF2K2VvtObVxCPGUkjukxL0amegJ8gutz9gnEFklqGwY7aV4epCh+BstUT/WXbkOPp9J30ZPiYeh6R/tOOWxV58fljft+RMV4N/XAidbxKoreZyKK94lqYl8lJR3f3U6IC6VjEOW57gIXHIfpBWWHOPvY7wf2URSf7D6DqymKjTVJMNEB3j+Zon4yx3H8KiZvnVB/cR+MYw/m9fsGaCs4BhXivCT5BOzC44ceimLNjBmgDDw25Li66DPE9RVw7zWUz3Gj5kK32a6kyy1/h+MIKy+M2zH3kQZ2TAAu80vUfVz+ibrrEkuHrMPFhUt/2yeupeJYwBWrrBroOUb+fxYVY5S54guwTBsoPz6CMslx6ZO40Ptepah9gf2Gwg46zhnPdbne1eUEX82Dj4CL8ZT8DQ9wwXxdT6U+MdQyghOOw+TqD9s+ktWxO8o0j/PWJ3CBfdyXQdy+4e7PgJMnRBa7XYOMl9LwlZPBBOzJ49zXLJ1ceuq+fi3ireJbPhwX1tXfw76XqLSsZQWop1YpfdJwUctyB5kfc8iMOuwLGnh/e7gBm3aJDojlqRP04zrKl3kJzOtyHDnMy9iJ+7/4nkkWfMSOp+8aayFuFcdkRR/Kh/pYP5d3xeGAP19Apc8GfAbK+EMiu11P6XJ3rkdcaE4Y6ywOtC53WLr6Cvgvxx9GfPC4vvzVHnKh5XmAdvyGGMrRN1TZNz9rrUu3ZXubi5s85ULLNLuMDi3W+T5Cf1/Nnh9H2dpAxTiJPpYt1Ffch0Jsc1c8wLs85wO23ZuKz7xdsVg7PddDy7bUkp1TuW+D+QK0bdNFXldfkfVCf8rnviL4aLU44IQyxmtd9vCYE+hws7I9ErjhmOzNnsqvgbJyKLm/0Yp+yolUWhZ9gu6X2Hxg+2GP5XeB25A+kV37exbme5O+FYjt+z2W3wZ8xDX3Bn1u8VifpG9pYnu5da7PgIz3JOjzgpzjo7/nlY/5Dn30N7yxBti3NjFvfKDMTxbZXevn+XnyuMDHsPJxcgIfvGbsADnPt/573vhAeee5uM0OTlBnzaBS/nyB5sOeh8Mcw8oM8QFOuNx/JHq4+lgdnupkf9dZP4PCmOpaT2WPA+qg1TF8cBlr81QnPIfi589Yy6S/Lcc+fyBla90lbMzf4eP6SY/Vwc3pco5v9ZUuT7yWaZOS/SfZp8/JAmDjixUHXL70N195vYbP7+7C3vztgfME+2WQC1uf26k0fUtRXzgLernk813mJKDs8zr5aygaI+6fMb3wfizWNdVanoEib2UsD8Aa+bysqw4IyCRaxzU0NMzp7FnQ0VmYP6Ozt9BR6J7cOtLsnHpDb6G90Wz08H8tixcWenY2G9O6FnZ1t17TNq/Q2sT7p0+1ThprNjij7oWFtusWdHb8f2R29+LCfwMzo0I=" // Uncompress;

Оно случайно вышло прозрачным, но это легко починить

dude // Colorize

Запишем вспомогательную функцию, которая трансформирует изображение в сетку, триангулирует ее, а затем извлекает вершины и грани

SetAttributes[UndirectedEdge, Orderless];

getMesh[img_Image] := With[{
  mesh = TriangulateMesh[ImageMesh[img], MaxCellMeasure->{"Area"->45}]
}, Module[{
  vertices = MeshCoordinates[mesh],
  edges = (Cases[mesh//MeshCells, _Polygon, 2] /. {
  Polygon[{a_, b_, c_}] :> Sequence[UndirectedEdge[a,b], UndirectedEdge[b,c], UndirectedEdge[a,c]]
} // DeleteDuplicates) /. {UndirectedEdge -> List}
},
  {vertices, edges}
] ]

Наверняка это можно сделать проще, автор не претендует на хорошие знания вычислительной геометрии и графов в Wolfram. MeshCells возвращает линии и полигоны, из которых потом правилами замен вытаскиваются рёбра

Применим это на нашего dude

dudeMesh = dude // getMesh;
Graphics[GraphicsComplex[%[[1]], {Line /@ %[[2]]}]]

Рендеринг

У традиционного подхода к анимации с помощью Graphics есть проблема при работе с большим количеством линий, так как рендерится оно с помощью векторной графики (SVG). Поэтому для отрисовки воспользуемся Canvas API

Needs["Canvas2D`"->"ctx`"] // Quiet;

Вспомогательная функция для отрисовки всех стержней и шариков

drawBonds[context_, vert_, edges_] := Module[{},

  ctx`ClearRect[context, {0,0}, {500,500}];

  ctx`BeginPath[context];
  ctx`SetLineWidth[context, 4];
  ctx`SetStrokeStyle[context, RGBColor[0.4, 0.6, 0.9]];

  Do[
    ctx`MoveTo[context, {0, 500} - {-5, 5} vert[[p[[1]]]]];
    ctx`LineTo[context, {0, 500} - {-5, 5} vert[[p[[2]]]]];
  , {p, edges}];

  ctx`Stroke[context];

  ctx`SetFillStyle[context, RGBColor[0.4, 0.9, 0.6]];

  (
    ctx`BeginPath[context];
    ctx`Arc[context, {0, 500} - {-5, 5} #, 6, 0, 2.0 Pi];
    ctx`Fill[context];
  ) &/@ vert;

  ctx`Dispatch[context];
];

Симуляция

Переходим к самому интересному — симуляции

generatorSimulation[finalVerts_, finalBonds_] := Module[{
  coords, coords3, coords2
},

coords = finalVerts;
coords2 = finalVerts;
coords3 = finalVerts;

Function[Null, Do[
    (* обычный метод Верле *)
  
    coords3 = coords2;
    coords2 = coords;
    
    Module[{
      new = 2 coords2 - coords3 + Table[{0,-1}, Length[coords]] 0.0001
    },
           
      (* применяем ограничения на связи *)
    
      MapThread[Function[{i,j,l}, With[{d = new[[i]]-new[[j]]}, With[{m = Min[(l/(Norm[d]+0.001)) - 1), 0.1]},
        new[[i]] +=  0.5 m d;
        new[[j]] -=  0.5 m d;
      ] ]], RandomSample[finalBonds]// Transpose]; 
      
      new = Map[Function[xy,
        If[xy[[2]] < 0.0, {xy[[1]], 0.0}, xy]
      ], new];
      
      coords = new;
    ];

    , {2 5}];

    coords
]]

{finalVerts, finalPairs} = dudeMesh;

finalBonds = {
    finalPairs[[All,1]], (* первые узлы *)
    finalPairs[[All,2]], (* вторые узлы *)
    (* начальные длины *)
    Norm[finalVerts[[#[[1]]]] - finalVerts[[#[[2]]]]] &/@ finalPairs
} // Transpose;

With[{
  context2D = ctx`Canvas2D[],
  frame = CreateUUID[],
  calc = generatorSimulation[{#[[1]], #[[2]] + 50.0} &/@ finalVerts, finalBonds]
},
     
  EventHandler[frame, Function[Null,
    drawBonds[context2D, calc[], finalPairs]
  ]];
     
  Image[context2D, ImageResolution->{500,500}, Epilog->{
    AnimationFrameListener[context2D, "Event"->frame]
  }]
]

Задача многих тел

Хорошо известная и широко распространённая проблема в вычислительной физике (и не только там) — это работа с системами, где множество объектов связаны попарной функцией.

Аналогично алгоритму со шпицрутенами, столкновения шаров или дисков можно разрешать простым смещением их позиций. Такой метод нельзя назвать полностью физически корректным, но он может приближённо описывать ситуацию, когда частицы сталкиваются неупруго, теряя импульс и энергию в тепло:

К сожалению, это также означает, что сложность нашего алгоритма растёт как
N², где N — количество частиц. В результате неизбежно возникают проблемы с производительностью. Существует множество способов справиться с этим, но самый простой — скомпилировать нашу небольшую подпрограмму.

К счастью, в Wolfram Language есть встроенное выражение Compile, которое позволяет напрямую транслировать код либо в C, либо во внутренний байткод. Разумеется, это накладывает свои ограничения — например, на использование строго типизированных структур данных и ряд других особенностей

compiled = Compile[{
  {p, _Real, 3}, {c, _Real, 2}, {target, _Real, 1}
}, Module[{b = p, \[Delta]t = 0.005, circ = c},
    Do[
      
      (* сглаживаем движение обода *)
      
      circ[[1]] = circ[[1]] + \[Delta]t (target - circ[[1]]);
      
      b[[{1,2}]] = MapThread[If[Norm[#1-circ[[2]]] > 1.0, Module[{
        u = {#1, #2},
        vel = (circ[[1]] - circ[[2]])
      },

        
        u[[1]] = u[[1]] - vel;
        
        
        With[{n = -Normalize[u[[1]]-circ[[2]]], delta = u[[1]]-u[[2]]},
          u[[2]] = u[[1]];
          u[[1]] = u[[1]] -  0.9 n (n.delta);
        ];

        u[[1]] = u[[1]] + vel;

        If[Norm[u[[1]]-circ[[2]]] > 1.0, (*clamp to the region *)
          u[[1]] += (1.0 - Norm[u[[1]]-circ[[2]]]) Normalize[u[[1]]-circ[[2]]];
        ];

        u
      ], {#1, #2}]&, {b[[1]], b[[2]]}] // Transpose;

      b[[3]] = b[[2]];
      b[[2]] = b[[1]];
      b[[1]] = 2 b[[2]] - b[[3]] + Table[{0,-1}, {Length[b[[1]]]}] \[Delta]t \[Delta]t;


    (* Столкновения шариков *)
      
    Do[
      Do[
        If[i < j,
          Module[{pi = b[[1, i]], pj = b[[1, j]], d, dist, n, overlap},
            d = pj - pi;
            dist = Norm[d];
            If[dist < 0.05, (* явный радиус *)
              n = Normalize[d];
              overlap = 0.05 - dist;
              b[[1, i]] -= 0.5 overlap n;
              b[[1, j]] += 0.5 overlap n;
            ];
          ]
        ],
        {j, Length[b[[1]]]}
      ],
      {i, Length[b[[1]]]}
    ];   


      circ[[2]] = circ[[1]];
      
    , {10}];

    (* записываем новые положения обода в конец *)
    b[[-1, 1]] = circ[[1]];
    b[[-1, 2]] = circ[[2]];

    b
  ], CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C","RuntimeOptions" -> "Speed"];

Остальное зависит уже от нас — нужно лишь подключить нашу функцию обратно в основной цикл:

With[{
  p = Unique["ppp"],
  toDraw = Unique["ppp"],
  circ = Unique["ppp"],
  target = Unique["ppp"],
  frame = CreateUUID[]
},
  p = Table[Exp[-i/80.0]{Sin[i], Cos[i]}, {3}, {i, 100}];
  toDraw = p[[1]];
  circ = {{0,0}, {0,0}};
  target = {0,0};

  EventHandler[frame, Function[Null,
      With[{
        c = compiled[Join[p, {p[[3]]}], circ, target]
      },
        (* распаковка данных *)   
        circ = c[[4, {1,2}]];
        p = c[[1;;3]];
        toDraw = p[[1]];
      ]
  ]];

  {
    Graphics[{
      PointSize[0.03],
      Circle[circ[[1]] // Offload, 1.0], 
      Pink, Point[toDraw // Offload],
      AnimationFrameListener[p//Offload, "Event"->frame]
    }, 
      PlotRange->{1.4{-1,1}, 1.4{-1,1}}, 
      AspectRatio->1, ImageSize->Medium, 
      TransitionType->None
    ],

    EventHandler[InputJoystick[], Function[xy, target = xy;]]
  } // Row
]

Как видите, с этим набором инструментов можно продвинуться довольно далеко.

Бонус: подключаем Nintendo JoyCon

Очевидно, что задача легко обобщается на 3 измерения. Можно пойти еще дальше, подключить Nintendo JoyCon как в этой статье и взять за положение обода (сферы) вектор ускорения контроллера

Бонус: поиграем в бильярд

Выкинем обод, отключим гравитацию и приведем код в порядок

compiled = Compile[{
  {p, _Real, 3}
}, Module[{b = p},
    Do[

      b[[3]] = b[[2]];
      b[[2]] = b[[1]];
      b[[1]] = 2 b[[2]] - b[[3]];

    (* Particle-particle collisions *)
    Do[
      Do[
        If[i < j,
          Module[{pi = b[[1, i]], pj = b[[1, j]], d, dist, n, overlap},
            d = pj - pi;
            dist = Norm[d];
            If[dist < 0.1,
              n = Normalize[d];
              overlap = 0.1 - dist;
              b[[1, i]] -= 0.5 overlap n;
              b[[1, j]] += 0.5 overlap n;
            ];
          ]
        ],
        {j, Length[b[[1]]]}
      ],
      {i, Length[b[[1]]]}
    ];   
    , {3}];

    b
  ], CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C","RuntimeOptions" -> "Speed"];

Процедурно расставим наши шары в правильный треугольник

p = With[{s = 
  Join @@ Table[Table[{1-y,x}, {x, -1.0 + y,1.0 - y, 0.1}], {y, 0,1.0, 0.1}]
}, {s,s,s}];

Кием будет выступать дополнительный шар, которому мы заранее зададим скорость и ускорение

p[[1, -1]] = {-2,0};
p[[2, -1]] = {-2.05,0};
p[[3, -1]] = {-2,0};

Собираем всё вместе

p = Module[{s = 
    Join @@ Table[Table[{1-y,x}, {x,-1.0 + y, 1.0 - y,0.1}], {y, 0,1.0,0.1}]
},
    s = Append[s, {0,0}];
    {s,s,s}
];

p[[1, -1]] = {-2,0};
p[[2, -1]] = {-2.05,0};
p[[3, -1]] = {-2,0};

pDraw = p[[1]];

EventHandler["frame", Function[Null,
    p = compiled[p];
    pDraw = p[[1]];
]];

Graphics[{
    PointSize[0.03],
    Pink, Point[pDraw // Offload],
    AnimationFrameListener[pDraw//Offload, "Event"->"frame"]
  }
, PlotRange->3{{-1,1}, {-1,1}}, 
  AspectRatio->1, ImageSize->Medium, 
  TransitionType->None
]

Видимо я что-то не понимаю в бильярде ?

Ах да, шары... Поправим пару строчек в начале

p = Module[{s = 
  With[{r=0.04}, {d=2r}, Flatten[
   Table[
     {d i, (j - i/2) d Sqrt[3]},   
     {i, 0, 8},                    
     {j, 0, i}                    
   ],
   1
  ]]
}, 
    s = Append[s, {0,0}];
    {s,s,s}
];

Так то лучше ?